diff --git a/src/tests/MuTransition/PSimulateMuTransition.cpp b/src/tests/MuTransition/PSimulateMuTransition.cpp index 8432dda8..113e0982 100644 --- a/src/tests/MuTransition/PSimulateMuTransition.cpp +++ b/src/tests/MuTransition/PSimulateMuTransition.cpp @@ -18,8 +18,8 @@ (analogous to MuBC in Si, B||(100)), a non-precessing signal, and two precessing states ("nu_12" and "nu_34"). Parameters: - 1) Precession frequencies of "nu_12" and "nu_34" - 2) fractions of nu_12, nu_34 + 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 @@ -91,8 +91,10 @@ PSimulateMuTransition::PSimulateMuTransition(UInt_t seed) } fNmuons = 100; // number of muons to simulate - fMuPrecFreq1 = 4463.; // vacuum Mu hyperfine coupling constant - fMuPrecFreq2 = 0.; // Mu precession frequency of a 2nd Mu transition + fMuPrecFreq34 = 4463.; // vacuum Mu hyperfine coupling constant + fMuPrecFreq12 = 0.; // Mu precession frequency of a 12 transition + fMuPrecFreq23 = 0.; // Mu precession frequency of a 23 transition + fMuPrecFreq14 = 0.; // Mu precession frequency of a 14 transition fBfield = 0.01; // magnetic field (T) fCaptureRate = 0.01; // Mu+ capture rate (MHz) fIonizationRate = 10.; // Mu0 ionization rate (MHz) @@ -128,8 +130,10 @@ PSimulateMuTransition::~PSimulateMuTransition() */ void PSimulateMuTransition::PrintSettings() const { - cout << endl << "Mu precession frequency state1 (MHz) = " << fMuPrecFreq1; - cout << endl << "Mu precession frequency state2 (MHz) = " << fMuPrecFreq2; + 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 << "B field (T) = " << fBfield; cout << endl << "Mu+ electron capture rate (MHz) = " << fCaptureRate; cout << endl << "Mu ionizatioan rate (MHz) = " << fIonizationRate; @@ -239,6 +243,10 @@ Double_t PSimulateMuTransition::PrecessionPhase(const Double_t &time, const Doub * at the capture event. Calculate muon spin precession. * 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 fMuFractionState1 + * ni_23 and nu_14 with equal probabilities, probability for both states fMuFractionState2 + * * \param muonString if eq. "Mu+" begin with Mu+ precession */ void PSimulateMuTransition::Event(const TString muonString) @@ -278,10 +286,18 @@ void PSimulateMuTransition::Event(const TString muonString) frac2 = 1. - fMuFractionState2; if ( rndm < frac1 ) muoniumPrecessionFreq = 0.; - else if (rndm >= frac1 && rndm <= frac2) - muoniumPrecessionFreq = fMuPrecFreq1; - else - muoniumPrecessionFreq = fMuPrecFreq2; + 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 << " Freq = " << muoniumPrecessionFreq << " Phase = " << fMuonPhase << endl; @@ -303,10 +319,18 @@ void PSimulateMuTransition::Event(const TString muonString) frac2 = 1. - fMuFractionState2; if ( rndm < frac1 ) muoniumPrecessionFreq = 0.; - else if (rndm >= frac1 && rndm <= frac2) - muoniumPrecessionFreq = fMuPrecFreq1; - else - muoniumPrecessionFreq = fMuPrecFreq2; + 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 << " Freq = " << muoniumPrecessionFreq diff --git a/src/tests/MuTransition/PSimulateMuTransition.h b/src/tests/MuTransition/PSimulateMuTransition.h index 00716683..cbe5cc66 100644 --- a/src/tests/MuTransition/PSimulateMuTransition.h +++ b/src/tests/MuTransition/PSimulateMuTransition.h @@ -49,8 +49,10 @@ class PSimulateMuTransition : public TObject virtual void SetNmuons(Int_t value) { fNmuons = value; } //!< number of muons virtual void SetDebugFlag(Bool_t value) { fDebugFlag = value; } //!< debug flag virtual void SetBfield(Double_t value) { fBfield = value; } //!< sets magnetic field (T) - virtual void SetMuPrecFreq1(Double_t value) { fMuPrecFreq1 = value; } //!< sets Mu hyperfine coupling (MHz) - virtual void SetMuPrecFreq2(Double_t value) { fMuPrecFreq2 = value; } //!< sets Mu hyperfine coupling (MHz) + virtual void SetMuPrecFreq12(Double_t value) { fMuPrecFreq12 = value; } //!< sets Mu transition frequency (MHz) + virtual void SetMuPrecFreq34(Double_t value) { fMuPrecFreq34 = value; } //!< sets Mu transition frequency (MHz) + virtual void SetMuPrecFreq23(Double_t value) { fMuPrecFreq23 = value; } //!< sets Mu transition frequency (MHz) + virtual void SetMuPrecFreq14(Double_t value) { fMuPrecFreq14 = value; } //!< sets Mu transition frequency (MHz) virtual void SetCaptureRate(Double_t value){ fCaptureRate = value; } //!< sets Mu+ electron capture 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 @@ -71,8 +73,10 @@ class PSimulateMuTransition : public TObject TRandom2 *fRandom; Double_t fBfield; //!< magnetic field (T) - Double_t fMuPrecFreq1; //!< Mu precession frequency of state 1 (MHz) - Double_t fMuPrecFreq2; //!< Mu precession frequency of state 2 (MHz) + Double_t fMuPrecFreq12; //!< Mu transition frequency 12 (MHz) + Double_t fMuPrecFreq34; //!< Mu transition frequency 34 (MHz) + Double_t fMuPrecFreq23; //!< Mu transition frequency 23 (MHz) + Double_t fMuPrecFreq14; //!< Mu transition frequency 14 (MHz) Double_t fCaptureRate; //!< Mu+ electron capture rate (MHz) Double_t fIonizationRate; //!< Mu0 ionization rate (MHz) Double_t fInitialPhase; //!< initial muon spin phase diff --git a/src/tests/MuTransition/runMuSimulation.C b/src/tests/MuTransition/runMuSimulation.C index 4e2ae2ca..b60c576c 100644 --- a/src/tests/MuTransition/runMuSimulation.C +++ b/src/tests/MuTransition/runMuSimulation.C @@ -43,14 +43,15 @@ void runMuSimulation() decayAnaModule = histosFolder->AddFolder("DecayAnaModule", "muSR decay histograms"); // feed run info header - UInt_t runNo = 9102; + UInt_t runNo = 9702; TString tstr; runInfo = gROOT->GetRootFolder()->AddFolder("RunInfo", "LEM RunInfo"); gROOT->GetListOfBrowsables()->Add(runInfo, "RunInfo"); header = new TLemRunHeader(); tstr = TString("0"); tstr += runNo; - tstr += TString(" - test"); + tstr += TString("Ge, Mu-frac 1.0, Mu12 737MHz (0.44), Mu34 -1622MHz(0.44), T=220K/EA=170meV, Cap.(200K) 0.0MHz, 100mT"); + header->SetRunTitle(tstr.Data()); header->SetLemSetup("trivial"); header->SetRunNumber(runNo); @@ -60,7 +61,7 @@ void runMuSimulation() header->SetSampleHV(0.0, 0.01); header->SetImpEnergy(31.8); header->SetSampleTemperature(0.2, 0.001); - header->SetSampleBField(-1.0, 0.1); + header->SetSampleBField(100.0, 0.1); header->SetTimeResolution(1.); header->SetNChannels(12001); header->SetNHist(2); @@ -86,15 +87,25 @@ void runMuSimulation() return; } - //prepare to run simulation - simulateMuTransition->SetMuPrecFreq1(41.); // MHz - simulateMuTransition->SetMuPrecFreq2(-35.); // MHz - simulateMuTransition->SetMuFraction(0.5); // initial Mu fraction - simulateMuTransition->SetMuFractionState1(0.42); // 100% of Mu in state 1 - simulateMuTransition->SetMuFractionState2(0.32); // 0% of Mu in state 2 - simulateMuTransition->SetBfield(0.01); // Tesla - simulateMuTransition->SetCaptureRate(1.5); // MHz - simulateMuTransition->SetIonizationRate(250.); // MHz + //prepare to run simulation; here: isotropic Mu in Germanium + Double_t ionRate; //assume Arrhenius behaviour ionRate = preFac*exp(-EA/kT) + Double_t capRate; //assume that capture rate varies as sqrt(T), capRate = sigma*v*p , v ~ sqrt(T) + Double_t EA, T; //activation energy (meV) and temperature (K) + EA = 170.; + T = 220.; + ionRate = 2.9e7 * exp(-EA/(0.08625*T)); // Ge: 2.9*10^7MHz "attempt" frequency; 1K = 0.08625 meV + capRate = 0.00001*sqrt(T/200.); + + simulateMuTransition->SetMuPrecFreq12(737.3); // MHz + simulateMuTransition->SetMuPrecFreq34(-1622.2); // MHz + simulateMuTransition->SetMuPrecFreq23(2051.6); // MHz + simulateMuTransition->SetMuPrecFreq14(4111.2); // MHz + simulateMuTransition->SetMuFraction(1.0); // initial Mu fraction + simulateMuTransition->SetMuFractionState1(0.88); // Mu in states 12, 34 + simulateMuTransition->SetMuFractionState2(0.12); // Mu in states 23, 14 + simulateMuTransition->SetBfield(0.1); // Tesla + simulateMuTransition->SetCaptureRate(capRate); // MHz + simulateMuTransition->SetIonizationRate(ionRate); // MHz simulateMuTransition->SetNmuons(1e7); simulateMuTransition->SetDecayAsymmetry(0.27); simulateMuTransition->SetDebugFlag(kFALSE); // to print time and phase during charge-changing cycle