Added two precessing and one non-precessing Mu state.

src/tests/MuTransition/PSimulateMuTransition.cpp

@@ -55,7 +55,8 @@ PSimulateMuTransition::PSimulateMuTransition(UInt_t seed)
   }
 
   fNmuons           = 100;   // number of muons to simulate
-  fMuCoupling     = 4463.; // vacuum Mu hyperfine coupling constant
+  fMuPrecFreq1      = 4463.; // vacuum Mu hyperfine coupling constant
+  fMuPrecFreq2      = 0.;    // Mu precession frequency of a 2nd Mu transition
   fBfield           = 0.01;  // magnetic field (T)
   fCaptureRate      = 0.01;  // Mu+ capture rate (MHz)
   fIonizationRate   = 10.;   // Mu0 ionization rate (MHz)
@@ -64,6 +65,8 @@ PSimulateMuTransition::PSimulateMuTransition(UInt_t seed)
   fMuonDecayTime    = 0.;
   fAsymmetry        = 0.27;
   fMuFraction       = 0.;
+  fMuFractionState1 = 0.;
+  fMuFractionState2 = 0.;
   fDebugFlag        = kFALSE;
 }
 
@@ -89,12 +92,15 @@ PSimulateMuTransition::~PSimulateMuTransition()
  */
 void PSimulateMuTransition::PrintSettings() const
 {
-  cout << endl << "Mu hyperfine couling (MHz)            = " << fMuCoupling;
+  cout << endl << "Mu precession frequency state1 (MHz)  = " << fMuPrecFreq1;
+  cout << endl << "Mu precession frequency state2 (MHz)  = " << fMuPrecFreq2;
   cout << endl << "B field (T)                           = " << fBfield;
   cout << endl << "Mu+ electron capture rate (MHz)       = " << fCaptureRate;
   cout << endl << "Mu ionizatioan rate (MHz)             = " << fIonizationRate;
   cout << endl << "Decay asymmetry                       = " << fAsymmetry;
   cout << endl << "Muonium fraction                      = " << fMuFraction;
+  cout << endl << "Muonium fraction state1               = " << fMuFractionState1;
+  cout << endl << "Muonium fraction state2               = " << fMuFractionState2;
   cout << endl << "Number of particles to simulate       = " << fNmuons;
   cout << endl << "Initial muon spin phase (degree)      = " << fInitialPhase;
   cout << endl << "Debug flag                            = " << fDebugFlag;
@@ -134,9 +140,9 @@ void PSimulateMuTransition::Run(TH1F *histoForward, TH1F *histoBackward)
     fMuonDecayTime     = NextEventTime(fMuonDecayRate);
     // initial muon state Mu+ or Mu0?
     if (fRandom->Rndm() <= 1.-fMuFraction) 
-     MuonEvent();
+     Event("Mu+");
     else
-     MuoniumEvent();
+     Event("");
 
     // fill 50% in "forward", and 50% in "backward" detector to get independent
     // events in "forward" and "backward" histograms. This allows "normal" uSR
@@ -184,27 +190,31 @@ Double_t PSimulateMuTransition::PrecessionPhase(const Double_t &time, const Doub
 }
 
 //--------------------------------------------------------------------------
-// MuonEvent (private)
+// Event (private)
 //--------------------------------------------------------------------------
 /**
  * <p> Generates "one muon event": simulate muon phase under free precession at
  * external field and Mu precession
- * initial muon state: Mu+
  *
+ * \param muonString if eq. "Mu+" begin with Mu+ precession
  */
-void PSimulateMuTransition::MuonEvent()
+void PSimulateMuTransition::Event(const TString muonString)
 {
   Double_t eventTime, eventDiffTime, captureTime, ionizationTime;
-  Double_t muonPrecessionFreq; // MHz
+  Double_t muonPrecessionFreq, muoniumPrecessionFreq; // MHz
+  Double_t rndm, frac1, frac2;
 
   muonPrecessionFreq = fMuonGyroRatio * fBfield; 
 
   // charge-exchange loop until muon decay
   eventTime     = 0.;
   eventDiffTime = 0.;
-  if (fDebugFlag) cout << "Muon Event, Decay time = " << fMuonDecayTime << endl;
+
+  if (fDebugFlag) cout << "Decay time = " << fMuonDecayTime << endl;
+  //cout << muonString << endl;
   while (1) {
-   // assume Mu+ as 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;
@@ -219,55 +229,54 @@ void PSimulateMuTransition::MuonEvent()
      // now, we have Mu0; get next ionization time
      ionizationTime = NextEventTime(fIonizationRate);
      eventTime += ionizationTime;
-   if (fDebugFlag) cout << "Ioniza. time = " << ionizationTime << " Phase = " << fMuonPhase << endl;
+     // 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)
+       muoniumPrecessionFreq = fMuPrecFreq1;
+     else
+       muoniumPrecessionFreq = fMuPrecFreq2;
+
+     if (fDebugFlag) cout << "Ioniza. time = " << ionizationTime << " Freq = " << muoniumPrecessionFreq 
+                          << " Phase = " << fMuonPhase << endl;
      if (eventTime < fMuonDecayTime)
-     fMuonPhase += PrecessionPhase(ionizationTime, fMuCoupling);
+       fMuonPhase += PrecessionPhase(ionizationTime, muoniumPrecessionFreq);
      else{ //muon decays; handle precession prior to muon decay
        eventDiffTime = fMuonDecayTime - (eventTime - ionizationTime);
-     fMuonPhase += PrecessionPhase(eventDiffTime, fMuCoupling);
+       fMuonPhase += PrecessionPhase(eventDiffTime, muoniumPrecessionFreq);
        break;
      }
    }
-
+   else{
-  if (fDebugFlag) cout << " Final Phase = " << fMuonPhase << endl;
+     // Mu0 as initial state; get next ionization time
-  //fMuonPhase = TMath::ACos(TMath::Cos(fMuonPhase))*360./TMath::TwoPi(); //transform back to [0, 180] degree interval
-  return;
-}
-
-//--------------------------------------------------------------------------
-// MuoniumEvent (private)
-//--------------------------------------------------------------------------
-/**
- * <p> Generates "one muonium event": simulate muon spin phase in Mu and as
- * free muon in external field
- * initial muon state: Mu0
- *
- */
-void PSimulateMuTransition::MuoniumEvent()
-{
-  Double_t eventTime, eventDiffTime, captureTime, ionizationTime;
-  Double_t muonPrecessionFreq; // MHz
-
-  muonPrecessionFreq = fMuonGyroRatio * fBfield; 
-
-  // charge-exchange loop until muon decay
-  eventTime     = 0.;
-  eventDiffTime = 0.;
-  if (fDebugFlag) cout << "Muonium event, Decay time = " << fMuonDecayTime << endl;
-  while (1) {
-   // we have Mu0 as initial state; get next ionization time
      ionizationTime = NextEventTime(fIonizationRate);
      eventTime += ionizationTime;
-   if (fDebugFlag) cout << "Ioniza. time = " << ionizationTime << " Phase = " << fMuonPhase << endl;
+     // 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)
+       muoniumPrecessionFreq = fMuPrecFreq1;
+     else
+       muoniumPrecessionFreq = fMuPrecFreq2;
+
+     if (fDebugFlag) 
+      cout << "Mu Ioniza. time = " << ionizationTime << " Freq = " << muoniumPrecessionFreq 
+         << " Phase = " << fMuonPhase << endl;
      if (eventTime < fMuonDecayTime)
-     fMuonPhase += PrecessionPhase(ionizationTime, fMuCoupling);
+       fMuonPhase += PrecessionPhase(ionizationTime, muoniumPrecessionFreq);
      else{ //muon decays; handle precession prior to muon decay
        eventDiffTime = fMuonDecayTime - (eventTime - ionizationTime);
-     fMuonPhase += PrecessionPhase(eventDiffTime, fMuCoupling);
+       fMuonPhase += PrecessionPhase(eventDiffTime, muoniumPrecessionFreq);
        break;
      }
 
-   // now we have Mu+, get next electron capture time
+     // Mu+ state; get next electron capture time
      captureTime = NextEventTime(fCaptureRate);
      eventTime += captureTime;
      if (fDebugFlag) cout << "Capture time = " << captureTime << " Phase = " << fMuonPhase << endl;
@@ -279,6 +288,7 @@ void PSimulateMuTransition::MuoniumEvent()
        break;
      }
    }
+  }
 
   if (fDebugFlag) cout << " Final Phase = " << fMuonPhase << endl;
   //fMuonPhase = TMath::ACos(TMath::Cos(fMuonPhase))*360./TMath::TwoPi(); //transform back to [0, 180] degree interval

src/tests/MuTransition/PSimulateMuTransition.h

@@ -49,17 +49,19 @@ 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 SetMuCoupling(Double_t value) { fMuCoupling = value; }        //!< sets Mu hyperfine coupling (MHz)
+    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 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
     virtual void SetMuFraction(Double_t value){ fMuFraction = value; }         //!< Muonium fraction
+    virtual void SetMuFractionState1(Double_t value){ fMuFractionState1 = value; }  
+    virtual void SetMuFractionState2(Double_t value){ fMuFractionState2 = value; }  
 
     virtual Bool_t IsValid() { return fValid; }
     virtual void   SetSeed(UInt_t seed);
 
     virtual Double_t GetBfield()      { return fBfield; }            //!< returns the magnetic field (T)
-    virtual Double_t GetMuCoupling()  { return fMuCoupling; }        //!< returns the Mu hyperfine coupling (MHz)
     virtual Double_t GetCaptureRate() { return fCaptureRate; }       //!< returns Mu+ electron capture rate (MHz)
     virtual Double_t GetIonizationRate() { return fIonizationRate; } //!< returns Mu0 ionization rate (MHz)
     virtual void     Run(TH1F *histoForward, TH1F *histoBackward);
@@ -69,21 +71,23 @@ class PSimulateMuTransition : public TObject
     TRandom2 *fRandom;
 
     Double_t  fBfield;          //!< magnetic field (T)
-    Double_t  fMuCoupling;      //!< Mu hyperfine coupling constant (MHz)
+    Double_t  fMuPrecFreq1;     //!< Mu precession frequency of state 1 (MHz)
+    Double_t  fMuPrecFreq2;     //!< Mu precession frequency of state 2 (MHz)
     Double_t  fCaptureRate;     //!< Mu+ electron capture rate (MHz)
     Double_t  fIonizationRate;  //!< Mu0 ionization rate (MHz)
     Double_t  fInitialPhase;    //!< initial muon spin phase
     Double_t  fMuonDecayTime;   //!< muon decay time (us)
     Double_t  fMuonPhase;       //!< phase of muon spin
     Double_t  fAsymmetry;       //!< muon decay asymmetry
-    Double_t  fMuFraction;      //!< Mu fraction [0,1]
+    Double_t  fMuFraction;      //!< total Mu fraction [0,1]
+    Double_t  fMuFractionState1; //!< fraction of Mu in state 1
+    Double_t  fMuFractionState2; //!< fraction of Mu in state 2
     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 void     MuonEvent();
+    virtual void     Event(const TString muonString);
-    virtual void     MuoniumEvent();
  
   ClassDef(PSimulateMuTransition, 0)
 };

