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

Calculate for Mu Px(t) by a superposition of the 4 Mu frequencies

Browse Source
This commit is contained in:
nemu
2016-01-14 11:55:43 +01:00
parent 18c90fef99
commit 8b2f0a16af
2 changed files with 97 additions and 76 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@ -14,9 +14,15 @@
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. Up to 3 Mu0 states are considered at the moment processes by a Monte-Carlo method. Consider transverse field geometry, and assume
(analogous to MuBC in Si, B||(100)), a non-precessing signal, and two precessing initial muon spin direction in x, and field applied along z. For PxMu(t) in
states ("nu_12" and "nu_34"). 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
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
@ -33,15 +39,14 @@
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 "3 Mu states" simulate muon spin phase under charge-exchange with "4 Mu transitions"
(as MuBC in Si, B || (100), for example): a non-precessing,
and two precessing signals (nu_12, nu_34).
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, also determine which Mu0 state has been formed 3) Determine next ionization time t_i; calculate Px(t_i) in Muonium; calculate the
at the capture event. Calculate muon spin precession. muon spin phase by acos(Px(t_i)).
4) get the next electron capture time, continue until t_d is reached. 4) get the next electron capture time, continue until t_d is reached; accumulate muon spin
phase.
***************************************************************************/ ***************************************************************************/
@ -95,6 +100,7 @@ PSimulateMuTransition::PSimulateMuTransition(UInt_t seed)
fMuPrecFreq12 = 0.; // Mu precession frequency of a 12 transition fMuPrecFreq12 = 0.; // Mu precession frequency of a 12 transition
fMuPrecFreq23 = 0.; // Mu precession frequency of a 23 transition fMuPrecFreq23 = 0.; // Mu precession frequency of a 23 transition
fMuPrecFreq14 = 0.; // Mu precession frequency of a 14 transition fMuPrecFreq14 = 0.; // Mu precession frequency of a 14 transition
fMuonPrecFreq = 0.; // muon precession frequency
fBfield = 0.01; // magnetic field (T) fBfield = 0.01; // magnetic field (T)
fCaptureRate = 0.01; // Mu+ capture rate (MHz) fCaptureRate = 0.01; // Mu+ capture rate (MHz)
fIonizationRate = 10.; // Mu0 ionization rate (MHz) fIonizationRate = 10.; // Mu0 ionization rate (MHz)
@ -103,8 +109,8 @@ PSimulateMuTransition::PSimulateMuTransition(UInt_t seed)
fMuonDecayTime = 0.; fMuonDecayTime = 0.;
fAsymmetry = 0.27; fAsymmetry = 0.27;
fMuFraction = 0.; fMuFraction = 0.;
fMuFractionState1 = 0.; fMuFractionState12 = 0.;
fMuFractionState2 = 0.; fMuFractionState23 = 0.;
fDebugFlag = kFALSE; fDebugFlag = kFALSE;
} }
@ -139,8 +145,8 @@ void PSimulateMuTransition::PrintSettings() const
cout << endl << "Mu ionizatioan rate (MHz) = " << fIonizationRate; cout << endl << "Mu ionizatioan rate (MHz) = " << fIonizationRate;
cout << endl << "Decay asymmetry = " << fAsymmetry; cout << endl << "Decay asymmetry = " << fAsymmetry;
cout << endl << "Muonium fraction = " << fMuFraction; cout << endl << "Muonium fraction = " << fMuFraction;
cout << endl << "Muonium fraction state1 = " << fMuFractionState1; cout << endl << "Muonium fraction state12 = " << fMuFractionState12;
cout << endl << "Muonium fraction state2 = " << fMuFractionState2; cout << endl << "Muonium fraction state23 = " << fMuFractionState23;
cout << endl << "Number of particles to simulate = " << fNmuons; cout << endl << "Number of particles to simulate = " << fNmuons;
cout << endl << "Initial muon spin phase (degree) = " << fInitialPhase; cout << endl << "Initial muon spin phase (degree) = " << fInitialPhase;
cout << endl << "Debug flag = " << fDebugFlag; cout << endl << "Debug flag = " << fDebugFlag;
@ -224,38 +230,53 @@ Double_t PSimulateMuTransition::NextEventTime(const Double_t &EventRate)
* \param time duration of precession (us); * \param time duration of precession (us);
* \param frequency muon spin precession frequency (MHz); * \param frequency muon spin precession frequency (MHz);
*/ */
Double_t PSimulateMuTransition::PrecessionPhase(const Double_t &time, const Double_t &frequency) Double_t PSimulateMuTransition::PrecessionPhase(const Double_t &time, const TString chargeState)
{ {
return TMath::TwoPi()*frequency*time; Double_t muonPhaseX;
Double_t muoniumPolX = 0;
if (chargeState == "Mu+")
muonPhaseX = TMath::TwoPi()*fMuonPrecFreq*time;
else if (chargeState == "Mu0"){
muoniumPolX = 0.5 *
(fMuFractionState12 * (TMath::Cos(TMath::TwoPi()*fMuPrecFreq12*time) + TMath::Cos(TMath::TwoPi()*fMuPrecFreq34*time)) +
fMuFractionState23 * (TMath::Cos(TMath::TwoPi()*fMuPrecFreq23*time) + TMath::Cos(TMath::TwoPi()*fMuPrecFreq14*time)));
muonPhaseX = TMath::ACos(muoniumPolX);
}
else
muonPhaseX = 0.;
return muonPhaseX;
} }
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------
// Event (private) // Event (private)
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------
/** /**
* <p> Generates "muon event": simulate muon spin phase under charge-exchange * <p> Generates "muon event": simulate muon spin phase under charge-exchange with
* with "3 Mu states" (as MuBC in Si, B || (100), for example): a non-precessing, * a neutral muonium state in transverse field, where the polarization evolution
* and two precessing signals (nu_12, nu_34). * 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 * 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, also determine which Mu0 state has been formed * 3) Determine next ionization time t_i; calculate Px(t_i) in Muonium; calculate the
* at the capture event. Calculate muon spin precession. * muon spin phase by acos(Px(t_i)).
* 4) get the next electron capture time, continue until t_d is reached. * 4) get the next electron capture time, continue until t_d is reached.
* *
* <p> For isotropic muonium, TF: * <p> For isotropic muonium, TF:
* nu_12 and nu_34 with equal probabilities, probability for both states fMuFractionState1 * 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 fMuFractionState2 * ni_23 and nu_14 with equal probabilities, probability for both states fMuFractionState23
* *
* \param muonString if eq. "Mu+" begin with Mu+ precession * \param muonString if eq. "Mu+" begin with Mu+ precession
*/ */
void PSimulateMuTransition::Event(const TString muonString) void PSimulateMuTransition::Event(const TString muonString)
{ {
Double_t eventTime, eventDiffTime, captureTime, ionizationTime; Double_t eventTime, eventDiffTime, captureTime, ionizationTime;
Double_t muonPrecessionFreq, muoniumPrecessionFreq; // MHz // Double_t muonPrecessionFreq, muoniumPrecessionFreq; // MHz
Double_t rndm, frac1, frac2; // Double_t rndm, frac1, frac2;
muonPrecessionFreq = fMuonGyroRatio * fBfield; fMuonPrecFreq = fMuonGyroRatio * fBfield;
// charge-exchange loop until muon decay // charge-exchange loop until muon decay
eventTime = 0.; eventTime = 0.;
@ -270,10 +291,10 @@ void PSimulateMuTransition::Event(const TString muonString)
eventTime += captureTime; eventTime += captureTime;
if (fDebugFlag) cout << "Capture time = " << captureTime << " Phase = " << fMuonPhase << endl; if (fDebugFlag) cout << "Capture time = " << captureTime << " Phase = " << fMuonPhase << endl;
if (eventTime < fMuonDecayTime) if (eventTime < fMuonDecayTime)
fMuonPhase += PrecessionPhase(captureTime, muonPrecessionFreq); fMuonPhase += PrecessionPhase(captureTime, "Mu+");
else{ //muon decays; handle precession prior to muon decay else{ //muon decays; handle precession prior to muon decay
eventDiffTime = fMuonDecayTime - (eventTime - captureTime); eventDiffTime = fMuonDecayTime - (eventTime - captureTime);
fMuonPhase += PrecessionPhase(eventDiffTime, muonPrecessionFreq); fMuonPhase += PrecessionPhase(eventDiffTime, "Mu+");
break; break;
} }
@ -281,31 +302,30 @@ void PSimulateMuTransition::Event(const TString muonString)
ionizationTime = NextEventTime(fIonizationRate); ionizationTime = NextEventTime(fIonizationRate);
eventTime += ionizationTime; eventTime += ionizationTime;
// determine Mu state // determine Mu state
rndm = fRandom->Rndm(); // rndm = fRandom->Rndm();
frac1 = 1. - fMuFractionState1 - fMuFractionState2; // non-precessing Mu states // frac1 = 1. - fMuFractionState1 - fMuFractionState2; // non-precessing Mu states
frac2 = 1. - fMuFractionState2; // frac2 = 1. - fMuFractionState2;
if ( rndm < frac1 ) // if ( rndm < frac1 )
muoniumPrecessionFreq = 0.; // muoniumPrecessionFreq = 0.;
else if (rndm >= frac1 && rndm <= frac2){ // else if (rndm >= frac1 && rndm <= frac2){
if (fRandom->Rndm() <= 0.5) // if (fRandom->Rndm() <= 0.5)
muoniumPrecessionFreq = fMuPrecFreq12; // muoniumPrecessionFreq = fMuPrecFreq12;
else // else
muoniumPrecessionFreq = fMuPrecFreq34; // muoniumPrecessionFreq = fMuPrecFreq34;
} // }
else{ // else{
if (fRandom->Rndm() <= 0.5) // if (fRandom->Rndm() <= 0.5)
muoniumPrecessionFreq = fMuPrecFreq23; // muoniumPrecessionFreq = fMuPrecFreq23;
else // else
muoniumPrecessionFreq = fMuPrecFreq14; // muoniumPrecessionFreq = fMuPrecFreq14;
} // }
if (fDebugFlag) cout << "Ioniza. time = " << ionizationTime << " Freq = " << muoniumPrecessionFreq if (fDebugFlag) cout << "Ioniza. time = " << ionizationTime << " Phase = " << fMuonPhase << endl;
<< " Phase = " << fMuonPhase << endl;
if (eventTime < fMuonDecayTime) if (eventTime < fMuonDecayTime)
fMuonPhase += PrecessionPhase(ionizationTime, muoniumPrecessionFreq); fMuonPhase += PrecessionPhase(ionizationTime, "Mu0");
else{ //muon decays; handle precession prior to muon decay else{ //muon decays; handle precession prior to muon decay
eventDiffTime = fMuonDecayTime - (eventTime - ionizationTime); eventDiffTime = fMuonDecayTime - (eventTime - ionizationTime);
fMuonPhase += PrecessionPhase(eventDiffTime, muoniumPrecessionFreq); fMuonPhase += PrecessionPhase(eventDiffTime, "Mu0");
break; break;
} }
} }
@ -314,32 +334,31 @@ void PSimulateMuTransition::Event(const TString muonString)
ionizationTime = NextEventTime(fIonizationRate); ionizationTime = NextEventTime(fIonizationRate);
eventTime += ionizationTime; eventTime += ionizationTime;
// determine Mu state // determine Mu state
rndm = fRandom->Rndm(); // rndm = fRandom->Rndm();
frac1 = 1. - fMuFractionState1 - fMuFractionState2; // non-precessing Mu states // frac1 = 1. - fMuFractionState1 - fMuFractionState2; // non-precessing Mu states
frac2 = 1. - fMuFractionState2; // frac2 = 1. - fMuFractionState2;
if ( rndm < frac1 ) // if ( rndm < frac1 )
muoniumPrecessionFreq = 0.; // muoniumPrecessionFreq = 0.;
else if (rndm >= frac1 && rndm <= frac2){ // else if (rndm >= frac1 && rndm <= frac2){
if (fRandom->Rndm() <= 0.5) // if (fRandom->Rndm() <= 0.5)
muoniumPrecessionFreq = fMuPrecFreq12; // muoniumPrecessionFreq = fMuPrecFreq12;
else // else
muoniumPrecessionFreq = fMuPrecFreq34; // muoniumPrecessionFreq = fMuPrecFreq34;
} // }
else{ // else{
if (fRandom->Rndm() <= 0.5) // if (fRandom->Rndm() <= 0.5)
muoniumPrecessionFreq = fMuPrecFreq23; // muoniumPrecessionFreq = fMuPrecFreq23;
else // else
muoniumPrecessionFreq = fMuPrecFreq14; // muoniumPrecessionFreq = fMuPrecFreq14;
} // }
if (fDebugFlag) if (fDebugFlag)
cout << "Mu Ioniza. time = " << ionizationTime << " Freq = " << muoniumPrecessionFreq cout << "Mu Ioniza. time = " << ionizationTime << " Phase = " << fMuonPhase << endl;
<< " Phase = " << fMuonPhase << endl;
if (eventTime < fMuonDecayTime) if (eventTime < fMuonDecayTime)
fMuonPhase += PrecessionPhase(ionizationTime, muoniumPrecessionFreq); fMuonPhase += PrecessionPhase(ionizationTime, "Mu0");
else{ //muon decays; handle precession prior to muon decay else{ //muon decays; handle precession prior to muon decay
eventDiffTime = fMuonDecayTime - (eventTime - ionizationTime); eventDiffTime = fMuonDecayTime - (eventTime - ionizationTime);
fMuonPhase += PrecessionPhase(eventDiffTime, muoniumPrecessionFreq); fMuonPhase += PrecessionPhase(eventDiffTime, "Mu0");
break; break;
} }
@ -348,10 +367,10 @@ void PSimulateMuTransition::Event(const TString muonString)
eventTime += captureTime; eventTime += captureTime;
if (fDebugFlag) cout << "Capture time = " << captureTime << " Phase = " << fMuonPhase << endl; if (fDebugFlag) cout << "Capture time = " << captureTime << " Phase = " << fMuonPhase << endl;
if (eventTime < fMuonDecayTime) if (eventTime < fMuonDecayTime)
fMuonPhase += PrecessionPhase(captureTime, muonPrecessionFreq); fMuonPhase += PrecessionPhase(captureTime, "Mu+");
else{ //muon decays; handle precession prior to muon decay else{ //muon decays; handle precession prior to muon decay
eventDiffTime = fMuonDecayTime - (eventTime - captureTime); eventDiffTime = fMuonDecayTime - (eventTime - captureTime);
fMuonPhase += PrecessionPhase(eventDiffTime, muonPrecessionFreq); fMuonPhase += PrecessionPhase(eventDiffTime, "Mu+");
break; break;
} }
} }

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

