From 07e6863939d1906ab81be277b8f92a37c0bf9c81 Mon Sep 17 00:00:00 2001 From: nemu Date: Mon, 1 Mar 2010 09:26:19 +0000 Subject: [PATCH] Added Mu0 as initial state with user defined fraction. --- .../MuTransition/PSimulateMuTransition.cpp | 76 ++++++++++++++++--- .../MuTransition/PSimulateMuTransition.h | 5 +- src/tests/MuTransition/runMuSimulation.C | 3 +- 3 files changed, 72 insertions(+), 12 deletions(-) diff --git a/src/tests/MuTransition/PSimulateMuTransition.cpp b/src/tests/MuTransition/PSimulateMuTransition.cpp index 34f0f4fe..f3a4e831 100644 --- a/src/tests/MuTransition/PSimulateMuTransition.cpp +++ b/src/tests/MuTransition/PSimulateMuTransition.cpp @@ -63,6 +63,7 @@ PSimulateMuTransition::PSimulateMuTransition(UInt_t seed) fMuonPhase = fInitialPhase; fMuonDecayTime = 0.; fAsymmetry = 0.27; + fMuFraction = 0.; fDebugFlag = kFALSE; } @@ -92,9 +93,11 @@ void PSimulateMuTransition::PrintSettings() const 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 << "Number of particles to simulate = " << fNmuons; cout << endl << "Initial muon spin phase (degree) = " << fInitialPhase; - cout << endl << "Decay asymmetry = " << fAsymmetry; + cout << endl << "Debug flag = " << fDebugFlag; cout << endl << endl; } @@ -127,8 +130,14 @@ void PSimulateMuTransition::Run(TH1F *histoForward, TH1F *histoBackward) return; for (i = 0; iRndm() <= 1.-fMuFraction) + MuonEvent(); + else + MuoniumEvent(); + // fill 50% in "forward", and 50% in "backward" detector to get independent // events in "forward" and "backward" histograms. This allows "normal" uSR // analysis of the data @@ -175,29 +184,25 @@ Double_t PSimulateMuTransition::PrecessionPhase(const Double_t &time, const Doub } //-------------------------------------------------------------------------- -// Event (private) +// MuonEvent (private) //-------------------------------------------------------------------------- /** *

Generates "one muon event": simulate muon phase under free precession at * external field and Mu precession * initial muon state: Mu+ * - * \param muonDecayTime muon decay time (us) - * \param muonPhase muon spin phase in [0,180] degree */ -void PSimulateMuTransition::Event() +void PSimulateMuTransition::MuonEvent() { Double_t eventTime, eventDiffTime, captureTime, ionizationTime; Double_t muonPrecessionFreq; // MHz muonPrecessionFreq = fMuonGyroRatio * fBfield; - fMuonPhase = TMath::TwoPi() * fInitialPhase/360.; // transform to radians - fMuonDecayTime = NextEventTime(fMuonDecayRate); // charge-exchange loop until muon decay eventTime = 0.; eventDiffTime = 0.; - if (fDebugFlag) cout << "Decay time = " << fMuonDecayTime << endl; + if (fDebugFlag) cout << "Muon Event, Decay time = " << fMuonDecayTime << endl; while (1) { // assume Mu+ as initial state; get next electron capture time captureTime = NextEventTime(fCaptureRate); @@ -228,3 +233,54 @@ void PSimulateMuTransition::Event() //fMuonPhase = TMath::ACos(TMath::Cos(fMuonPhase))*360./TMath::TwoPi(); //transform back to [0, 180] degree interval return; } + +//-------------------------------------------------------------------------- +// MuoniumEvent (private) +//-------------------------------------------------------------------------- +/** + *

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; + if (eventTime < fMuonDecayTime) + fMuonPhase += PrecessionPhase(ionizationTime, fMuCoupling); + else{ //muon decays; handle precession prior to muon decay + eventDiffTime = fMuonDecayTime - (eventTime - ionizationTime); + fMuonPhase += PrecessionPhase(eventDiffTime, fMuCoupling); + break; + } + + // now we have Mu+, get next electron capture time + captureTime = NextEventTime(fCaptureRate); + eventTime += captureTime; + if (fDebugFlag) cout << "Capture time = " << captureTime << " Phase = " << fMuonPhase << endl; + if (eventTime < fMuonDecayTime) + fMuonPhase += PrecessionPhase(captureTime, muonPrecessionFreq); + else{ //muon decays; handle precession prior to muon decay + eventDiffTime = fMuonDecayTime - (eventTime - captureTime); + fMuonPhase += PrecessionPhase(eventDiffTime, muonPrecessionFreq); + break; + } + } + + if (fDebugFlag) cout << " Final Phase = " << fMuonPhase << endl; + //fMuonPhase = TMath::ACos(TMath::Cos(fMuonPhase))*360./TMath::TwoPi(); //transform back to [0, 180] degree interval + return; +} diff --git a/src/tests/MuTransition/PSimulateMuTransition.h b/src/tests/MuTransition/PSimulateMuTransition.h index 4c3034cd..a0b7bf59 100644 --- a/src/tests/MuTransition/PSimulateMuTransition.h +++ b/src/tests/MuTransition/PSimulateMuTransition.h @@ -53,6 +53,7 @@ class PSimulateMuTransition : public TObject 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 Bool_t IsValid() { return fValid; } virtual void SetSeed(UInt_t seed); @@ -75,12 +76,14 @@ class PSimulateMuTransition : public TObject 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] 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 Event(); + virtual void MuonEvent(); + virtual void MuoniumEvent(); ClassDef(PSimulateMuTransition, 0) }; diff --git a/src/tests/MuTransition/runMuSimulation.C b/src/tests/MuTransition/runMuSimulation.C index 9536e3fc..625b47b1 100644 --- a/src/tests/MuTransition/runMuSimulation.C +++ b/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 = 9015; + UInt_t runNo = 9016; TString tstr; runInfo = gROOT->GetRootFolder()->AddFolder("RunInfo", "LEM RunInfo"); gROOT->GetListOfBrowsables()->Add(runInfo, "RunInfo"); @@ -92,6 +92,7 @@ void runMuSimulation() simulateMuTransition->SetIonizationRate(250.0); // MHz simulateMuTransition->SetNmuons(1e6); simulateMuTransition->SetDecayAsymmetry(0.27); + simulateMuTransition->SetMuFraction(0.5); simulateMuTransition->SetDebugFlag(kFALSE); // to print time and phase during charge-changing cycle simulateMuTransition->PrintSettings();