src/tests/MuTransition/runMuSimulation.C

@@ -43,7 +43,7 @@ void runMuSimulation()
   decayAnaModule = histosFolder->AddFolder("DecayAnaModule", "muSR decay histograms");
 
   // feed run info header
-  UInt_t runNo = 9016;
+  UInt_t runNo = 9100;
   TString tstr;
   runInfo = gROOT->GetRootFolder()->AddFolder("RunInfo", "LEM RunInfo");  
   gROOT->GetListOfBrowsables()->Add(runInfo, "RunInfo");
@@ -61,7 +61,7 @@ void runMuSimulation()
   header->SetImpEnergy(31.8);
   header->SetSampleTemperature(0.2, 0.001);
   header->SetSampleBField(-1.0, 0.1);
-  header->SetTimeResolution(0.1953125);
+  header->SetTimeResolution(1.);
   header->SetNChannels(12001);
   header->SetNHist(2);
   header->SetCuts("none");
@@ -86,13 +86,16 @@ void runMuSimulation()
   }
 
   //prepare to run simulation
-  simulateMuTransition->SetMuCoupling(35.);      // MHz
+  simulateMuTransition->SetMuPrecFreq1(51.);   // MHz
+  simulateMuTransition->SetMuPrecFreq2(-27.);  // MHz
+  simulateMuTransition->SetMuFraction(0.5);    // initial Mu fraction
+  simulateMuTransition->SetMuFractionState1(1.0);  // 100% of Mu in state 1
+  simulateMuTransition->SetMuFractionState2(0.0);  // 0% of Mu in state 2
   simulateMuTransition->SetBfield(0.1);           // Tesla
-  simulateMuTransition->SetCaptureRate(1.0);       // MHz
+  simulateMuTransition->SetCaptureRate(1.5);       // MHz
-  simulateMuTransition->SetIonizationRate(250.0);    // MHz
+  simulateMuTransition->SetIonizationRate(250.);    // MHz
-  simulateMuTransition->SetNmuons(1e6);
+  simulateMuTransition->SetNmuons(1e7);
   simulateMuTransition->SetDecayAsymmetry(0.27);
-  simulateMuTransition->SetMuFraction(0.5);
   simulateMuTransition->SetDebugFlag(kFALSE); // to print time and phase during charge-changing cycle
 
   simulateMuTransition->PrintSettings();