diff --git a/src/tests/MuTransition/PSimulateMuTransition.cpp b/src/tests/MuTransition/PSimulateMuTransition.cpp index d84dfe5d..3ec4def7 100644 --- a/src/tests/MuTransition/PSimulateMuTransition.cpp +++ b/src/tests/MuTransition/PSimulateMuTransition.cpp @@ -100,10 +100,13 @@ PSimulateMuTransition::PSimulateMuTransition(UInt_t seed) } fNmuons = 100; // number of muons to simulate + fNshowProgress = 100; // print progress on screen every fNshowProgress events 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 + fMuPrecFreq13 = 0.; // Mu precession frequency of a 13 transition + fMuPrecFreq24 = 0.; // Mu precession frequency of a 24 transition fMuonPrecFreq = 0.; // muon precession frequency fBfield = 0.01; // magnetic field (T) fCaptureRate = 0.01; // Mu+ capture rate (MHz) @@ -118,6 +121,9 @@ PSimulateMuTransition::PSimulateMuTransition(UInt_t seed) fMuFractionState34 = 0.25; fMuFractionState23 = 0.25; fMuFractionState14 = 0.25; + fMuFractionState13 = 0.; + fMuFractionState24 = 0.; + fDebugFlag = kFALSE; } @@ -147,6 +153,8 @@ void PSimulateMuTransition::PrintSettings() const cout << endl << "Mu0 precession frequency 34 (MHz) = " << fMuPrecFreq34; cout << endl << "Mu0 precession frequency 23 (MHz) = " << fMuPrecFreq23; cout << endl << "Mu0 precession frequency 14 (MHz) = " << fMuPrecFreq14; + cout << endl << "Mu0 precession frequency 13 (MHz) = " << fMuPrecFreq13; + cout << endl << "Mu0 precession frequency 24 (MHz) = " << fMuPrecFreq24; cout << endl << "Mu+ precession frequency (MHz) = " << fMuonGyroRatio * fBfield; cout << endl << "B field (T) = " << fBfield; cout << endl << "Mu+ electron capture rate (MHz) = " << fCaptureRate; @@ -165,6 +173,7 @@ void PSimulateMuTransition::PrintSettings() const cout << endl << "Muonium fraction state23 = " << fMuFractionState23; cout << endl << "Muonium fraction state14 = " << fMuFractionState14; cout << endl << "Number of particles to simulate = " << fNmuons; + cout << endl << "Print progress on screen frequency = " << fNshowProgress; cout << endl << "Initial muon spin phase (degree) = " << fInitialPhase; cout << endl << "Debug flag = " << fDebugFlag; cout << endl << endl; @@ -224,7 +233,7 @@ void PSimulateMuTransition::Run(TH1F *histoForward, TH1F *histoBackward) else histoBackward->Fill(fMuonDecayTime*1000., 1. - fAsymmetry*TMath::Cos(fMuonPhase)); - if ( (i%100000) == 0) cout << "number of events processed: " << i << endl; + if ( (i%fNshowProgress) == 0) cout << "number of events processed: " << i << endl; } cout << "number of events processed: " << i << endl; return; @@ -291,11 +300,14 @@ TComplex PSimulateMuTransition::GTFunction(const Double_t &time, const TString c complexPol = TComplex::Exp(-TComplex::I()*twoPi*fMuonPrecFreq*time); else{ complexPol = - (fMuFractionState12 * TComplex::Exp(TComplex::I()*twoPi*fMuPrecFreq12*time) + + (fMuFractionState12 * TComplex::Exp(+TComplex::I()*twoPi*fMuPrecFreq12*time) + fMuFractionState34 * TComplex::Exp(-TComplex::I()*twoPi*fMuPrecFreq34*time)) + - (fMuFractionState23 * TComplex::Exp(TComplex::I()*twoPi*fMuPrecFreq23*time) + - fMuFractionState14 * TComplex::Exp(TComplex::I()*twoPi*fMuPrecFreq14*time)); + (fMuFractionState23 * TComplex::Exp(+TComplex::I()*twoPi*fMuPrecFreq23*time) + + fMuFractionState14 * TComplex::Exp(+TComplex::I()*twoPi*fMuPrecFreq14*time)) + + + (fMuFractionState13 * TComplex::Exp(+TComplex::I()*twoPi*fMuPrecFreq13*time) + + fMuFractionState24 * TComplex::Exp(+TComplex::I()*twoPi*fMuPrecFreq24*time)); } return complexPol; @@ -354,9 +366,6 @@ Double_t PSimulateMuTransition::GTSpinFlip(const Double_t &time) * after ionization process is given by Px(t_i+1)*Px(t_i). * 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 fMuFractionState12 - * ni_23 and nu_14 with equal probabilities, probability for both states fMuFractionState23 * *

