diff --git a/src/tests/MuTransition/PSimulateMuTransition.cpp b/src/tests/MuTransition/PSimulateMuTransition.cpp index f3a4e831..f6826902 100644 --- a/src/tests/MuTransition/PSimulateMuTransition.cpp +++ b/src/tests/MuTransition/PSimulateMuTransition.cpp @@ -54,17 +54,20 @@ PSimulateMuTransition::PSimulateMuTransition(UInt_t seed) fValid = false; } - fNmuons = 100; // number of muons to simulate - fMuCoupling = 4463.; // vacuum Mu hyperfine coupling constant - fBfield = 0.01; // magnetic field (T) - fCaptureRate = 0.01; // Mu+ capture rate (MHz) - fIonizationRate = 10.; // Mu0 ionization rate (MHz) - fInitialPhase = 0.; - fMuonPhase = fInitialPhase; - fMuonDecayTime = 0.; - fAsymmetry = 0.27; - fMuFraction = 0.; - fDebugFlag = kFALSE; + fNmuons = 100; // number of muons to simulate + 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) + fInitialPhase = 0.; + fMuonPhase = fInitialPhase; + 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,99 +190,103 @@ Double_t PSimulateMuTransition::PrecessionPhase(const Double_t &time, const Doub } //-------------------------------------------------------------------------- -// MuonEvent (private) +// Event (private) //-------------------------------------------------------------------------- /** *

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; - while (1) { - // assume Mu+ as initial state; 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; - } - // now, we have Mu0; 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; - } - } - - if (fDebugFlag) cout << " Final Phase = " << fMuonPhase << endl; - //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 << "Decay time = " << fMuonDecayTime << endl; + //cout << muonString << endl; + while (1) { + 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 (eventTime < fMuonDecayTime) + fMuonPhase += PrecessionPhase(captureTime, muonPrecessionFreq); + else{ //muon decays; handle precession prior to muon decay + eventDiffTime = fMuonDecayTime - (eventTime - captureTime); + fMuonPhase += PrecessionPhase(eventDiffTime, muonPrecessionFreq); + 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) + muoniumPrecessionFreq = fMuPrecFreq1; + else + muoniumPrecessionFreq = fMuPrecFreq2; + + if (fDebugFlag) cout << "Ioniza. time = " << ionizationTime << " Freq = " << muoniumPrecessionFreq + << " Phase = " << fMuonPhase << endl; + if (eventTime < fMuonDecayTime) + fMuonPhase += PrecessionPhase(ionizationTime, muoniumPrecessionFreq); + else{ //muon decays; handle precession prior to muon decay + eventDiffTime = fMuonDecayTime - (eventTime - ionizationTime); + fMuonPhase += PrecessionPhase(eventDiffTime, muoniumPrecessionFreq); + break; + } + } + 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) + muoniumPrecessionFreq = fMuPrecFreq1; + else + muoniumPrecessionFreq = fMuPrecFreq2; + + if (fDebugFlag) + cout << "Mu Ioniza. time = " << ionizationTime << " Freq = " << muoniumPrecessionFreq + << " Phase = " << fMuonPhase << endl; + if (eventTime < fMuonDecayTime) + fMuonPhase += PrecessionPhase(ionizationTime, muoniumPrecessionFreq); + else{ //muon decays; handle precession prior to muon decay + eventDiffTime = fMuonDecayTime - (eventTime - ionizationTime); + fMuonPhase += PrecessionPhase(eventDiffTime, muoniumPrecessionFreq); + break; + } + + // Mu+ state; 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; + } } } diff --git a/src/tests/MuTransition/PSimulateMuTransition.h b/src/tests/MuTransition/PSimulateMuTransition.h index a0b7bf59..00716683 100644 --- a/src/tests/MuTransition/PSimulateMuTransition.h +++ b/src/tests/MuTransition/PSimulateMuTransition.h @@ -46,20 +46,22 @@ class PSimulateMuTransition : public TObject virtual ~PSimulateMuTransition(); virtual void PrintSettings() const; - 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 SetCaptureRate(Double_t value){ fCaptureRate = value; } //!< sets Mu+ electron capture rate (MHz) + 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 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,22 +71,24 @@ 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 MuoniumEvent(); - + virtual void Event(const TString muonString); + ClassDef(PSimulateMuTransition, 0) }; diff --git a/src/tests/MuTransition/runMuSimulation.C b/src/tests/MuTransition/runMuSimulation.C index 625b47b1..12d18388 100644 --- a/src/tests/MuTransition/runMuSimulation.C +++ b/src/tests/MuTransition/runMuSimulation.C @@ -32,7 +32,7 @@ void runMuSimulation() { // load library gSystem->Load("$ROOTSYS/lib/libPSimulateMuTransition"); - + // generate data TFolder *histosFolder; TFolder *decayAnaModule; @@ -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->SetIonizationRate(250.0); // MHz - simulateMuTransition->SetNmuons(1e6); + simulateMuTransition->SetCaptureRate(1.5); // MHz + simulateMuTransition->SetIonizationRate(250.); // MHz + 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();