@ -57,8 +57,8 @@ class PSimulateMuTransition : public TObject
virtual void SetIonizationRate(Double_t value){ fIonizationRate = value; } //!< sets Mu0 ionization rate (MHz) virtual void SetIonizationRate(Double_t value){ fIonizationRate = value; } //!< sets Mu0 ionization rate (MHz)
virtual void SetDecayAsymmetry(Double_t value){ fAsymmetry = value; } //!< muon decay asymmetry virtual void SetDecayAsymmetry(Double_t value){ fAsymmetry = value; } //!< muon decay asymmetry
virtual void SetMuFraction(Double_t value){ fMuFraction = value; } //!< Muonium fraction virtual void SetMuFraction(Double_t value){ fMuFraction = value; } //!< Muonium fraction
virtual void SetMuFractionState1(Double_t value){ fMuFractionState1 = value; } virtual void SetMuFractionState12(Double_t value){ fMuFractionState12 = value; }
virtual void SetMuFractionState2(Double_t value){ fMuFractionState2 = value; } virtual void SetMuFractionState23(Double_t value){ fMuFractionState23 = value; }
virtual Bool_t IsValid() { return fValid; } virtual Bool_t IsValid() { return fValid; }
virtual void SetSeed(UInt_t seed); virtual void SetSeed(UInt_t seed);
@ -77,6 +77,7 @@ class PSimulateMuTransition : public TObject
Double_t fMuPrecFreq34; //!< Mu transition frequency 34 (MHz) Double_t fMuPrecFreq34; //!< Mu transition frequency 34 (MHz)
Double_t fMuPrecFreq23; //!< Mu transition frequency 23 (MHz) Double_t fMuPrecFreq23; //!< Mu transition frequency 23 (MHz)
Double_t fMuPrecFreq14; //!< Mu transition frequency 14 (MHz) Double_t fMuPrecFreq14; //!< Mu transition frequency 14 (MHz)
Double_t fMuonPrecFreq; //!< muon precession frequency (MHz)
Double_t fCaptureRate; //!< Mu+ electron capture rate (MHz) Double_t fCaptureRate; //!< Mu+ electron capture rate (MHz)
Double_t fIonizationRate; //!< Mu0 ionization rate (MHz) Double_t fIonizationRate; //!< Mu0 ionization rate (MHz)
Double_t fInitialPhase; //!< initial muon spin phase Double_t fInitialPhase; //!< initial muon spin phase
@ -84,13 +85,14 @@ class PSimulateMuTransition : public TObject
Double_t fMuonPhase; //!< phase of muon spin Double_t fMuonPhase; //!< phase of muon spin
Double_t fAsymmetry; //!< muon decay asymmetry Double_t fAsymmetry; //!< muon decay asymmetry
Double_t fMuFraction; //!< total Mu fraction [0,1] Double_t fMuFraction; //!< total Mu fraction [0,1]
Double_t fMuFractionState1; //!< fraction of Mu in state 1 Double_t fMuFractionState12; //!< fraction of Mu in state 12, 34
Double_t fMuFractionState2; //!< fraction of Mu in state 2 Double_t fMuFractionState23; //!< fraction of Mu in state 23, 14
Int_t fNmuons; //!< number of muons to simulate Int_t fNmuons; //!< number of muons to simulate
Bool_t fDebugFlag; //!< debug flag Bool_t fDebugFlag; //!< debug flag
virtual Double_t NextEventTime(const Double_t &EventRate); 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 Double_t &frequency);
virtual Double_t PrecessionPhase(const Double_t &time, const TString chargeState);
virtual void Event(const TString muonString); virtual void Event(const TString muonString);
ClassDef(PSimulateMuTransition, 0) ClassDef(PSimulateMuTransition, 0)