From 908f728744b956b7a186a1acd631eb14c2647509 Mon Sep 17 00:00:00 2001 From: Thomas Prokscha Date: Thu, 14 Apr 2016 15:28:33 +0200 Subject: [PATCH 1/2] Changed calculation of charge-exchange processes to a product of complex polarization functions, as in the paper of M. Senba --- .../MuTransition/PSimulateMuTransition.cpp | 250 ++++++++---------- .../MuTransition/PSimulateMuTransition.h | 7 +- src/tests/MuTransition/runMuSimulation.C | 3 +- 3 files changed, 120 insertions(+), 140 deletions(-) diff --git a/src/tests/MuTransition/PSimulateMuTransition.cpp b/src/tests/MuTransition/PSimulateMuTransition.cpp index a983b3a5..e5b00c34 100644 --- a/src/tests/MuTransition/PSimulateMuTransition.cpp +++ b/src/tests/MuTransition/PSimulateMuTransition.cpp @@ -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). ***************************************************************************/ @@ -138,15 +141,21 @@ 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; @@ -181,23 +190,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 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 +245,29 @@ Double_t PSimulateMuTransition::NextEventTime(const Double_t &EventRate) //-------------------------------------------------------------------------- // Phase (private) //-------------------------------------------------------------------------- -/** - *

Determines phase of the muon spin +// /** +/* *

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 +277,26 @@ 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 = + 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)); + } 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 +318,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 +340,100 @@ Double_t PSimulateMuTransition::GTSpinFlip(const Double_t &time) // Event (private) //-------------------------------------------------------------------------- /** - *

Generates "muon event": simulate muon spin phase under charge-exchange with + *

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. * *

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 * + *

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; } diff --git a/src/tests/MuTransition/PSimulateMuTransition.h b/src/tests/MuTransition/PSimulateMuTransition.h index 1a8d0f31..a5a44f4f 100644 --- a/src/tests/MuTransition/PSimulateMuTransition.h +++ b/src/tests/MuTransition/PSimulateMuTransition.h @@ -94,11 +94,10 @@ class PSimulateMuTransition : public TObject 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) }; diff --git a/src/tests/MuTransition/runMuSimulation.C b/src/tests/MuTransition/runMuSimulation.C index 07d66475..c88074b0 100644 --- a/src/tests/MuTransition/runMuSimulation.C +++ b/src/tests/MuTransition/runMuSimulation.C @@ -69,6 +69,7 @@ void runMuSimulation() Double_t MuFrac23 = 2*0.013; //Mu in states 23 and 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; @@ -97,7 +98,7 @@ void runMuSimulation() 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"); From b65e07b8c8b8a38de096200756bb9f7bb81f4b5e Mon Sep 17 00:00:00 2001 From: nemu Date: Fri, 15 Apr 2016 19:06:10 +0200 Subject: [PATCH 2/2] added fractions for 34 and 14 lines for anisotropic Mu --- .../MuTransition/PSimulateMuTransition.cpp | 18 ++++++++++-------- src/tests/MuTransition/PSimulateMuTransition.h | 8 ++++++-- src/tests/MuTransition/runMuSimulation.C | 14 ++++++++++---- 3 files changed, 26 insertions(+), 14 deletions(-) diff --git a/src/tests/MuTransition/PSimulateMuTransition.cpp b/src/tests/MuTransition/PSimulateMuTransition.cpp index e5b00c34..d84dfe5d 100644 --- a/src/tests/MuTransition/PSimulateMuTransition.cpp +++ b/src/tests/MuTransition/PSimulateMuTransition.cpp @@ -114,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; } @@ -159,7 +161,9 @@ void PSimulateMuTransition::PrintSettings() const 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; @@ -287,13 +291,11 @@ TComplex PSimulateMuTransition::GTFunction(const Double_t &time, const TString c complexPol = TComplex::Exp(-TComplex::I()*twoPi*fMuonPrecFreq*time); else{ complexPol = - 0.5 * fMuFractionState12 * - (TComplex::Exp(TComplex::I()*twoPi*fMuPrecFreq12*time) + - TComplex::Exp(-TComplex::I()*twoPi*fMuPrecFreq34*time)) + (fMuFractionState12 * TComplex::Exp(TComplex::I()*twoPi*fMuPrecFreq12*time) + + fMuFractionState34 * TComplex::Exp(-TComplex::I()*twoPi*fMuPrecFreq34*time)) + - 0.5 * fMuFractionState23 * - (TComplex::Exp(TComplex::I()*twoPi*fMuPrecFreq23*time) + - TComplex::Exp(TComplex::I()*twoPi*fMuPrecFreq14*time)); + (fMuFractionState23 * TComplex::Exp(TComplex::I()*twoPi*fMuPrecFreq23*time) + + fMuFractionState14 * TComplex::Exp(TComplex::I()*twoPi*fMuPrecFreq14*time)); } return complexPol; diff --git a/src/tests/MuTransition/PSimulateMuTransition.h b/src/tests/MuTransition/PSimulateMuTransition.h index a5a44f4f..552ff132 100644 --- a/src/tests/MuTransition/PSimulateMuTransition.h +++ b/src/tests/MuTransition/PSimulateMuTransition.h @@ -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,8 +90,10 @@ 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 diff --git a/src/tests/MuTransition/runMuSimulation.C b/src/tests/MuTransition/runMuSimulation.C index c88074b0..80e0e589 100644 --- a/src/tests/MuTransition/runMuSimulation.C +++ b/src/tests/MuTransition/runMuSimulation.C @@ -65,8 +65,10 @@ 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 @@ -90,8 +92,10 @@ 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 @@ -169,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);