From 1c31fc88d069dd8b192d8e9285c1e87d9e580f35 Mon Sep 17 00:00:00 2001 From: Thomas Prokscha Date: Wed, 17 Feb 2016 14:22:21 +0100 Subject: [PATCH 1/5] Added Mu0 spin-flip process --- .../MuTransition/PSimulateMuTransition.cpp | 824 ++++++++++-------- .../MuTransition/PSimulateMuTransition.h | 206 ++--- src/tests/MuTransition/runMuSimulation.C | 153 ++++ 3 files changed, 700 insertions(+), 483 deletions(-) diff --git a/src/tests/MuTransition/PSimulateMuTransition.cpp b/src/tests/MuTransition/PSimulateMuTransition.cpp index 15e558c1..8da00b00 100644 --- a/src/tests/MuTransition/PSimulateMuTransition.cpp +++ b/src/tests/MuTransition/PSimulateMuTransition.cpp @@ -1,382 +1,442 @@ -/*************************************************************************** - - PSimulateMuTransition.cpp - - Author: Thomas Prokscha - Date: 25-Feb-2010 - - $Id$ - - Use root macros runMuSimulation.C and testAnalysis.C to run the simulation - and to get a quick look on the data. Data are saved to a root histogram file - with a structure similar to LEM histogram files; musrfit can be used to - analyze the simulated data. - - Description: - Root class to simulate muon spin phase under successive Mu+/Mu0 charge-exchange - processes by a Monte-Carlo method. Consider transverse field geometry, and assume - initial muon spin direction in x, and field applied along z. For PxMu(t) in - muonium use the equation 8.22 of the muSR book of Yaounc and Dalmas de Réotier, in - slightly modified form (see Senba, J. Phys. B 23, 1545 (1990)); note that PxMu(t) - is given by a superposition of the four frequencies "nu_12", "nu_34", "nu_23", "nu_14". - These frequencies and the corresponding probabilities ("SetMuFractionState12" for - transitions 12 and 34, "SetMuFractionState23" for states 23 and 14) can be calculated - for a given field with the root macro AnisotropicMu.C - - Parameters: - 1) Precession frequencies of "nu_12", "nu_34", "nu_23", "nu_14" - 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 - 6) initial muon spin phase - 7) total muon decay asymmetry - 8) number of muon decays to be generated. - 9) debug flag: if TRUE print capture/ionization events on screen - - Output: - Two histograms ("forward" and "backward") are written to a root file. - - The muon event simulation with a sequence of charge-changing processes is - done in Event(): - simulate muon spin phase under charge-exchange with "4 Mu transitions" - 1) according to Mu+/Mu0 fraction begin either with a Mu+ state or Mu state - 2) Mu+: determine next electron-capture time t_c. If t_c is larger than decay time t_d - calculate muon spin precession for t_d; else calculate spin precession for t_c. - 3) Determine next ionization time t_i; calculate Px(t_i) in Muonium; calculate the - muon spin phase by acos(Px(t_i)). - 4) get the next electron capture time, continue until t_d is reached; accumulate muon spin - phase. - -***************************************************************************/ - -/*************************************************************************** - * Copyright (C) 2010 by Thomas Prokscha, Paul Scherrer Institut * - * * - * This program is free software; you can redistribute it and/or modify * - * it under the terms of the GNU General Public License as published by * - * the Free Software Foundation; either version 2 of the License, or * - * (at your option) any later version. * - * * - * This program is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * - * GNU General Public License for more details. * - * * - * You should have received a copy of the GNU General Public License * - * along with this program; if not, write to the * - * Free Software Foundation, Inc., * - * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * - ***************************************************************************/ - -#include -using namespace std; - -#include - -#include "PSimulateMuTransition.h" - -ClassImp(PSimulateMuTransition) - -//-------------------------------------------------------------------------- -// Constructor -//-------------------------------------------------------------------------- -/** - *