Calculates Mu0 polarization in x direction during cyclic charge exchange. * See M. Senba, J.Phys. B23, 1545 (1990), equations (9), (11) diff --git a/src/tests/MuTransition/PSimulateMuTransition.h b/src/tests/MuTransition/PSimulateMuTransition.h index 552ff132..aab0e8d3 100644 --- a/src/tests/MuTransition/PSimulateMuTransition.h +++ b/src/tests/MuTransition/PSimulateMuTransition.h @@ -47,13 +47,16 @@ 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 SetNmuons(Int_t value) { fNmuons = value; } //!< number of muons + virtual void SetNshowProgress(Int_t value) { fNshowProgress = value; } //!< frequency of output on screen how many muons have been processed + virtual void SetDebugFlag(Bool_t value) { fDebugFlag = value; } //!< debug flag + virtual void SetBfield(Double_t value) { fBfield = value; } //!< sets magnetic field (T) 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 SetMuPrecFreq13(Double_t value) { fMuPrecFreq13 = value; } //!< sets Mu transition frequency (MHz) + virtual void SetMuPrecFreq24(Double_t value) { fMuPrecFreq24 = 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 SetSpinFlipRate(Double_t value){ fSpinFlipRate = value; } //!< sets Mu0 spin flip rate (MHz) @@ -63,6 +66,8 @@ class PSimulateMuTransition : public TObject virtual void SetMuFractionState34(Double_t value){ fMuFractionState34 = value; } virtual void SetMuFractionState23(Double_t value){ fMuFractionState23 = value; } virtual void SetMuFractionState14(Double_t value){ fMuFractionState14 = value; } + virtual void SetMuFractionState13(Double_t value){ fMuFractionState13 = value; } + virtual void SetMuFractionState24(Double_t value){ fMuFractionState24 = value; } virtual Bool_t IsValid() { return fValid; } virtual void SetSeed(UInt_t seed); @@ -81,6 +86,8 @@ class PSimulateMuTransition : public TObject 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 fMuPrecFreq13; //!< Mu transition frequency 13 (MHz) + Double_t fMuPrecFreq24; //!< Mu transition frequency 24 (MHz) Double_t fMuonPrecFreq; //!< muon precession frequency (MHz) Double_t fCaptureRate; //!< Mu+ electron capture rate (MHz) Double_t fIonizationRate; //!< Mu0 ionization rate (MHz) @@ -94,8 +101,11 @@ class PSimulateMuTransition : public TObject Double_t fMuFractionState34; //!< fraction of Mu in state 34 Double_t fMuFractionState23; //!< fraction of Mu in state 23 Double_t fMuFractionState14; //!< fraction of Mu in state 14 - Int_t fNmuons; //!< number of muons to simulate - Bool_t fDebugFlag; //!< debug flag + Double_t fMuFractionState13; //!< fraction of Mu in state 13 + Double_t fMuFractionState24; //!< fraction of Mu in state 24 + Int_t fNmuons; //!< number of muons to simulate + Int_t fNshowProgress; //!< output on screen how many muons have been processed + Bool_t fDebugFlag; //!< debug flag virtual Double_t NextEventTime(const Double_t &EventRate); // virtual Double_t PrecessionPhase(const Double_t &time, const TString chargeState); diff --git a/src/tests/MuTransition/runMuSimulation.C b/src/tests/MuTransition/runMuSimulation.C index bb2b3606..748ba674 100644 --- a/src/tests/MuTransition/runMuSimulation.C +++ b/src/tests/MuTransition/runMuSimulation.C @@ -63,29 +63,37 @@ void runMuSimulation() Double_t T = 300.; //temperature Double_t EA = 100; //activation energy (meV) Double_t spinFlipRate = 0.01; //if spinFlipRate > 0.001 only spin-flip processes will be simulated - Double_t capRate = 0.0001;//*sqrt(T/200.); //assume that capture rate varies as sqrt(T), capRate = sigma*v*p , v ~ sqrt(T) + Double_t capRate = 0.001;//*sqrt(T/200.); //assume that capture rate varies as sqrt(T), capRate = sigma*v*p , v ~ sqrt(T) + Double_t preFac = 6.7e7; Double_t ionRate; //assume Arrhenius behaviour ionRate = preFac*exp(-EA/kT) - ionRate = 0.1; //2.9e7 * exp(-EA/(0.08625*T)); // Ge: 2.9*10^7MHz "attempt" frequency; 1K = 0.08625 meV - Double_t B = 106.5; //field in G - Double_t Bvar = 0.; //field variance - Double_t Freq12 = 40.433; //Mu freq of the 12 transition - Double_t Freq34 = 59.567; //Mu freq of the 34 transition - Double_t Freq23 = 256.245; //Mu freq of the 23 transition - Double_t Freq14 = 356.245; //Mu freq of the 14 transition - Double_t MuFrac = 1.0; //total Mu fraction - Double_t MuFrac12 = 0.487; //weight of transition 12 - Double_t MuFrac34 = 0.487; //weight of transition 34 - Double_t MuFrac23 = 0.013; //weight of transition 23 - Double_t MuFrac14 = 0.013; //weight of transition 14 - Int_t Nmuons = 5e6; //number of muons - Double_t Asym = 0.27; //muon decay asymmetry - Int_t debugFlag = 0; //print debug information on screen + ionRate = 0.001; //preFac * exp(-EA/(0.08625*T)); // Ge: 2.9*10^7MHz "attempt" frequency; 1K = 0.08625 meV + Double_t B = 100.0; //field in G + Double_t Bvar = 0.; //field variance + Double_t Freq12 = 40.023; //Mu freq of the 12 transition + Double_t Freq34 = 59.977; //Mu freq of the 34 transition + Double_t Freq23 = 238.549; //Mu freq of the 23 transition + Double_t Freq14 = 338.549; //Mu freq of the 14 transition + Double_t Freq13 = 278.571; //Mu freq of the 23 transition + Double_t Freq24 = 325.165; //Mu freq of the 14 transition + + Double_t MuFrac = 1.0; //total Mu fraction + Double_t MuFrac12 = 0.486; //weight of transition 12 + Double_t MuFrac34 = 0.486; //weight of transition 34 + Double_t MuFrac23 = 0.014; //weight of transition 23 + Double_t MuFrac14 = 0.014; //weight of transition 14 + Double_t MuFrac13 = 0.0; //weight of transition 13 + Double_t MuFrac24 = 0.0; //weight of transition 24 + + Int_t Nmuons = 5e6; //number of muons + Int_t NshowProgress = 1e4; //frequency to show progress on screen + Double_t Asym = 0.27; //muon decay asymmetry + Int_t debugFlag = 0; //print debug information on screen histogramFileName = TString("0"); histogramFileName += runNo; histogramFileName += TString(".root"); - sprintf(titleStr,"- complexMuPol, A0 100MHz, Mu-frac %3.2f, Mu12 %6.2f MHz(%3.2f), Mu23 %6.2f MHz(%3.2f), ionRate %8.3f MHz, capRate %6.3f MHz, SF rate %6.3f MHz, %5.1f G", MuFrac, Freq12, MuFrac12/2, Freq23, MuFrac23/2, ionRate, capRate, spinFlipRate, B); + sprintf(titleStr,"- complexMuPol, A0 100MHz, Mu-frac %3.2f, Mu12 %6.2f MHz(%3.2f), Mu23 %6.2f MHz(%3.2f), ionRate %8.3f MHz, capRate %6.3f MHz, SF rate %6.3f MHz, %5.1f G", MuFrac, Freq12, MuFrac12, Freq23, MuFrac23, ionRate, capRate, spinFlipRate, B); runTitle = TString("0"); runTitle += runNo; runTitle += TString(titleStr); @@ -99,15 +107,20 @@ void runMuSimulation() simulateMuTransition->SetMuPrecFreq34(Freq34); // MHz simulateMuTransition->SetMuPrecFreq23(Freq23); // MHz simulateMuTransition->SetMuPrecFreq14(Freq14); // MHz + simulateMuTransition->SetMuPrecFreq13(Freq13); // MHz + simulateMuTransition->SetMuPrecFreq24(Freq24); // MHz simulateMuTransition->SetMuFraction(MuFrac); // initial Mu fraction simulateMuTransition->SetMuFractionState12(MuFrac12); simulateMuTransition->SetMuFractionState34(MuFrac34); simulateMuTransition->SetMuFractionState23(MuFrac23); - simulateMuTransition->SetMuFractionState14(MuFrac14); + simulateMuTransition->SetMuFractionState14(MuFrac14); + simulateMuTransition->SetMuFractionState13(MuFrac13); + simulateMuTransition->SetMuFractionState24(MuFrac24); simulateMuTransition->SetBfield(B/10000.); // Tesla simulateMuTransition->SetCaptureRate(capRate); // MHz simulateMuTransition->SetIonizationRate(ionRate); // MHz simulateMuTransition->SetSpinFlipRate(spinFlipRate); // MHz + simulateMuTransition->SetNshowProgress(NshowProgress); simulateMuTransition->SetNmuons(Nmuons); simulateMuTransition->SetDecayAsymmetry(Asym); simulateMuTransition->SetDebugFlag(debugFlag); // to print time and phase during charge-changing cycle @@ -179,12 +192,18 @@ void runMuSimulation() header->Set("Simulation/Mu0 Precession frequency 34", Freq34); header->Set("Simulation/Mu0 Precession frequency 23", Freq23); header->Set("Simulation/Mu0 Precession frequency 14", Freq14); + header->Set("Simulation/Mu0 Precession frequency 13", Freq13); + header->Set("Simulation/Mu0 Precession frequency 24", Freq24); header->Set("Simulation/Mu0 Fraction", MuFrac); header->Set("Simulation/Mu0 Fraction 12", MuFrac12); header->Set("Simulation/Mu0 Fraction 34", MuFrac34); header->Set("Simulation/Mu0 Fraction 23", MuFrac23); header->Set("Simulation/Mu0 Fraction 14", MuFrac14); - header->Set("Simulation/muon Capture Rate", capRate); + header->Set("Simulation/Mu0 Fraction 13", MuFrac13); + header->Set("Simulation/Mu0 Fraction 24", MuFrac24); + header->Set("Simulation/Mu0 Activation Energy", EA); + header->Set("Simulation/Mu0 Activation PreFactor", preFac); + header->Set("Simulation/Mux Capture Rate", capRate); header->Set("Simulation/Mu0 Ionization Rate", ionRate); header->Set("Simulation/Mu0 Spin Flip Rate", spinFlipRate); header->Set("Simulation/Number of Muons", Nmuons); @@ -206,7 +225,7 @@ void runMuSimulation() histo[i] = new TH1F(str, str, 18001, -0.5, 18000.5); } - for (i=0; iAdd(histo[i]); // run simulation