Modified for isotropic Mu in TF.

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"
+  1) Precession frequencies of "nu_12", "nu_34", "nu_23", "nu_14"
-  2) fractions of nu_12, nu_34
+  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
+  fMuPrecFreq34     = 4463.; // vacuum Mu hyperfine coupling constant
-  fMuPrecFreq2      = 0.;    // Mu precession frequency of a 2nd Mu transition
+  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 12 (MHz)  = " << fMuPrecFreq12;
-  cout << endl << "Mu precession frequency state2 (MHz)  = " << fMuPrecFreq2;
+  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.
  *
+ * <p> 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)
+     else if (rndm >= frac1 && rndm <= frac2){
-       muoniumPrecessionFreq = fMuPrecFreq1;
+       if (fRandom->Rndm() <= 0.5) 
-     else
+	muoniumPrecessionFreq = fMuPrecFreq12;
-       muoniumPrecessionFreq = fMuPrecFreq2;
+       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)
+     else if (rndm >= frac1 && rndm <= frac2){
-       muoniumPrecessionFreq = fMuPrecFreq1;
+       if (fRandom->Rndm() <= 0.5) 
-     else
+	muoniumPrecessionFreq = fMuPrecFreq12;
-       muoniumPrecessionFreq = fMuPrecFreq2;
+       else
+        muoniumPrecessionFreq = fMuPrecFreq34;
+     }
+     else{
+       if (fRandom->Rndm() <= 0.5)
+        muoniumPrecessionFreq = fMuPrecFreq23;
+       else
+        muoniumPrecessionFreq = fMuPrecFreq14;
+     }
 
      if (fDebugFlag) 
       cout << "Mu Ioniza. time = " << ionizationTime << " Freq = " << muoniumPrecessionFreq 

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 SetMuPrecFreq12(Double_t value) { fMuPrecFreq12 = value; } //!< sets Mu transition frequency (MHz)
-    virtual void SetMuPrecFreq2(Double_t value) { fMuPrecFreq2 = value; } //!< sets Mu hyperfine coupling (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  fMuPrecFreq12;    //!< Mu transition frequency 12 (MHz)
-    Double_t  fMuPrecFreq2;     //!< Mu precession frequency of state 2 (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

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
+  //prepare to run simulation; here: isotropic Mu in Germanium
-  simulateMuTransition->SetMuPrecFreq1(41.);   // MHz
+  Double_t ionRate;  //assume Arrhenius behaviour ionRate = preFac*exp(-EA/kT)
-  simulateMuTransition->SetMuPrecFreq2(-35.);  // MHz
+  Double_t capRate;  //assume that capture rate varies as sqrt(T), capRate = sigma*v*p , v ~ sqrt(T) 
-  simulateMuTransition->SetMuFraction(0.5);    // initial Mu fraction
+  Double_t EA, T;     //activation energy (meV) and temperature (K) 
-  simulateMuTransition->SetMuFractionState1(0.42);  // 100% of Mu in state 1
+  EA = 170.;
-  simulateMuTransition->SetMuFractionState2(0.32);  // 0% of Mu in state 2
+  T  = 220.;
-  simulateMuTransition->SetBfield(0.01);           // Tesla
+  ionRate = 2.9e7 * exp(-EA/(0.08625*T)); // Ge: 2.9*10^7MHz "attempt" frequency; 1K = 0.08625 meV 
-  simulateMuTransition->SetCaptureRate(1.5);       // MHz
+  capRate = 0.00001*sqrt(T/200.); 
-  simulateMuTransition->SetIonizationRate(250.);    // MHz
+
+  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