Constructor. - * - * \param seed for the random number generator - */ -PSimulateMuTransition::PSimulateMuTransition(UInt_t seed) -{ - fValid = true; - - fRandom = new TRandom2(seed); - if (fRandom == 0) { - fValid = false; - } - - fNmuons = 100; // number of muons to simulate - 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 - fMuonPrecFreq = 0.; // muon precession frequency - 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.; - fMuFractionState12 = 0.; - fMuFractionState23 = 0.; - fDebugFlag = kFALSE; -} - -//-------------------------------------------------------------------------- -// Destructor -//-------------------------------------------------------------------------- -/** - *

Destructor. - * - */ -PSimulateMuTransition::~PSimulateMuTransition() -{ - if (fRandom) { - delete fRandom; - fRandom = 0; - } -} -//-------------------------------------------------------------------------- -// Output of current settings -//-------------------------------------------------------------------------- -/*! - *

Prints the current settings onto std output. - */ -void PSimulateMuTransition::PrintSettings() const -{ - cout << endl << "Mu precession frequency 12 (MHz) = " << fMuPrecFreq12; - 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; - cout << endl << "Decay asymmetry = " << fAsymmetry; - cout << endl << "Muonium fraction = " << fMuFraction; - cout << endl << "Muonium fraction state12 = " << fMuFractionState12; - cout << endl << "Muonium fraction state23 = " << fMuFractionState23; - cout << endl << "Number of particles to simulate = " << fNmuons; - cout << endl << "Initial muon spin phase (degree) = " << fInitialPhase; - cout << endl << "Debug flag = " << fDebugFlag; - cout << endl << endl; -} - -//-------------------------------------------------------------------------- -// SetSeed (public) -//-------------------------------------------------------------------------- -/** - *

Sets the seed of the random generator. - * - * \param seed for the random number generator - */ -void PSimulateMuTransition::SetSeed(UInt_t seed) -{ - if (!fValid) - return; - - fRandom->SetSeed(seed); -} - -//-------------------------------------------------------------------------- -// Run (public) -//-------------------------------------------------------------------------- -/** - * \param histoForward - */ -void PSimulateMuTransition::Run(TH1F *histoForward, TH1F *histoBackward) -{ - Int_t i; - if (histoForward == 0 || histoBackward == 0) - return; - - for (i = 0; iRndm() <= 1.-fMuFraction) - Event("Mu+"); - else - Event(""); - - // 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 - // change muon decay time to ns - if (fRandom->Rndm() <= 0.5) - histoForward->Fill(fMuonDecayTime*1000., 1. + fAsymmetry*TMath::Cos(fMuonPhase)); - else - histoBackward->Fill(fMuonDecayTime*1000., 1. - fAsymmetry*TMath::Cos(fMuonPhase)); - - if ( (i%100000) == 0) cout << "number of events processed: " << i << endl; - } - cout << "number of events processed: " << i << endl; - return; -} - -//-------------------------------------------------------------------------- -// NextEventTime (private) -//-------------------------------------------------------------------------- -/** - *

Determine time of next event, assuming "Poisson" distribution in time - * - * \param EventRate event rate in MHz; returns next event time in micro-seconds - */ -Double_t PSimulateMuTransition::NextEventTime(const Double_t &EventRate) -{ - if (EventRate <= 0.) - return -1.; // signal error - - return -1./EventRate * TMath::Log(fRandom->Rndm()); -} - -//-------------------------------------------------------------------------- -// Phase (private) -//-------------------------------------------------------------------------- -/** - *

Determines phase of the muon spin - * - * \param time duration of precession (us); - * \param frequency muon spin precession frequency (MHz); - */ -Double_t PSimulateMuTransition::PrecessionPhase(const Double_t &time, const TString chargeState) -{ - Double_t muonPhaseX; - Double_t muoniumPolX = 0; - - if (chargeState == "Mu+") - muonPhaseX = TMath::TwoPi()*fMuonPrecFreq*time; - else if (chargeState == "Mu0"){ - muoniumPolX = 0.5 * - (fMuFractionState12 * (TMath::Cos(TMath::TwoPi()*fMuPrecFreq12*time) + TMath::Cos(TMath::TwoPi()*fMuPrecFreq34*time)) + - fMuFractionState23 * (TMath::Cos(TMath::TwoPi()*fMuPrecFreq23*time) + TMath::Cos(TMath::TwoPi()*fMuPrecFreq14*time))); - muonPhaseX = TMath::ACos(muoniumPolX); - } - else - muonPhaseX = 0.; - - return muonPhaseX; -} - -//-------------------------------------------------------------------------- -// Event (private) -//-------------------------------------------------------------------------- -/** - *

Generates "muon event": simulate muon spin phase under charge-exchange with - * a neutral muonium state in transverse field, where the polarization evolution - * PxMu(t) of the muon spin in muonium is determined by a superposition of the - * four "Mu transitions" nu_12, nu_34, nu_23, and nu_14. - * 1) according to Mu+/Mu0 fraction begin either with a Mu+ state or Mu state - * 2) Mu+: determine next electron-capture time t_c. If t_c is larger than decay time t_d - * calculate muon spin precession for t_d; else calculate spin precession for t_c. - * 3) Determine next ionization time t_i; calculate Px(t_i) in Muonium; calculate the - * muon spin phase by acos(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 - * - * \param muonString if eq. "Mu+" begin with Mu+ precession - */ -void PSimulateMuTransition::Event(const TString muonString) -{ - Double_t eventTime, eventDiffTime, captureTime, ionizationTime; -// Double_t muonPrecessionFreq, muoniumPrecessionFreq; // MHz -// Double_t rndm, frac1, frac2; - - fMuonPrecFreq = fMuonGyroRatio * fBfield; - - // charge-exchange loop until muon decay - eventTime = 0.; - eventDiffTime = 0.; - - 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, "Mu+"); - else{ //muon decays; handle precession prior to muon decay - eventDiffTime = fMuonDecayTime - (eventTime - captureTime); - fMuonPhase += PrecessionPhase(eventDiffTime, "Mu+"); - 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){ -// if (fRandom->Rndm() <= 0.5) -// muoniumPrecessionFreq = fMuPrecFreq12; -// else -// muoniumPrecessionFreq = fMuPrecFreq34; -// } -// else{ -// if (fRandom->Rndm() <= 0.5) -// muoniumPrecessionFreq = fMuPrecFreq23; -// else -// muoniumPrecessionFreq = fMuPrecFreq14; -// } - - if (fDebugFlag) cout << "Ioniza. time = " << ionizationTime << " Phase = " << fMuonPhase << endl; - if (eventTime < fMuonDecayTime) - fMuonPhase += PrecessionPhase(ionizationTime, "Mu0"); - else{ //muon decays; handle precession prior to muon decay - eventDiffTime = fMuonDecayTime - (eventTime - ionizationTime); - fMuonPhase += PrecessionPhase(eventDiffTime, "Mu0"); - 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){ -// if (fRandom->Rndm() <= 0.5) -// muoniumPrecessionFreq = fMuPrecFreq12; -// else -// muoniumPrecessionFreq = fMuPrecFreq34; -// } -// else{ -// if (fRandom->Rndm() <= 0.5) -// muoniumPrecessionFreq = fMuPrecFreq23; -// else -// muoniumPrecessionFreq = fMuPrecFreq14; -// } - - if (fDebugFlag) - cout << "Mu Ioniza. time = " << ionizationTime << " Phase = " << fMuonPhase << endl; - if (eventTime < fMuonDecayTime) - fMuonPhase += PrecessionPhase(ionizationTime, "Mu0"); - else{ //muon decays; handle precession prior to muon decay - eventDiffTime = fMuonDecayTime - (eventTime - ionizationTime); - fMuonPhase += PrecessionPhase(eventDiffTime, "Mu0"); - 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, "Mu+"); - else{ //muon decays; handle precession prior to muon decay - eventDiffTime = fMuonDecayTime - (eventTime - captureTime); - fMuonPhase += PrecessionPhase(eventDiffTime, "Mu+"); - 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; -} +/*************************************************************************** + + PSimulateMuTransition.cpp + + Author: Thomas Prokscha + Date: 25-Feb-2010 + + $Id$ + + Use root macros runMuSimulation.C and testAnalysis.C to run the simulation + and to get a quick look on the data. Data are saved to a root histogram file + with a structure similar to LEM histogram files; musrfit can be used to + analyze the simulated data. + + Description: + Root class to simulate muon spin phase under successive Mu+/Mu0 charge-exchange + processes by a Monte-Carlo method. Consider transverse field geometry, and assume + initial muon spin direction in x, and field applied along z. For PxMu(t) in + muonium use the equation 8.22 of the muSR book of Yaounc and Dalmas de Réotier, in + slightly modified form (see Senba, J. Phys. B 23, 1545 (1990)); note that PxMu(t) + is given by a superposition of the four frequencies "nu_12", "nu_34", "nu_23", "nu_14". + These frequencies and the corresponding probabilities ("SetMuFractionState12" for + transitions 12 and 34, "SetMuFractionState23" for states 23 and 14) can be calculated + for a given field with the root macro AnisotropicMu.C + + Parameters: + 1) Precession frequencies of "nu_12", "nu_34", "nu_23", "nu_14" + 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 + 6) initial muon spin phase + 7) total muon decay asymmetry + 8) number of muon decays to be generated. + 9) debug flag: if TRUE print capture/ionization events on screen + + Output: + Two histograms ("forward" and "backward") are written to a root file. + + The muon event simulation with a sequence of charge-changing processes is + done in Event(): + simulate muon spin phase under charge-exchange with "4 Mu transitions" + 1) according to Mu+/Mu0 fraction begin either with a Mu+ state or Mu state + 2) Mu+: determine next electron-capture time t_c. If t_c is larger than decay time t_d + calculate muon spin precession for t_d; else calculate spin precession for t_c. + 3) Determine next ionization time t_i; calculate Px(t_i) in Muonium; calculate the + muon spin phase by acos(Px(t_i)). + 4) get the next electron capture time, continue until t_d is reached; accumulate muon spin + phase. + +***************************************************************************/ + +/*************************************************************************** + * Copyright (C) 2010 by Thomas Prokscha, Paul Scherrer Institut * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program; if not, write to the * + * Free Software Foundation, Inc., * + * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * + ***************************************************************************/ + +#include +using namespace std; + +#include + +#include "PSimulateMuTransition.h" + +ClassImp(PSimulateMuTransition) + +//-------------------------------------------------------------------------- +// Constructor +//-------------------------------------------------------------------------- +/** + *

Constructor. + * + * \param seed for the random number generator + */ +PSimulateMuTransition::PSimulateMuTransition(UInt_t seed) +{ + fValid = true; + + fRandom = new TRandom2(seed); + if (fRandom == 0) { + fValid = false; + } + + fNmuons = 100; // number of muons to simulate + 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 + fMuonPrecFreq = 0.; // muon precession frequency + fBfield = 0.01; // magnetic field (T) + fCaptureRate = 0.01; // Mu+ capture rate (MHz) + fIonizationRate = 10.; // Mu0 ionization rate (MHz) + fSpinFlipRate = 0.001; // Mu0 spin flip rate (MHz) + fInitialPhase = 0.; + fMuonPhase = fInitialPhase; + fMuonDecayTime = 0.; + fAsymmetry = 0.27; + fMuFraction = 0.; + fMuFractionState12 = 0.; + fMuFractionState23 = 0.; + fDebugFlag = kFALSE; +} + +//-------------------------------------------------------------------------- +// Destructor +//-------------------------------------------------------------------------- +/** + *

Destructor. + * + */ +PSimulateMuTransition::~PSimulateMuTransition() +{ + if (fRandom) { + delete fRandom; + fRandom = 0; + } +} +//-------------------------------------------------------------------------- +// Output of current settings +//-------------------------------------------------------------------------- +/*! + *

Prints the current settings onto std output. + */ +void PSimulateMuTransition::PrintSettings() const +{ + cout << endl << "Mu precession frequency 12 (MHz) = " << fMuPrecFreq12; + 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 << "Mu0 ionizatioan rate (MHz) = " << fIonizationRate; + cout << endl << "Mu0 spin-flip rate (MHz) = " << fSpinFlipRate; + cout << endl << "!!! Note: if spin-flip rate > 0.001 only spin-flip process is considered!!!"; + cout << endl << "Decay asymmetry = " << fAsymmetry; + cout << endl << "Muonium fraction = " << fMuFraction; + cout << endl << "Muonium fraction state12 = " << fMuFractionState12; + cout << endl << "Muonium fraction state23 = " << fMuFractionState23; + cout << endl << "Number of particles to simulate = " << fNmuons; + cout << endl << "Initial muon spin phase (degree) = " << fInitialPhase; + cout << endl << "Debug flag = " << fDebugFlag; + cout << endl << endl; +} + +//-------------------------------------------------------------------------- +// SetSeed (public) +//-------------------------------------------------------------------------- +/** + *

Sets the seed of the random generator. + * + * \param seed for the random number generator + */ +void PSimulateMuTransition::SetSeed(UInt_t seed) +{ + if (!fValid) + return; + + fRandom->SetSeed(seed); +} + +//-------------------------------------------------------------------------- +// Run (public) +//-------------------------------------------------------------------------- +/** + * \param histoForward + */ +void PSimulateMuTransition::Run(TH1F *histoForward, TH1F *histoBackward) +{ + Int_t i; + if (histoForward == 0 || histoBackward == 0) + return; + + for (i = 0; i 0.001){// consider only Mu0 spin-flip in this case + fMuonPhase = TMath::ACos(GTSpinFlip(fMuonDecayTime)); + } + else{ + // initial muon state Mu+ or Mu0? + if (fRandom->Rndm() <= 1.-fMuFraction) + Event("Mu+"); + else + Event(""); + } + // 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 + // change muon decay time to ns + if (fRandom->Rndm() <= 0.5) + histoForward->Fill(fMuonDecayTime*1000., 1. + fAsymmetry*TMath::Cos(fMuonPhase)); + else + histoBackward->Fill(fMuonDecayTime*1000., 1. - fAsymmetry*TMath::Cos(fMuonPhase)); + + if ( (i%100000) == 0) cout << "number of events processed: " << i << endl; + } + cout << "number of events processed: " << i << endl; + return; +} + +//-------------------------------------------------------------------------- +// NextEventTime (private) +//-------------------------------------------------------------------------- +/** + *

Determine time of next event, assuming "Poisson" distribution in time + * + * \param EventRate event rate in MHz; returns next event time in micro-seconds + */ +Double_t PSimulateMuTransition::NextEventTime(const Double_t &EventRate) +{ + if (EventRate <= 0.) + return -1.; // signal error + + return -1./EventRate * TMath::Log(fRandom->Rndm()); +} + +//-------------------------------------------------------------------------- +// Phase (private) +//-------------------------------------------------------------------------- +/** + *

Determines phase of the muon spin + * + * \param time duration of precession (us); + * \param chargeState charge state of Mu ("Mu+" or "Mu0") + */ +Double_t PSimulateMuTransition::PrecessionPhase(const Double_t &time, const TString chargeState) +{ + Double_t muonPhaseX; + Double_t muoniumPolX = 0; + + if (chargeState == "Mu+") + muonPhaseX = TMath::TwoPi()*fMuonPrecFreq*time; + else if (chargeState == "Mu0"){ + muoniumPolX = GTFunction(time); + muonPhaseX = TMath::ACos(muoniumPolX); + } + else + muonPhaseX = 0.; + + return muonPhaseX; +} + +//-------------------------------------------------------------------------- +// Mu0 transverse field polarization function (private) +//-------------------------------------------------------------------------- +/** + *

Calculates Mu0 polarization in x direction by superposition of four Mu0 frequencies + * + * \param time (us); + */ +Double_t PSimulateMuTransition::GTFunction(const Double_t &time) +{ + Double_t muoniumPolX = 0; + + muoniumPolX = 0.5 * + (fMuFractionState12 * (TMath::Cos(TMath::TwoPi()*fMuPrecFreq12*time) + TMath::Cos(TMath::TwoPi()*fMuPrecFreq34*time)) + + fMuFractionState23 * (TMath::Cos(TMath::TwoPi()*fMuPrecFreq23*time) + TMath::Cos(TMath::TwoPi()*fMuPrecFreq14*time))); + + return muoniumPolX; +} + +//-------------------------------------------------------------------------- +// Mu0 transverse field polarization function after n spin-flip collisions (private) +//-------------------------------------------------------------------------- +/** + *

Calculates Mu0 polarization in x direction after n spin flip collisions. + * See M. Senba, J.Phys. B24, 3531 (1991), equation (17) + * + * \param time (us); + */ +Double_t PSimulateMuTransition::GTSpinFlip(const Double_t &time) +{ + Double_t muoniumPolX = 1.0; //initial polarization in x direction + Double_t eventTime = 0; + Double_t eventDiffTime = 0; + Double_t lastEventTime = 0; + + eventTime += NextEventTime(fSpinFlipRate); + if (eventTime >= time){ + muoniumPolX = GTFunction(time); + } + else{ + while (eventTime < time){ + eventDiffTime = eventTime - lastEventTime; + muoniumPolX = muoniumPolX * GTFunction(eventDiffTime); + lastEventTime = eventTime; + eventTime += NextEventTime(fSpinFlipRate); + } + // calculate for the last collision + eventDiffTime = time - lastEventTime; + muoniumPolX = muoniumPolX * GTFunction(eventDiffTime); + } + + return muoniumPolX; +} + +//-------------------------------------------------------------------------- +// Event (private) +//-------------------------------------------------------------------------- +/** + *

Generates "muon event": simulate muon spin phase under charge-exchange with + * a neutral muonium state in transverse field, where the polarization evolution + * PxMu(t) of the muon spin in muonium is determined by a superposition of the + * four "Mu transitions" nu_12, nu_34, nu_23, and nu_14. + * 1) according to Mu+/Mu0 fraction begin either with a Mu+ state or Mu state + * 2) Mu+: determine next electron-capture time t_c. If t_c is larger than decay time t_d + * calculate muon spin precession for t_d; else calculate spin precession for t_c. + * 3) Determine next ionization time t_i; calculate Px(t_i) in Muonium; calculate the + * muon spin phase by acos(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 + * + * \param muonString if eq. "Mu+" begin with Mu+ precession + */ +void PSimulateMuTransition::Event(const TString muonString) +{ + Double_t eventTime, eventDiffTime, captureTime, ionizationTime; +// Double_t muonPrecessionFreq, muoniumPrecessionFreq; // MHz +// Double_t rndm, frac1, frac2; + + fMuonPrecFreq = fMuonGyroRatio * fBfield; + + // charge-exchange loop until muon decay + eventTime = 0.; + eventDiffTime = 0.; + + 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, "Mu+"); + else{ //muon decays; handle precession prior to muon decay + eventDiffTime = fMuonDecayTime - (eventTime - captureTime); + fMuonPhase += PrecessionPhase(eventDiffTime, "Mu+"); + 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){ +// if (fRandom->Rndm() <= 0.5) +// muoniumPrecessionFreq = fMuPrecFreq12; +// else +// muoniumPrecessionFreq = fMuPrecFreq34; +// } +// else{ +// if (fRandom->Rndm() <= 0.5) +// muoniumPrecessionFreq = fMuPrecFreq23; +// else +// muoniumPrecessionFreq = fMuPrecFreq14; +// } + + if (fDebugFlag) cout << "Ioniza. time = " << ionizationTime << " Phase = " << fMuonPhase << endl; + if (eventTime < fMuonDecayTime) + fMuonPhase += PrecessionPhase(ionizationTime, "Mu0"); + else{ //muon decays; handle precession prior to muon decay + eventDiffTime = fMuonDecayTime - (eventTime - ionizationTime); + fMuonPhase += PrecessionPhase(eventDiffTime, "Mu0"); + 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){ +// if (fRandom->Rndm() <= 0.5) +// muoniumPrecessionFreq = fMuPrecFreq12; +// else +// muoniumPrecessionFreq = fMuPrecFreq34; +// } +// else{ +// if (fRandom->Rndm() <= 0.5) +// muoniumPrecessionFreq = fMuPrecFreq23; +// else +// muoniumPrecessionFreq = fMuPrecFreq14; +// } + + if (fDebugFlag) + cout << "Mu Ioniza. time = " << ionizationTime << " Phase = " << fMuonPhase << endl; + if (eventTime < fMuonDecayTime) + fMuonPhase += PrecessionPhase(ionizationTime, "Mu0"); + else{ //muon decays; handle precession prior to muon decay + eventDiffTime = fMuonDecayTime - (eventTime - ionizationTime); + fMuonPhase += PrecessionPhase(eventDiffTime, "Mu0"); + 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, "Mu+"); + else{ //muon decays; handle precession prior to muon decay + eventDiffTime = fMuonDecayTime - (eventTime - captureTime); + fMuonPhase += PrecessionPhase(eventDiffTime, "Mu+"); + 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 6640fd3b..4caac1c1 100644 --- a/src/tests/MuTransition/PSimulateMuTransition.h +++ b/src/tests/MuTransition/PSimulateMuTransition.h @@ -1,101 +1,105 @@ -/*************************************************************************** - - PSimulateMuTransition.h - - Author: Thomas Prokscha - Date: 25-Feb-2010 - - $Id$ - -***************************************************************************/ - -/*************************************************************************** - * Copyright (C) 2010 by Thomas Prokscha, Paul Scherrer Institut * - * * - * This program is free software; you can redistribute it and/or modify * - * it under the terms of the GNU General Public License as published by * - * the Free Software Foundation; either version 2 of the License, or * - * (at your option) any later version. * - * * - * This program is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * - * GNU General Public License for more details. * - * * - * You should have received a copy of the GNU General Public License * - * along with this program; if not, write to the * - * Free Software Foundation, Inc., * - * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * - ***************************************************************************/ - -#ifndef _PSIMULATEMUTRANSITION_H_ -#define _PSIMULATEMUTRANSITION_H_ - -#include -#include -#include - -// global constants -const Double_t fMuonGyroRatio = 135.54; //!< muon gyromagnetic ratio (MHz/T) -const Double_t fMuonDecayRate = 0.4551; //!< muon decay rate (1/tau_mu, MHz) - -class PSimulateMuTransition : public TObject -{ - public: - PSimulateMuTransition(UInt_t seed = 0); - 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 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 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 SetMuFractionState12(Double_t value){ fMuFractionState12 = value; } - virtual void SetMuFractionState23(Double_t value){ fMuFractionState23 = 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 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); - - private: - Bool_t fValid; - TRandom2 *fRandom; - - Double_t fBfield; //!< magnetic field (T) - Double_t fMuPrecFreq12; //!< Mu transition frequency 12 (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 fMuonPrecFreq; //!< muon precession frequency (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; //!< total Mu fraction [0,1] - Double_t fMuFractionState12; //!< fraction of Mu in state 12, 34 - Double_t fMuFractionState23; //!< fraction of Mu in state 23, 14 - 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 Double_t PrecessionPhase(const Double_t &time, const TString chargeState); - virtual void Event(const TString muonString); - - ClassDef(PSimulateMuTransition, 0) -}; - -#endif // _PSIMULATEMUTRANSITION_H_ +/*************************************************************************** + + PSimulateMuTransition.h + + Author: Thomas Prokscha + Date: 25-Feb-2010 + + $Id$ + +***************************************************************************/ + +/*************************************************************************** + * Copyright (C) 2010 by Thomas Prokscha, Paul Scherrer Institut * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program; if not, write to the * + * Free Software Foundation, Inc., * + * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * + ***************************************************************************/ + +#ifndef _PSIMULATEMUTRANSITION_H_ +#define _PSIMULATEMUTRANSITION_H_ + +#include +#include +#include + +// global constants +const Double_t fMuonGyroRatio = 135.54; //!< muon gyromagnetic ratio (MHz/T) +const Double_t fMuonDecayRate = 0.4551; //!< muon decay rate (1/tau_mu, MHz) + +class PSimulateMuTransition : public TObject +{ + public: + PSimulateMuTransition(UInt_t seed = 0); + 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 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 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) + virtual void SetDecayAsymmetry(Double_t value){ fAsymmetry = value; } //!< muon decay asymmetry + virtual void SetMuFraction(Double_t value){ fMuFraction = value; } //!< Muonium fraction + virtual void SetMuFractionState12(Double_t value){ fMuFractionState12 = value; } + virtual void SetMuFractionState23(Double_t value){ fMuFractionState23 = 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 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); + + private: + Bool_t fValid; + TRandom2 *fRandom; + + Double_t fBfield; //!< magnetic field (T) + Double_t fMuPrecFreq12; //!< Mu transition frequency 12 (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 fMuonPrecFreq; //!< muon precession frequency (MHz) + Double_t fCaptureRate; //!< Mu+ electron capture rate (MHz) + Double_t fIonizationRate; //!< Mu0 ionization rate (MHz) + Double_t fSpinFlipRate; //!< Mu0 spin-flip 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; //!< total Mu fraction [0,1] + Double_t fMuFractionState12; //!< fraction of Mu in state 12, 34 + Double_t fMuFractionState23; //!< fraction of Mu in state 23, 14 + 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 Double_t PrecessionPhase(const Double_t &time, const TString chargeState); + virtual Double_t GTFunction(const Double_t &time); //!< transverse field polarization function of Mu0 + virtual Double_t GTSpinFlip(const Double_t &time); //!< transverse field polarization function after spin-flip collisions + virtual void Event(const TString muonString); + + ClassDef(PSimulateMuTransition, 0) +}; + +#endif // _PSIMULATEMUTRANSITION_H_ diff --git a/src/tests/MuTransition/runMuSimulation.C b/src/tests/MuTransition/runMuSimulation.C index 4e2a55d2..433d6f67 100644 --- a/src/tests/MuTransition/runMuSimulation.C +++ b/src/tests/MuTransition/runMuSimulation.C @@ -1,3 +1,155 @@ +<<<<<<< HEAD +/*************************************************************************** + + runMuSimulation.C + + Author: Thomas Prokscha + Date: 25-Feb-2010 + + $Id$ + +***************************************************************************/ + +/*************************************************************************** + * Copyright (C) 2010 by Thomas Prokscha, Paul Scherrer Institut * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program; if not, write to the * + * Free Software Foundation, Inc., * + * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * + ***************************************************************************/ + +void runMuSimulation() +{ + // load library + gSystem->Load("$ROOTSYS/lib/libPSimulateMuTransition"); + + // generate data + TFolder *histosFolder; + TFolder *decayAnaModule; + TFolder *runInfo; + + histosFolder = gROOT->GetRootFolder()->AddFolder("histos", "Histograms"); + gROOT->GetListOfBrowsables()->Add(histosFolder, "histos"); + decayAnaModule = histosFolder->AddFolder("DecayAnaModule", "muSR decay histograms"); + + //prepare to run simulation; here: isotropic Mu in Germanium + UInt_t runNo = 9903; + Double_t T = 300.; //temperature + Double_t capRate = 1.0;//*sqrt(T/200.); + Double_t spinFlipRate = 0.001; + //assume that capture rate varies as sqrt(T), capRate = sigma*v*p , v ~ sqrt(T) + Double_t ionRate; //assume Arrhenius behaviour ionRate = preFac*exp(-EA/kT) + Double_t EA = 100.; //activation energy (meV) + ionRate = 2.9e7 * exp(-EA/(0.08625*T)); // Ge: 2.9*10^7MHz "attempt" frequency; 1K = 0.08625 meV + Double_t B = 100.; //field in G + Double_t Freq12 = 4463; //Mu freq of the 12 transition + Double_t Freq34 = 4463; //Mu freq of the 34 transition + Double_t Freq23 = 4463; //Mu freq of the 23 transition + Double_t Freq14 = 4463; //Mu freq of the 14 transition + Double_t MuFrac = 1.0; //total Mu fraction + Double_t MuFrac12 = 0.5; //Mu in states 12 and 34 + Double_t MuFrac23 = 0.5; //Mu in states 23 and 14 + Int_t Nmuons = 1e7; //number of muons + Double_t Asym = 0.27; //muon decay asymmetry + + // feed run info header + TString tstr; + runInfo = gROOT->GetRootFolder()->AddFolder("RunInfo", "LEM RunInfo"); + gROOT->GetListOfBrowsables()->Add(runInfo, "RunInfo"); + header = new TLemRunHeader(); + tstr = TString("0"); + tstr += runNo; + tstr += TString(" - Mu-frac 1.0, Mu12 -4463MHz (0.5), Mu34 -4463MHz(0.5), T=300K/EA=100meV, Cap. 1.0MHz, 10mT"); + + header->SetRunTitle(tstr.Data()); + header->SetLemSetup("trivial"); + header->SetRunNumber(runNo); + header->SetStartTime(0); + header->SetStopTime(1); + header->SetModeratorHV(32.0, 0.01); + header->SetSampleHV(0.0, 0.01); + header->SetImpEnergy(31.8); + header->SetSampleTemperature(T, 0.001); + header->SetSampleBField(B, 0.1); + header->SetTimeResolution(1.); + header->SetNChannels(12001); + header->SetNHist(2); + header->SetOffsetPPCHistograms(20); + header->SetCuts("none"); + header->SetModerator("none"); + Double_t tt0[2] = {0., 0.}; + header->SetTimeZero(tt0); + runInfo->Add(header); //add header to RunInfo folder + + TH1F *histo[4]; + char str[128]; + for (UInt_t i=0; i<2; i++) { + sprintf(str, "hDecay0%d", (Int_t)i); + histo[i] = new TH1F(str, str, 12001, -0.5, 12000.5); + sprintf(str, "hDecay2%d", (Int_t)i); + histo[i+2] = new TH1F(str, str, 12001, -0.5, 12000.5); + } + + PSimulateMuTransition *simulateMuTransition = new PSimulateMuTransition(); + if (!simulateMuTransition->IsValid()) { + cerr << endl << "**ERROR** while invoking PSimulateTransition" << endl; + return; + } + + simulateMuTransition->SetMuPrecFreq12(Freq12); // MHz + simulateMuTransition->SetMuPrecFreq34(Freq34); // MHz + simulateMuTransition->SetMuPrecFreq23(Freq23); // MHz + simulateMuTransition->SetMuPrecFreq14(Freq14); // MHz + simulateMuTransition->SetMuFraction(MuFrac); // initial Mu fraction + simulateMuTransition->SetMuFractionState12(MuFrac12); // Mu in states 12, 34 + simulateMuTransition->SetMuFractionState23(MuFrac23); // Mu in states 23, 14 + simulateMuTransition->SetBfield(B/10000.); // Tesla + simulateMuTransition->SetCaptureRate(capRate); // MHz + simulateMuTransition->SetIonizationRate(ionRate); // MHz + simulateMuTransition->SetSpinFlipRate(spinFlipRate); // MHz + simulateMuTransition->SetNmuons(Nmuons); + simulateMuTransition->SetDecayAsymmetry(Asym); + simulateMuTransition->SetDebugFlag(kFALSE); // to print time and phase during charge-changing cycle + + simulateMuTransition->PrintSettings(); + + simulateMuTransition->Run(histo[0], histo[1]); + + for (UInt_t i=0; i<4; i++) + decayAnaModule->Add(histo[i]); + + // write file + tstr = TString("0"); + tstr += runNo; + tstr += TString(".root"); + TFile *fout = new TFile(tstr.Data(), "RECREATE", "Midas Fake Histograms"); + if (fout == 0) { + cout << endl << "**ERROR** Couldn't create ROOT file"; + cout << endl << endl; + exit(0); + } + + fout->cd(); + runInfo->Write(); + histosFolder->Write(); + fout->Close(); + cout << "Histograms written to " << tstr.Data() << endl; + delete fout; + + delete [] histo; +} +======= /*************************************************************************** runMuSimulation.C @@ -146,3 +298,4 @@ void runMuSimulation() delete [] histo; } +>>>>>>> 4fec25e423493c58fa21cedec161430f90ffc10d From 96ad763940735abfe33db1a51f02a0182a97966c Mon Sep 17 00:00:00 2001 From: Thomas Prokscha Date: Wed, 17 Feb 2016 14:27:11 +0100 Subject: [PATCH 2/5] Removed faulty lines at the beginning and end of file --- src/tests/MuTransition/runMuSimulation.C | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/tests/MuTransition/runMuSimulation.C b/src/tests/MuTransition/runMuSimulation.C index 433d6f67..4c9f8395 100644 --- a/src/tests/MuTransition/runMuSimulation.C +++ b/src/tests/MuTransition/runMuSimulation.C @@ -1,4 +1,3 @@ -<<<<<<< HEAD /*************************************************************************** runMuSimulation.C @@ -298,4 +297,3 @@ void runMuSimulation() delete [] histo; } ->>>>>>> 4fec25e423493c58fa21cedec161430f90ffc10d From 1c5c94a86b478f2a72c3d40e217c3332170cd3a6 Mon Sep 17 00:00:00 2001 From: Thomas Prokscha Date: Wed, 17 Feb 2016 20:42:13 +0100 Subject: [PATCH 3/5] Use Musr-Root run header --- src/tests/MuTransition/9900.msr | 67 +++++ src/tests/MuTransition/runMuSimulation.C | 337 ++++++++--------------- 2 files changed, 187 insertions(+), 217 deletions(-) create mode 100644 src/tests/MuTransition/9900.msr diff --git a/src/tests/MuTransition/9900.msr b/src/tests/MuTransition/9900.msr new file mode 100644 index 00000000..6b0e8c7c --- /dev/null +++ b/src/tests/MuTransition/9900.msr @@ -0,0 +1,67 @@ +09900- Mu-frac 1.00, Mu12 134.86 MHz(0.27), Mu23 143.71 MHz(0.23), ionRate 608086.30 MHz, capRate 1.00 MHz, SF rate 0.00, 100 G +############################################################### +FITPARAMETER +# Nr. Name Value Step Pos_Error Boundaries + 1 alpha 1.0008 -0.0021 0.0021 0 none + 2 asy 0.2717 -0.0014 0.0014 0 0.33 + 3 phase 1.78 -0.46 0.46 + 4 field 100.418 -0.035 0.035 0 none + 5 rate 0.0000000072 -0.0000000072 0.0013386264 0 100 + 6 asyMu 0 0 none + 7 phaseMu 0 0 none + 8 freqMu 35 0 none + 9 rateMu 0 0 none + +############################################################### +THEORY +asymmetry 2 +TFieldCos 3 fun1 (phase frequency) +simplExpo 5 (rate) ++ +asymmetry 6 +TFieldCos 7 8 (phase frequency) +simplExpo 9 (rate) + +############################################################### +FUNCTIONS +fun1 = par4 * gamma_mu + +############################################################### +RUN 09900 MUE4 PSI MUSR-ROOT (name beamline institute data-file-format) +fittype 2 (asymmetry fit) +alpha 1 +map 0 0 0 0 0 0 0 0 0 0 +forward 1 +backward 2 +backgr.fix 0 0 +data 1 12000 1 12000 +t0 0.0 0.0 +fit 0 8 +packing 5 + +############################################################### +COMMANDS +SET BATCH +MINIMIZE +MINOS +SAVE +END RETURN + +############################################################### +FOURIER +units MHz # units either 'Gauss', 'Tesla', 'MHz', or 'Mc/s' +fourier_power 12 +apodization STRONG # NONE, WEAK, MEDIUM, STRONG +plot POWER # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE +phase 8 +#range_for_phase_correction 50.0 70.0 +range 0 200 + +############################################################### +PLOT 2 (asymmetry plot) +runs 1 +range 0 8 -0.35 0.35 + +############################################################### +STATISTIC --- 2016-02-17 20:38:53 + chisq = 1457.5, NDF = 1595, chisq/NDF = 0.913780 diff --git a/src/tests/MuTransition/runMuSimulation.C b/src/tests/MuTransition/runMuSimulation.C index 4c9f8395..96b3f5d4 100644 --- a/src/tests/MuTransition/runMuSimulation.C +++ b/src/tests/MuTransition/runMuSimulation.C @@ -28,84 +28,62 @@ * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * ***************************************************************************/ +#include "/apps/cern/root-git/include/TMusrRunHeader.h" +#define NDECAYHISTS 2 + void runMuSimulation() { // load library gSystem->Load("$ROOTSYS/lib/libPSimulateMuTransition"); + + // load TMusrRunHeader class if not already done during root startup + if (!TClass::GetDict("TMusrRunHeader")) { + gROOT->LoadMacro("$(ROOTSYS)/lib/libTMusrRunHeader.so"); + } - // generate data + char titleStr[256]; TFolder *histosFolder; TFolder *decayAnaModule; - TFolder *runInfo; - - histosFolder = gROOT->GetRootFolder()->AddFolder("histos", "Histograms"); - gROOT->GetListOfBrowsables()->Add(histosFolder, "histos"); - decayAnaModule = histosFolder->AddFolder("DecayAnaModule", "muSR decay histograms"); - + TFolder *gRunHeader; + TString runTitle; + TString histogramFileName; + TObjArray Slist(0); + TMusrRunPhysicalQuantity prop; + //prepare to run simulation; here: isotropic Mu in Germanium - UInt_t runNo = 9903; - Double_t T = 300.; //temperature - Double_t capRate = 1.0;//*sqrt(T/200.); - Double_t spinFlipRate = 0.001; - //assume that capture rate varies as sqrt(T), capRate = sigma*v*p , v ~ sqrt(T) + UInt_t runNo = 9900; + Double_t T = 300.; //temperature + Double_t spinFlipRate = 0.001; //if spinFlipRate > 0.001 only spin-flip processes will be simulated + Double_t capRate = 1.0;//*sqrt(T/200.); //assume that capture rate varies as sqrt(T), capRate = sigma*v*p , v ~ sqrt(T) Double_t ionRate; //assume Arrhenius behaviour ionRate = preFac*exp(-EA/kT) Double_t EA = 100.; //activation energy (meV) ionRate = 2.9e7 * exp(-EA/(0.08625*T)); // Ge: 2.9*10^7MHz "attempt" frequency; 1K = 0.08625 meV Double_t B = 100.; //field in G - Double_t Freq12 = 4463; //Mu freq of the 12 transition - Double_t Freq34 = 4463; //Mu freq of the 34 transition - Double_t Freq23 = 4463; //Mu freq of the 23 transition - Double_t Freq14 = 4463; //Mu freq of the 14 transition + Double_t Bvar = 0.; //field variance + Double_t Freq12 = 134.858; //Mu freq of the 12 transition + Double_t Freq34 = 4328.142; //Mu freq of the 34 transition + Double_t Freq23 = 143.713; //Mu freq of the 23 transition + Double_t Freq14 = 4606.713; //Mu freq of the 14 transition Double_t MuFrac = 1.0; //total Mu fraction - Double_t MuFrac12 = 0.5; //Mu in states 12 and 34 - Double_t MuFrac23 = 0.5; //Mu in states 23 and 14 - Int_t Nmuons = 1e7; //number of muons + Double_t MuFrac12 = 2*0.266; //Mu in states 12 and 34 + Double_t MuFrac23 = 2*0.234; //Mu in states 23 and 14 + Int_t Nmuons = 1e6; //number of muons Double_t Asym = 0.27; //muon decay asymmetry - // feed run info header - TString tstr; - runInfo = gROOT->GetRootFolder()->AddFolder("RunInfo", "LEM RunInfo"); - gROOT->GetListOfBrowsables()->Add(runInfo, "RunInfo"); - header = new TLemRunHeader(); - tstr = TString("0"); - tstr += runNo; - tstr += TString(" - Mu-frac 1.0, Mu12 -4463MHz (0.5), Mu34 -4463MHz(0.5), T=300K/EA=100meV, Cap. 1.0MHz, 10mT"); + histogramFileName = TString("0"); + histogramFileName += runNo; + histogramFileName += TString(".root"); - header->SetRunTitle(tstr.Data()); - header->SetLemSetup("trivial"); - header->SetRunNumber(runNo); - header->SetStartTime(0); - header->SetStopTime(1); - header->SetModeratorHV(32.0, 0.01); - header->SetSampleHV(0.0, 0.01); - header->SetImpEnergy(31.8); - header->SetSampleTemperature(T, 0.001); - header->SetSampleBField(B, 0.1); - header->SetTimeResolution(1.); - header->SetNChannels(12001); - header->SetNHist(2); - header->SetOffsetPPCHistograms(20); - header->SetCuts("none"); - header->SetModerator("none"); - Double_t tt0[2] = {0., 0.}; - header->SetTimeZero(tt0); - runInfo->Add(header); //add header to RunInfo folder - - TH1F *histo[4]; - char str[128]; - for (UInt_t i=0; i<2; i++) { - sprintf(str, "hDecay0%d", (Int_t)i); - histo[i] = new TH1F(str, str, 12001, -0.5, 12000.5); - sprintf(str, "hDecay2%d", (Int_t)i); - histo[i+2] = new TH1F(str, str, 12001, -0.5, 12000.5); - } + sprintf(titleStr,"- Mu-frac %3.2f, Mu12 %6.2f MHz(%3.2f), Mu23 %6.2f MHz(%3.2f), ionRate %8.2f MHz, capRate %6.2f MHz, SF rate %6.2f MHz, %5.0f G", MuFrac, Freq12, MuFrac12/2, Freq23, MuFrac23/2, ionRate, capRate, spinFlipRate, B); + runTitle = TString("0"); + runTitle += runNo; + runTitle += TString(titleStr); PSimulateMuTransition *simulateMuTransition = new PSimulateMuTransition(); - if (!simulateMuTransition->IsValid()) { + if (!simulateMuTransition->IsValid()){ cerr << endl << "**ERROR** while invoking PSimulateTransition" << endl; return; - } - + } simulateMuTransition->SetMuPrecFreq12(Freq12); // MHz simulateMuTransition->SetMuPrecFreq34(Freq34); // MHz simulateMuTransition->SetMuPrecFreq23(Freq23); // MHz @@ -120,19 +98,92 @@ void runMuSimulation() simulateMuTransition->SetNmuons(Nmuons); simulateMuTransition->SetDecayAsymmetry(Asym); simulateMuTransition->SetDebugFlag(kFALSE); // to print time and phase during charge-changing cycle + + // feed run info header + gRunHeader = gROOT->GetRootFolder()->AddFolder("RunHeader", "MuTransition Simulation Header Info"); + gROOT->GetListOfBrowsables()->Add(gRunHeader, "RunHeader"); +// header = new TLemRunHeader(); + header = new TMusrRunHeader(true); + header->FillFolder(gRunHeader); + gRunHeader->Add(&Slist); + Slist.SetName("RunSummary"); + + header->Set("RunInfo/Generic Validator URL", "http://lmu.web.psi.ch/facilities/software/MusrRoot/validation/MusrRoot.xsd"); + header->Set("RunInfo/Specific Validator URL", "http://lmu.web.psi.ch/facilities/software/MusrRoot/validation/MusrRootLEM.xsd"); + header->Set("RunInfo/Generator", "runMuSimulation"); + + header->Set("RunInfo/File Name", histogramFileName.Data()); + header->Set("RunInfo/Run Title", runTitle.Data()); + header->Set("RunInfo/Run Number", (Int_t) runNo); + header->Set("RunInfo/Run Start Time", "0"); + header->Set("RunInfo/Run Stop Time", "1"); + prop.Set("Run Duration", 1.0, "sec"); + header->Set("RunInfo/Run Duration", prop); + header->Set("RunInfo/Laboratory", "PSI"); + header->Set("RunInfo/Instrument", "MC-Simulation"); + prop.Set("Muon Beam Momentum", 0.0, "MeV/c"); + header->Set("RunInfo/Muon Beam Momentum", prop); + header->Set("RunInfo/Muon Species", "positive muon and muonium"); + header->Set("RunInfo/Muon Source", "Simulation"); + header->Set("RunInfo/Setup", "Monte-Carlo setup"); + header->Set("RunInfo/Comment", "Testing effect of charge-exchange or Mu0 spin flip processes on uSR signal"); + header->Set("RunInfo/Sample Name", "Monte-Carlo"); + header->Set("RunInfo/Sample Temperature", 300); + prop.Set("Sample Magnetic Field", MRH_UNDEFINED, B, Bvar, "G"); + header->Set("RunInfo/Sample Magnetic Field", prop); + header->Set("RunInfo/No of Histos", 2); + prop.Set("Time Resolution", 1.0, "ns", "Simulation"); + header->Set("RunInfo/Time Resolution", prop); + + header->Set("DetectorInfo/Detector001/Name", "hDecay001"); + header->Set("DetectorInfo/Detector001/Histo Number", 1); + header->Set("DetectorInfo/Detector001/Histo Length", 12001); + header->Set("DetectorInfo/Detector001/Time Zero Bin", 0); + header->Set("DetectorInfo/Detector001/First Good Bin", 1); + header->Set("DetectorInfo/Detector001/Last Good Bin", 12001); + + header->Set("DetectorInfo/Detector002/Name", "hDecay002"); + header->Set("DetectorInfo/Detector002/Histo Number", 1); + header->Set("DetectorInfo/Detector002/Histo Length", 12001); + header->Set("DetectorInfo/Detector002/Time Zero Bin", 0); + header->Set("DetectorInfo/Detector002/First Good Bin", 1); + header->Set("DetectorInfo/Detector002/Last Good Bin", 12001); + + // simulation parameters + header->Set("Simulation/Mu Precession frequency 12", Freq12); + header->Set("Simulation/Mu Precession frequency 34", Freq34); + header->Set("Simulation/Mu Precession frequency 23", Freq23); + header->Set("Simulation/Mu Precession frequency 14", Freq14); + header->Set("Simulation/Mu Fraction", MuFrac); + header->Set("Simulation/Mu Fraction 12", MuFrac12); + header->Set("Simulation/Mu Fraction 23", MuFrac23); + header->Set("Simulation/Mu+ 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); + header->Set("Simulation/Decay Asymmetry", Asym); + + histosFolder = gROOT->GetRootFolder()->AddFolder("histos", "Histograms"); + gROOT->GetListOfBrowsables()->Add(histosFolder, "histos"); + decayAnaModule = histosFolder->AddFolder("DecayAnaModule", "muSR decay histograms"); + + TH1F *histo[NDECAYHISTS]; + char str[128]; + for (UInt_t i=0; iAdd(histo[i]); + + // run simulation simulateMuTransition->PrintSettings(); - simulateMuTransition->Run(histo[0], histo[1]); - for (UInt_t i=0; i<4; i++) - decayAnaModule->Add(histo[i]); - // write file - tstr = TString("0"); - tstr += runNo; - tstr += TString(".root"); - TFile *fout = new TFile(tstr.Data(), "RECREATE", "Midas Fake Histograms"); + TFile *fout = new TFile(histogramFileName.Data(), "RECREATE", "Midas MC Histograms"); if (fout == 0) { cout << endl << "**ERROR** Couldn't create ROOT file"; cout << endl << endl; @@ -140,160 +191,12 @@ void runMuSimulation() } fout->cd(); - runInfo->Write(); + header->FillFolder(gRunHeader); + gRunHeader->Write(); histosFolder->Write(); fout->Close(); - cout << "Histograms written to " << tstr.Data() << endl; + cout << "Histograms written to " << histogramFileName.Data() << endl; delete fout; delete [] histo; } -======= -/*************************************************************************** - - runMuSimulation.C - - Author: Thomas Prokscha - Date: 25-Feb-2010 - - $Id$ - -***************************************************************************/ - -/*************************************************************************** - * Copyright (C) 2010 by Thomas Prokscha, Paul Scherrer Institut * - * * - * This program is free software; you can redistribute it and/or modify * - * it under the terms of the GNU General Public License as published by * - * the Free Software Foundation; either version 2 of the License, or * - * (at your option) any later version. * - * * - * This program is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * - * GNU General Public License for more details. * - * * - * You should have received a copy of the GNU General Public License * - * along with this program; if not, write to the * - * Free Software Foundation, Inc., * - * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * - ***************************************************************************/ - -void runMuSimulation() -{ - // load library - gSystem->Load("$ROOTSYS/lib/libPSimulateMuTransition"); - - // generate data - TFolder *histosFolder; - TFolder *decayAnaModule; - TFolder *runInfo; - - histosFolder = gROOT->GetRootFolder()->AddFolder("histos", "Histograms"); - gROOT->GetListOfBrowsables()->Add(histosFolder, "histos"); - decayAnaModule = histosFolder->AddFolder("DecayAnaModule", "muSR decay histograms"); - - //prepare to run simulation; here: isotropic Mu in Germanium - UInt_t runNo = 9903; - Double_t T = 300.; //temperature - Double_t capRate = 1.0;//*sqrt(T/200.); - //assume that capture rate varies as sqrt(T), capRate = sigma*v*p , v ~ sqrt(T) - Double_t ionRate; //assume Arrhenius behaviour ionRate = preFac*exp(-EA/kT) - Double_t EA = 100.; //activation energy (meV) - ionRate = 2.9e7 * exp(-EA/(0.08625*T)); // Ge: 2.9*10^7MHz "attempt" frequency; 1K = 0.08625 meV - Double_t B = 100.; //field in G - Double_t Freq12 = 4463; //Mu freq of the 12 transition - Double_t Freq34 = 4463; //Mu freq of the 34 transition - Double_t Freq23 = 4463; //Mu freq of the 23 transition - Double_t Freq14 = 4463; //Mu freq of the 14 transition - Double_t MuFrac = 1.0; //total Mu fraction - Double_t MuFrac12 = 0.5; //Mu in states 12 and 34 - Double_t MuFrac23 = 0.5; //Mu in states 23 and 14 - Int_t Nmuons = 1e7; //number of muons - Double_t Asym = 0.27; //muon decay asymmetry - - // feed run info header - TString tstr; - runInfo = gROOT->GetRootFolder()->AddFolder("RunInfo", "LEM RunInfo"); - gROOT->GetListOfBrowsables()->Add(runInfo, "RunInfo"); - header = new TLemRunHeader(); - tstr = TString("0"); - tstr += runNo; - tstr += TString(" - Mu-frac 1.0, Mu12 -4463MHz (0.5), Mu34 -4463MHz(0.5), T=300K/EA=100meV, Cap. 1.0MHz, 10mT"); - - header->SetRunTitle(tstr.Data()); - header->SetLemSetup("trivial"); - header->SetRunNumber(runNo); - header->SetStartTime(0); - header->SetStopTime(1); - header->SetModeratorHV(32.0, 0.01); - header->SetSampleHV(0.0, 0.01); - header->SetImpEnergy(31.8); - header->SetSampleTemperature(T, 0.001); - header->SetSampleBField(B, 0.1); - header->SetTimeResolution(1.); - header->SetNChannels(12001); - header->SetNHist(2); - header->SetOffsetPPCHistograms(20); - header->SetCuts("none"); - header->SetModerator("none"); - Double_t tt0[2] = {0., 0.}; - header->SetTimeZero(tt0); - runInfo->Add(header); //add header to RunInfo folder - - TH1F *histo[4]; - char str[128]; - for (UInt_t i=0; i<2; i++) { - sprintf(str, "hDecay0%d", (Int_t)i); - histo[i] = new TH1F(str, str, 12001, -0.5, 12000.5); - sprintf(str, "hDecay2%d", (Int_t)i); - histo[i+2] = new TH1F(str, str, 12001, -0.5, 12000.5); - } - - PSimulateMuTransition *simulateMuTransition = new PSimulateMuTransition(); - if (!simulateMuTransition->IsValid()) { - cerr << endl << "**ERROR** while invoking PSimulateTransition" << endl; - return; - } - - simulateMuTransition->SetMuPrecFreq12(Freq12); // MHz - simulateMuTransition->SetMuPrecFreq34(Freq34); // MHz - simulateMuTransition->SetMuPrecFreq23(Freq23); // MHz - simulateMuTransition->SetMuPrecFreq14(Freq14); // MHz - simulateMuTransition->SetMuFraction(MuFrac); // initial Mu fraction - simulateMuTransition->SetMuFractionState12(MuFrac12); // Mu in states 12, 34 - simulateMuTransition->SetMuFractionState23(MuFrac23); // Mu in states 23, 14 - simulateMuTransition->SetBfield(B/10000.); // Tesla - simulateMuTransition->SetCaptureRate(capRate); // MHz - simulateMuTransition->SetIonizationRate(ionRate); // MHz - simulateMuTransition->SetNmuons(Nmuons); - simulateMuTransition->SetDecayAsymmetry(Asym); - simulateMuTransition->SetDebugFlag(kFALSE); // to print time and phase during charge-changing cycle - - simulateMuTransition->PrintSettings(); - - simulateMuTransition->Run(histo[0], histo[1]); - - for (UInt_t i=0; i<4; i++) - decayAnaModule->Add(histo[i]); - - // write file - tstr = TString("0"); - tstr += runNo; - tstr += TString(".root"); - TFile *fout = new TFile(tstr.Data(), "RECREATE", "Midas Fake Histograms"); - if (fout == 0) { - cout << endl << "**ERROR** Couldn't create ROOT file"; - cout << endl << endl; - exit(0); - } - - fout->cd(); - runInfo->Write(); - histosFolder->Write(); - fout->Close(); - cout << "Histograms written to " << tstr.Data() << endl; - delete fout; - - delete [] histo; -} From 7ae1b9ab5e4e715fff74518e184c58c23c6758aa Mon Sep 17 00:00:00 2001 From: Thomas Prokscha Date: Mon, 22 Feb 2016 22:47:58 +0100 Subject: [PATCH 4/5] changed to complex Mu polarization function --- .../MuTransition/PSimulateMuTransition.cpp | 39 +++++++++++++------ .../MuTransition/PSimulateMuTransition.h | 3 +- 2 files changed, 30 insertions(+), 12 deletions(-) diff --git a/src/tests/MuTransition/PSimulateMuTransition.cpp b/src/tests/MuTransition/PSimulateMuTransition.cpp index 8da00b00..a983b3a5 100644 --- a/src/tests/MuTransition/PSimulateMuTransition.cpp +++ b/src/tests/MuTransition/PSimulateMuTransition.cpp @@ -73,6 +73,7 @@ using namespace std; #include +#include #include "PSimulateMuTransition.h" @@ -246,7 +247,7 @@ Double_t PSimulateMuTransition::PrecessionPhase(const Double_t &time, const TStr if (chargeState == "Mu+") muonPhaseX = TMath::TwoPi()*fMuonPrecFreq*time; else if (chargeState == "Mu0"){ - muoniumPolX = GTFunction(time); + muoniumPolX = GTFunction(time).Re(); muonPhaseX = TMath::ACos(muoniumPolX); } else @@ -263,15 +264,29 @@ Double_t PSimulateMuTransition::PrecessionPhase(const Double_t &time, const TStr * * \param time (us); */ -Double_t PSimulateMuTransition::GTFunction(const Double_t &time) +TComplex PSimulateMuTransition::GTFunction(const Double_t &time) { - Double_t muoniumPolX = 0; + Double_t twoPi = TMath::TwoPi(); + + TComplex complexPol = 0; + complexPol = + 0.5 * fMuFractionState12 * + (TComplex::Exp(TComplex::I()*twoPi*fMuPrecFreq12*time) + + TComplex::Exp(-TComplex::I()*twoPi*fMuPrecFreq34*time)) + + + 0.5 * fMuFractionState23 * + (TComplex::Exp(TComplex::I()*twoPi*fMuPrecFreq23*time) + + TComplex::Exp(TComplex::I()*twoPi*fMuPrecFreq14*time)); + + return complexPol; + +// Double_t muoniumPolX = 0; +// muoniumPolX = 0.5 * +// (fMuFractionState12 * (TMath::Cos(twoPi*fMuPrecFreq12*time) + TMath::Cos(twoPi*fMuPrecFreq34*time)) + +// fMuFractionState23 * (TMath::Cos(twoPi*fMuPrecFreq23*time) + TMath::Cos(twoPi*fMuPrecFreq14*time))); +// +// return muoniumPolX; - muoniumPolX = 0.5 * - (fMuFractionState12 * (TMath::Cos(TMath::TwoPi()*fMuPrecFreq12*time) + TMath::Cos(TMath::TwoPi()*fMuPrecFreq34*time)) + - fMuFractionState23 * (TMath::Cos(TMath::TwoPi()*fMuPrecFreq23*time) + TMath::Cos(TMath::TwoPi()*fMuPrecFreq14*time))); - - return muoniumPolX; } //-------------------------------------------------------------------------- @@ -285,6 +300,7 @@ Double_t PSimulateMuTransition::GTFunction(const Double_t &time) */ Double_t PSimulateMuTransition::GTSpinFlip(const Double_t &time) { + TComplex complexPolX = 1.0; Double_t muoniumPolX = 1.0; //initial polarization in x direction Double_t eventTime = 0; Double_t eventDiffTime = 0; @@ -292,18 +308,19 @@ Double_t PSimulateMuTransition::GTSpinFlip(const Double_t &time) eventTime += NextEventTime(fSpinFlipRate); if (eventTime >= time){ - muoniumPolX = GTFunction(time); + muoniumPolX = GTFunction(time).Re(); } else{ while (eventTime < time){ eventDiffTime = eventTime - lastEventTime; - muoniumPolX = muoniumPolX * GTFunction(eventDiffTime); + complexPolX = complexPolX * GTFunction(eventDiffTime); lastEventTime = eventTime; eventTime += NextEventTime(fSpinFlipRate); } // calculate for the last collision eventDiffTime = time - lastEventTime; - muoniumPolX = muoniumPolX * GTFunction(eventDiffTime); + complexPolX = complexPolX * GTFunction(eventDiffTime); + muoniumPolX = complexPolX.Re(); } return muoniumPolX; diff --git a/src/tests/MuTransition/PSimulateMuTransition.h b/src/tests/MuTransition/PSimulateMuTransition.h index 4caac1c1..1a8d0f31 100644 --- a/src/tests/MuTransition/PSimulateMuTransition.h +++ b/src/tests/MuTransition/PSimulateMuTransition.h @@ -34,6 +34,7 @@ #include #include #include +#include // global constants const Double_t fMuonGyroRatio = 135.54; //!< muon gyromagnetic ratio (MHz/T) @@ -95,7 +96,7 @@ class PSimulateMuTransition : public TObject virtual Double_t NextEventTime(const Double_t &EventRate); // virtual Double_t PrecessionPhase(const Double_t &time, const Double_t &frequency); virtual Double_t PrecessionPhase(const Double_t &time, const TString chargeState); - virtual Double_t GTFunction(const Double_t &time); //!< transverse field polarization function of Mu0 + virtual TComplex GTFunction(const Double_t &time); //!< transverse field polarization function of Mu0 virtual Double_t GTSpinFlip(const Double_t &time); //!< transverse field polarization function after spin-flip collisions virtual void Event(const TString muonString); From fe48da74a14963c3dc98709c5136d8cdab904691 Mon Sep 17 00:00:00 2001 From: Suter Andreas Date: Tue, 23 Feb 2016 13:22:41 +0100 Subject: [PATCH 5/5] It is now possible to export the averaged data/Fourier --- ChangeLog | 2 +- src/classes/PMusrCanvas.cpp | 754 +++++++++++++----------------------- src/include/PMusrCanvas.h | 5 +- 3 files changed, 274 insertions(+), 487 deletions(-) diff --git a/ChangeLog b/ChangeLog index dd3ab127..4bbe7137 100644 --- a/ChangeLog +++ b/ChangeLog @@ -4,7 +4,7 @@ changes since 0.17.0 =================================== - +NEW 2016-02-23 It is now possible to export the averaged data/Fourier changes since 0.16.0 =================================== diff --git a/src/classes/PMusrCanvas.cpp b/src/classes/PMusrCanvas.cpp index ec99e93a..54bcc51c 100644 --- a/src/classes/PMusrCanvas.cpp +++ b/src/classes/PMusrCanvas.cpp @@ -1516,7 +1516,6 @@ void PMusrCanvas::ExportData(const Char_t *fileName) Int_t xmaxBin; Double_t xmin; Double_t xmax; - Double_t time, freq; Double_t xval, yval; switch (fPlotType) { @@ -1533,29 +1532,12 @@ void PMusrCanvas::ExportData(const Char_t *fileName) xmax = fHistoFrame->GetXaxis()->GetBinCenter(xmaxBin); // fill ascii dump data - for (UInt_t i=0; iGetNbinsX(); j++) { - // get time - time = fData[i].diff->GetBinCenter(j); - // check if time is in the current range - if ((time >= xmin) && (time <= xmax)) { - dump.dataX.push_back(time); - dump.data.push_back(fData[i].diff->GetBinContent(j)); - dump.dataErr.push_back(fData[i].diff->GetBinError(j)); - } + if (fAveragedView) { + GetExportDataSet(fDataAvg.diff, xmin, xmax, dumpVector); + } else { // go through all the histogramms + for (UInt_t i=0; i 0) - dumpVector.push_back(dump); } break; case PV_FOURIER_REAL: @@ -1566,28 +1548,12 @@ void PMusrCanvas::ExportData(const Char_t *fileName) xmax = fData[0].diffFourierRe->GetXaxis()->GetBinCenter(xmaxBin); // fill ascii dump data - for (UInt_t i=0; iGetNbinsX(); j++) { - // get frequency - freq = fData[i].diffFourierRe->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.dataX.push_back(freq); - dump.data.push_back(fData[i].diffFourierRe->GetBinContent(j)); - } + if (fAveragedView) { + GetExportDataSet(fDataAvg.diffFourierRe, xmin, xmax, dumpVector, false); + } else { // go through all the histogramms + for (UInt_t i=0; i 0) - dumpVector.push_back(dump); } break; case PV_FOURIER_IMAG: @@ -1598,28 +1564,12 @@ void PMusrCanvas::ExportData(const Char_t *fileName) xmax = fData[0].diffFourierIm->GetXaxis()->GetBinCenter(xmaxBin); // fill ascii dump data - for (UInt_t i=0; iGetNbinsX(); j++) { - // get frequency - freq = fData[i].diffFourierIm->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.dataX.push_back(freq); - dump.data.push_back(fData[i].diffFourierIm->GetBinContent(j)); - } + if (fAveragedView) { + GetExportDataSet(fDataAvg.diffFourierIm, xmin, xmax, dumpVector, false); + } else { // go through all the histogramms + for (UInt_t i=0; i 0) - dumpVector.push_back(dump); } break; case PV_FOURIER_REAL_AND_IMAG: @@ -1630,37 +1580,14 @@ void PMusrCanvas::ExportData(const Char_t *fileName) xmax = fData[0].diffFourierRe->GetXaxis()->GetBinCenter(xmaxBin); // fill ascii dump data - for (UInt_t i=0; iGetNbinsX(); j++) { - // get frequency - freq = fData[i].diffFourierRe->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.dataX.push_back(freq); - dump.data.push_back(fData[i].diffFourierRe->GetBinContent(j)); - } + if (fAveragedView) { + GetExportDataSet(fDataAvg.diffFourierRe, xmin, xmax, dumpVector, false); + GetExportDataSet(fDataAvg.diffFourierIm, xmin, xmax, dumpVector, false); + } else { // go through all the histogramms + for (UInt_t i=0; iGetNbinsX(); j++) { - // get frequency - freq = fData[i].diffFourierIm->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.dataX.push_back(freq); - dump.data.push_back(fData[i].diffFourierIm->GetBinContent(j)); - } - } - - // if anything found keep it - if (dump.dataX.size() > 0) - dumpVector.push_back(dump); } break; case PV_FOURIER_PWR: @@ -1671,28 +1598,12 @@ void PMusrCanvas::ExportData(const Char_t *fileName) xmax = fData[0].diffFourierPwr->GetXaxis()->GetBinCenter(xmaxBin); // fill ascii dump data - for (UInt_t i=0; iGetNbinsX(); j++) { - // get frequency - freq = fData[i].diffFourierPwr->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.dataX.push_back(freq); - dump.data.push_back(fData[i].diffFourierPwr->GetBinContent(j)); - } + if (fAveragedView) { + GetExportDataSet(fDataAvg.diffFourierPwr, xmin, xmax, dumpVector, false); + } else { // go through all the histogramms + for (UInt_t i=0; i 0) - dumpVector.push_back(dump); } break; case PV_FOURIER_PHASE: @@ -1703,28 +1614,12 @@ void PMusrCanvas::ExportData(const Char_t *fileName) xmax = fData[0].diffFourierPhase->GetXaxis()->GetBinCenter(xmaxBin); // fill ascii dump data - for (UInt_t i=0; iGetNbinsX(); j++) { - // get frequency - freq = fData[i].diffFourierPhase->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.dataX.push_back(freq); - dump.data.push_back(fData[i].diffFourierPhase->GetBinContent(j)); - } + if (fAveragedView) { + GetExportDataSet(fDataAvg.diffFourierPhase, xmin, xmax, dumpVector, false); + } else { // go through all the histogramms + for (UInt_t i=0; i 0) - dumpVector.push_back(dump); } break; default: @@ -1740,40 +1635,14 @@ void PMusrCanvas::ExportData(const Char_t *fileName) xmax = fHistoFrame->GetXaxis()->GetBinCenter(xmaxBin); // fill ascii dump data - for (UInt_t i=0; iGetNbinsX(); j++) { - // get time - time = fData[i].data->GetBinCenter(j); - // check if time is in the current range - if ((time >= xmin) && (time <= xmax)) { - dump.dataX.push_back(time); - dump.data.push_back(fData[i].data->GetBinContent(j)); - dump.dataErr.push_back(fData[i].data->GetBinError(j)); - } + if (fAveragedView) { + GetExportDataSet(fDataAvg.data, xmin, xmax, dumpVector); + GetExportDataSet(fDataAvg.theory, xmin, xmax, dumpVector, false); + } else { // go through all the histogramms + for (UInt_t i=0; iGetNbinsX(); j++) { - // get time - time = fData[i].theory->GetBinCenter(j); - // check if time is in the current range - if ((time >= xmin) && (time <= xmax)) { - dump.theoryX.push_back(time); - dump.theory.push_back(fData[i].theory->GetBinContent(j)); - } - } - - // if anything found keep it - if (dump.dataX.size() > 0) - dumpVector.push_back(dump); } break; @@ -1785,39 +1654,14 @@ void PMusrCanvas::ExportData(const Char_t *fileName) xmax = fData[0].dataFourierRe->GetXaxis()->GetBinCenter(xmaxBin); // fill ascii dump data - for (UInt_t i=0; iGetNbinsX(); j++) { - // get frequency - freq = fData[i].dataFourierRe->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.dataX.push_back(freq); - dump.data.push_back(fData[i].dataFourierRe->GetBinContent(j)); - } + if (fAveragedView) { + GetExportDataSet(fDataAvg.dataFourierRe, xmin, xmax, dumpVector, false); + GetExportDataSet(fDataAvg.theoryFourierRe, xmin, xmax, dumpVector, false); + } else { // go through all the histogramms + for (UInt_t i=0; iGetNbinsX(); j++) { - // get frequency - freq = fData[i].theoryFourierRe->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.theoryX.push_back(freq); - dump.theory.push_back(fData[i].theoryFourierRe->GetBinContent(j)); - } - } - - // if anything found keep it - if (dump.dataX.size() > 0) - dumpVector.push_back(dump); } break; case PV_FOURIER_IMAG: @@ -1828,39 +1672,14 @@ void PMusrCanvas::ExportData(const Char_t *fileName) xmax = fData[0].dataFourierIm->GetXaxis()->GetBinCenter(xmaxBin); // fill ascii dump data - for (UInt_t i=0; iGetNbinsX(); j++) { - // get frequency - freq = fData[i].dataFourierIm->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.dataX.push_back(freq); - dump.data.push_back(fData[i].dataFourierIm->GetBinContent(j)); - } + if (fAveragedView) { + GetExportDataSet(fDataAvg.dataFourierIm, xmin, xmax, dumpVector, false); + GetExportDataSet(fDataAvg.theoryFourierIm, xmin, xmax, dumpVector, false); + } else { // go through all the histogramms + for (UInt_t i=0; iGetNbinsX(); j++) { - // get frequency - freq = fData[i].theoryFourierIm->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.theoryX.push_back(freq); - dump.theory.push_back(fData[i].theoryFourierIm->GetBinContent(j)); - } - } - - // if anything found keep it - if (dump.dataX.size() > 0) - dumpVector.push_back(dump); } break; case PV_FOURIER_REAL_AND_IMAG: @@ -1870,80 +1689,18 @@ void PMusrCanvas::ExportData(const Char_t *fileName) xmin = fData[0].dataFourierRe->GetXaxis()->GetBinCenter(xminBin); xmax = fData[0].dataFourierRe->GetXaxis()->GetBinCenter(xmaxBin); - // fill ascii dump data - for (UInt_t i=0; iGetNbinsX(); j++) { - // get frequency - freq = fData[i].dataFourierRe->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.dataX.push_back(freq); - dump.data.push_back(fData[i].dataFourierRe->GetBinContent(j)); - } + if (fAveragedView) { + GetExportDataSet(fDataAvg.dataFourierRe, xmin, xmax, dumpVector, false); + GetExportDataSet(fDataAvg.theoryFourierRe, xmin, xmax, dumpVector, false); + GetExportDataSet(fDataAvg.dataFourierIm, xmin, xmax, dumpVector, false); + GetExportDataSet(fDataAvg.theoryFourierIm, xmin, xmax, dumpVector, false); + } else { // go through all the histogramms + for (UInt_t i=0; iGetNbinsX(); j++) { - // get frequency - freq = fData[i].theoryFourierRe->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.theoryX.push_back(freq); - dump.theory.push_back(fData[i].theoryFourierRe->GetBinContent(j)); - } - } - - // if anything found keep it - if (dump.dataX.size() > 0) - dumpVector.push_back(dump); - - //----------------------------- - // Im - //----------------------------- - // clean up dump - dump.dataX.clear(); - dump.data.clear(); - dump.dataErr.clear(); - dump.theoryX.clear(); - dump.theory.clear(); - - // go through all data bins - for (Int_t j=1; jGetNbinsX(); j++) { - // get frequency - freq = fData[i].dataFourierIm->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.dataX.push_back(freq); - dump.data.push_back(fData[i].dataFourierIm->GetBinContent(j)); - } - } - - // go through all theory bins - for (Int_t j=1; jGetNbinsX(); j++) { - // get frequency - freq = fData[i].theoryFourierIm->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.theoryX.push_back(freq); - dump.theory.push_back(fData[i].theoryFourierIm->GetBinContent(j)); - } - } - - // if anything found keep it - if (dump.dataX.size() > 0) - dumpVector.push_back(dump); - } break; case PV_FOURIER_PWR: @@ -1954,39 +1711,14 @@ void PMusrCanvas::ExportData(const Char_t *fileName) xmax = fData[0].dataFourierPwr->GetXaxis()->GetBinCenter(xmaxBin); // fill ascii dump data - for (UInt_t i=0; iGetNbinsX(); j++) { - // get frequency - freq = fData[i].dataFourierPwr->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.dataX.push_back(freq); - dump.data.push_back(fData[i].dataFourierPwr->GetBinContent(j)); - } + if (fAveragedView) { + GetExportDataSet(fDataAvg.dataFourierPwr, xmin, xmax, dumpVector, false); + GetExportDataSet(fDataAvg.theoryFourierPwr, xmin, xmax, dumpVector, false); + } else { // go through all the histogramms + for (UInt_t i=0; iGetNbinsX(); j++) { - // get frequency - freq = fData[i].theoryFourierPwr->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.theoryX.push_back(freq); - dump.theory.push_back(fData[i].theoryFourierPwr->GetBinContent(j)); - } - } - - // if anything found keep it - if (dump.dataX.size() > 0) - dumpVector.push_back(dump); } break; case PV_FOURIER_PHASE: @@ -1997,39 +1729,14 @@ void PMusrCanvas::ExportData(const Char_t *fileName) xmax = fData[0].dataFourierPhase->GetXaxis()->GetBinCenter(xmaxBin); // fill ascii dump data - for (UInt_t i=0; iGetNbinsX(); j++) { - // get frequency - freq = fData[i].dataFourierPhase->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.dataX.push_back(freq); - dump.data.push_back(fData[i].dataFourierPhase->GetBinContent(j)); - } + if (fAveragedView) { + GetExportDataSet(fDataAvg.dataFourierPhase, xmin, xmax, dumpVector, false); + GetExportDataSet(fDataAvg.theoryFourierPhase, xmin, xmax, dumpVector, false); + } else { // go through all the histogramms + for (UInt_t i=0; iGetNbinsX(); j++) { - // get frequency - freq = fData[i].theoryFourierPhase->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.theoryX.push_back(freq); - dump.theory.push_back(fData[i].theoryFourierPhase->GetBinContent(j)); - } - } - - // if anything found keep it - if (dump.dataX.size() > 0) - dumpVector.push_back(dump); } break; default: @@ -2053,8 +1760,6 @@ void PMusrCanvas::ExportData(const Char_t *fileName) dump.dataX.clear(); dump.data.clear(); dump.dataErr.clear(); - dump.theoryX.clear(); - dump.theory.clear(); // go through all data bins for (Int_t j=0; jGetN(); j++) { @@ -2102,8 +1807,6 @@ void PMusrCanvas::ExportData(const Char_t *fileName) dump.dataX.clear(); dump.data.clear(); dump.dataErr.clear(); - dump.theoryX.clear(); - dump.theory.clear(); // go through all data bins for (Int_t j=0; jGetN(); j++) { @@ -2123,8 +1826,8 @@ void PMusrCanvas::ExportData(const Char_t *fileName) fNonMusrData[i].theory->GetPoint(j,xval,yval); // check if time is in the current range if ((xval >= xmin) && (xval <= xmax)) { - dump.theoryX.push_back(xval); - dump.theory.push_back(yval); + dump.dataX.push_back(xval); + dump.data.push_back(yval); } } @@ -2164,72 +1867,95 @@ void PMusrCanvas::ExportData(const Char_t *fileName) } // find out what is the longest data/theory vector - UInt_t maxDataLength = 0; - UInt_t maxTheoryLength = 0; + UInt_t maxLength = 0; for (UInt_t i=0; i maxTheoryLength) - maxLength = maxDataLength; - else - maxLength = maxTheoryLength; - // write data and theory for (UInt_t i=0; i extract data for export. + * + * \param data + * \param xmin + * \param xmax + * \param dumpData + * \param hasError + */ +void PMusrCanvas::GetExportDataSet(const TH1F *data, const Double_t xmin, const Double_t xmax, + PMusrCanvasAsciiDumpVector &dumpData, const Bool_t hasError) +{ + PMusrCanvasAsciiDump dump; + Double_t x=0.0; + + // go through all difference data bins + for (Int_t j=1; jGetNbinsX(); j++) { + // get time/freq + x = data->GetBinCenter(j); + // check if x is in the current range + if ((x >= xmin) && (x <= xmax)) { + dump.dataX.push_back(x); + dump.data.push_back(data->GetBinContent(j)); + if (hasError) + dump.dataErr.push_back(data->GetBinError(j)); + } + } + + // if anything found keep it + if (dump.dataX.size() > 0) + dumpData.push_back(dump); +} + //-------------------------------------------------------------------------- // CreateStyle (private) //-------------------------------------------------------------------------- diff --git a/src/include/PMusrCanvas.h b/src/include/PMusrCanvas.h index 735af4f4..42e73dcc 100644 --- a/src/include/PMusrCanvas.h +++ b/src/include/PMusrCanvas.h @@ -181,8 +181,6 @@ typedef struct { PDoubleVector dataX; ///< x-axis data set PDoubleVector data; ///< y-axis data set PDoubleVector dataErr; ///< error of the y-axis data set - PDoubleVector theoryX; ///< x-axis theory set - PDoubleVector theory; ///< y-axis theory set } PMusrCanvasAsciiDump; //------------------------------------------------------------------------ @@ -333,6 +331,9 @@ class PMusrCanvas : public TObject, public TQObject virtual Double_t GetInterpolatedValue(TH1F* histo, Double_t xVal); + virtual void GetExportDataSet(const TH1F *data, const Double_t xmin, const Double_t xmax, + PMusrCanvasAsciiDumpVector &dumpData, const Bool_t hasError=true); + ClassDef(PMusrCanvas, 1) };