From 341fa9f77cad94980a0ed0f7f0a949e44141d17d Mon Sep 17 00:00:00 2001 From: prokscha Date: Tue, 20 Jul 2021 10:51:10 +0200 Subject: [PATCH] Removed MuTransistion from tests directory. --- src/tests/MuTransition/7700.msr | 78 --- src/tests/MuTransition/9102.msr | 67 --- src/tests/MuTransition/9900.msr | 67 --- .../Makefile.PSimulateMuTransition | 105 ---- .../MuTransition/PSimulateMuTransition.cpp | 452 ------------------ .../MuTransition/PSimulateMuTransition.h | 119 ----- .../PSimulateMuTransitionLinkDef.h | 40 -- src/tests/MuTransition/runMuSimulation.C | 270 ----------- src/tests/MuTransition/testAnalysis.C | 59 --- 9 files changed, 1257 deletions(-) delete mode 100644 src/tests/MuTransition/7700.msr delete mode 100644 src/tests/MuTransition/9102.msr delete mode 100644 src/tests/MuTransition/9900.msr delete mode 100644 src/tests/MuTransition/Makefile.PSimulateMuTransition delete mode 100644 src/tests/MuTransition/PSimulateMuTransition.cpp delete mode 100644 src/tests/MuTransition/PSimulateMuTransition.h delete mode 100644 src/tests/MuTransition/PSimulateMuTransitionLinkDef.h delete mode 100644 src/tests/MuTransition/runMuSimulation.C delete mode 100644 src/tests/MuTransition/testAnalysis.C diff --git a/src/tests/MuTransition/7700.msr b/src/tests/MuTransition/7700.msr deleted file mode 100644 index ee61c73f..00000000 --- a/src/tests/MuTransition/7700.msr +++ /dev/null @@ -1,78 +0,0 @@ -07700- complexMuPol, A0 100MHz, Mu-frac 1.00, Mu12 40.43 MHz(0.49), Mu23 256.25 MHz(0.01), ionRate 0.000 MHz, capRate 0.000 MHz, SF rate 0.010 MHz, 106.5 G -############################################################### -FITPARAMETER -# Nr. Name Value Step Pos_Error Boundaries - 1 alpha 0.99961 -0.00091 0.00091 0 none - 2 asy 0 0 none 0 0.33 - 3 phase 0 0 none - 4 field 106.5 0 none 0 none - 5 rate 0 0 none 0 150 - 6 asyMu12 0.135 0 none - 7 phaseMu12 0.68 -0.42 0.42 - 8 freqMu12 40.43248 -0.00045 0.00045 - 9 rateMu12 0.0104 -0.0019 0.0019 0 100 - 10 asyMu34 0.135 0 none - 11 phaseMu34 -0.01 -0.43 0.43 - 12 freqMu34 59.56737 -0.00046 0.00046 - 13 rateMu34 0.0164 -0.0019 0.0019 0 100 - -############################################################### -THEORY -asymmetry 2 -TFieldCos 3 fun1 (phase frequency) -simplExpo 5 (rate) -+ -asymmetry 6 -TFieldCos 7 8 (phase frequency) -simplExpo 9 (rate) -+ -asymmetry 10 -TFieldCos 11 12 (phase frequency) -simplExpo 13 (rate) - -############################################################### -FUNCTIONS -fun1 = par4 * gamma_mu -fun2 = par8 * gamma_mu - -############################################################### -GLOBAL -fittype 2 (asymmetry fit) -data 1 18000 1 18000 -fit 0.005 8 -packing 1 - -############################################################### -RUN 07700 MUE4 PSI MUSR-ROOT (name beamline institute data-file-format) -#ADDRUN 07701 MUE4 PSI MUSR-ROOT (name beamline institute data-file-format) -alpha 1 -map 0 0 0 0 0 0 0 0 0 0 -forward 1 -backward 2 -backgr.fix 0.000000 0.000000 - -############################################################### -COMMANDS -MINIMIZE -MINOS -SAVE - -############################################################### -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.000000 -#range_for_phase_correction 50.0 70.0 -#range 0 200 - -############################################################### -PLOT 2 (asymmetry plot) -runs 1 -range 0 8 -view_packing 2 - -############################################################### -STATISTIC --- 2016-03-03 17:25:40 - chisq = 8091.8, NDF = 7988, chisq/NDF = 1.012996 diff --git a/src/tests/MuTransition/9102.msr b/src/tests/MuTransition/9102.msr deleted file mode 100644 index bde974af..00000000 --- a/src/tests/MuTransition/9102.msr +++ /dev/null @@ -1,67 +0,0 @@ -09102, Mu-frac 0.5, Mu 41MHz (0.42), Mu -35MHz(0.32), Ion. 250MHz, Capt. 1.5MHz, 0.1kG, -############################################################### -FITPARAMETER -# Nr. Name Value Step Pos_Error Boundaries - 1 alpha 0.9998 0.000655606 none 0 none - 2 asy 0.222118 0.00114398 none 0 0.33 - 3 phase 1.95348 0.308634 none - 4 field 100.706 0.0663548 none 0 none - 5 rate 0.526345 0.00550619 none 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 09102 MUE4 PSI ROOT-NPP (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 1 1 -fit 0.00 8.00 -packing 5 - -############################################################### -COMMANDS -SET BATCH -MINIMIZE -#MINOS -SAVE -END RETURN - -############################################################### -FOURIER -units MHz # units either 'Gauss', 'MHz', or 'Mc/s' -fourier_power 12 -apodization STRONG # NONE, WEAK, MEDIUM, STRONG -plot POWER # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE -phase 8.50 -#range_for_phase_correction 50.0 70.0 -range 0.00 200.00 - -############################################################### -PLOT 2 (asymmetry plot) -runs 1 -range 0.00 8.00 -0.30 0.30 - -############################################################### -STATISTIC --- 2010-03-06 18:19:16 - chisq = 1666.28948, NDF = 1595, chisq/NDF = 1.0446956 diff --git a/src/tests/MuTransition/9900.msr b/src/tests/MuTransition/9900.msr deleted file mode 100644 index 6b0e8c7c..00000000 --- a/src/tests/MuTransition/9900.msr +++ /dev/null @@ -1,67 +0,0 @@ -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/Makefile.PSimulateMuTransition b/src/tests/MuTransition/Makefile.PSimulateMuTransition deleted file mode 100644 index 93d1a9ae..00000000 --- a/src/tests/MuTransition/Makefile.PSimulateMuTransition +++ /dev/null @@ -1,105 +0,0 @@ -#--------------------------------------------------- -# Makefile.PSimulateMuTransition -# -# Author: Thomas Prokscha -# Date: 25-Feb-2010 -# -# $Id:$ -# -#--------------------------------------------------- - -#--------------------------------------------------- -# get compilation and library flags from root-config - -ROOTCFLAGS = $(shell $(ROOTSYS)/bin/root-config --cflags) -ROOTLIBS = $(shell $(ROOTSYS)/bin/root-config --libs) -ROOTGLIBS = $(shell $(ROOTSYS)/bin/root-config --glibs) - -#--------------------------------------------------- -# depending on the architecture, choose the compiler, -# linker, and the flags to use -# - -OSTYPE = $(shell uname) - -ifeq ($(OSTYPE),Linux) -OS = LINUX -endif -ifeq ($(OSTYPE),Linux-gnu) -OS = LINUX -endif -ifeq ($(OSTYPE),darwin) -OS = DARWIN -endif - -# -- Linux -ifeq ($(OS),LINUX) -CXX = g++ -CXXFLAGS = -Wall -Wno-trigraphs -fPIC -INCLUDES = -I../include -LD = g++ -LDFLAGS = -g -SOFLAGS = -O -shared -endif - -# -- Darwin -ifeq ($(OS),DARWIN) -CXX = g++ -CXXFLAGS = -Wall -Wno-trigraphs -fPIC -INCLUDES = -I../include -LD = g++ -LDFLAGS = -g -SOFLAGS = -dynamic -endif - -# installation directory -INSTALL_DIR = /usr/local - -# the output from the root-config script: -CXXFLAGS += $(ROOTCFLAGS) -LDFLAGS += - -# the ROOT libraries (G = graphic) -LIBS = $(ROOTLIBS) -lXMLParser -GLIBS = $(ROOTGLIBS) -lXMLParser - -# some definitions: headers (used to generate *Dict* stuff), sources, objects,... -OBJS = -OBJS += PSimulateMuTransition.o PSimulateMuTransitionDict.o - -SHLIB = libPSimulateMuTransition.so - -# make the shared lib: -# -all: $(SHLIB) - -$(SHLIB): $(OBJS) - @echo "---> Building shared library $(SHLIB) ..." - /bin/rm -f $(SHLIB) - $(LD) $(OBJS) $(SOFLAGS) -o $(SHLIB) $(LIBS) - @echo "done" - -# clean up: remove all object file (and core files) -# semicolon needed to tell make there is no source -# for this target! -# -clean:; @rm -f $(OBJS) *Dict* core* - @echo "---> removing $(OBJS)" - @echo "---> removing *~" - @rm -f *~ -# -$(OBJS): %.o: %.cpp - $(CXX) $(INCLUDES) $(CXXFLAGS) -c $< - -PSimulateMuTransitionDict.cpp: PSimulateMuTransition.h PSimulateMuTransitionLinkDef.h - @echo "Generating dictionary $@..." - rootcint -f $@ -c -p $^ - -install: all - @echo "Installing shared lib: libPSimulateMuTransition.so" -ifeq ($(OS),LINUX) - cp -pv $(SHLIB) $(INSTALL_DIR)/lib - cp -pv PSimulateMuTransition.h $(INSTALL_DIR)/include -# for root6 - cp -pv PSimulateMuTransitionDict_rdict.pcm $(INSTALL_DIR)/lib -endif diff --git a/src/tests/MuTransition/PSimulateMuTransition.cpp b/src/tests/MuTransition/PSimulateMuTransition.cpp deleted file mode 100644 index 3b453aea..00000000 --- a/src/tests/MuTransition/PSimulateMuTransition.cpp +++ /dev/null @@ -1,452 +0,0 @@ -/*************************************************************************** - - PSimulateMuTransition.cpp - - Author: Thomas Prokscha - Date: 25-Feb-2010, 14-Apr-2016 - - 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 polarization under successive Mu+/Mu0 charge-exchange - or Mu0 spin-flip 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 complex expression of - equation (4) in the paper of M. Senba, J. Phys. B 23, 1545 (1990), or - equation (7) in the paper of M. Senba, J. Phys. B 24, 3531 (1991); - 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) Mu+ electron-capture rate - 5) Mu0 ionization rate - 6) Mu0 spin-flip rate - 7) initial muon spin phase - 9) total muon decay asymmetry - 9) number of muon decays to be generated. - 10) 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 total - muon spin polarization Px(t_i)*Px(t_c). - 4) get the next electron capture time, continue until t_d is reached, and calculate - the resulting polarization. - - The Mu0 spin-flip processes are calculated in GTSpinFlip(), using eq. (17) of - M. Senba, J. Phys. B 24, 3531 (1991). - -***************************************************************************/ - -/*************************************************************************** - * 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 - -#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 - 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) - 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.25; - fMuFractionState34 = 0.25; - fMuFractionState23 = 0.25; - fMuFractionState14 = 0.25; - fMuFractionState13 = 0.; - fMuFractionState24 = 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 << "Mu0 precession frequency 12 (MHz) = " << fMuPrecFreq12; - 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; - cout << endl << "Mu0 ionizatioan rate (MHz) = " << fIonizationRate; - cout << endl << "Mu0 spin-flip rate (MHz) = " << fSpinFlipRate; - if (fSpinFlipRate > 0.001) - cout << endl << "!!! Note: spin-flip rate > 0.001 only spin-flip processes are considered!!!"; - else{ - cout << endl << "!!! spin-flip rate <= 0.001: only charge-exchange cycles are considered!!!"; - cout << endl << "!!! if spin-flip rate > 0.001, only spin-flip processes are considered!!!"; - } - cout << endl << "Decay asymmetry = " << fAsymmetry; - cout << endl << "Muonium fraction = " << fMuFraction; - cout << endl << "Muonium fraction state12 = " << fMuFractionState12; - cout << endl << "Muonium fraction state34 = " << fMuFractionState34; - cout << endl << "Muonium fraction state23 = " << fMuFractionState23; - cout << endl << "Muonium fraction state14 = " << fMuFractionState14; - cout << endl << "Muonium fraction state13 = " << fMuFractionState13; - cout << endl << "Muonium fraction state24 = " << fMuFractionState24; - 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; -} - -//-------------------------------------------------------------------------- -// 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) -{ -// Double_t muoniumPolX = 1.0; //polarization in x direction - Int_t i; - if (histoForward == 0 || histoBackward == 0) - return; - - fMuonPrecFreq = fMuonGyroRatio * fBfield; - - 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) - fMuonPhase += TMath::ACos(Event("Mu+")); - else - fMuonPhase += TMath::ACos(Event("Mu0")); - } - // 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%fNshowProgress) == 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).Re(); -// if (fDebugFlag) cout << "muoniumPolX = " << muoniumPolX << endl; -// 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); - */ -TComplex PSimulateMuTransition::GTFunction(const Double_t &time, const TString chargeState) -{ - Double_t twoPi = TMath::TwoPi(); - - TComplex complexPol = 0; - - if (chargeState == "Mu+") - complexPol = TComplex::Exp(-TComplex::I()*twoPi*fMuonPrecFreq*time); - else{ - complexPol = - (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)) - + - (fMuFractionState13 * TComplex::Exp(+TComplex::I()*twoPi*fMuPrecFreq13*time) + - fMuFractionState24 * TComplex::Exp(+TComplex::I()*twoPi*fMuPrecFreq24*time)); - } - - return complexPol; -} - -//-------------------------------------------------------------------------- -// 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) -{ - TComplex complexPolX = 1.0; - 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, "Mu0").Re(); - } - else{ - while (eventTime < time){ - eventDiffTime = eventTime - lastEventTime; - complexPolX = complexPolX * GTFunction(eventDiffTime, "Mu0"); - lastEventTime = eventTime; - eventTime += NextEventTime(fSpinFlipRate); - } - // calculate for the last collision - eventDiffTime = time - lastEventTime; - complexPolX = complexPolX * GTFunction(eventDiffTime, "Mu0"); - muoniumPolX = complexPolX.Re(); - } - - return muoniumPolX; -} - -//-------------------------------------------------------------------------- -// Event (private) -//-------------------------------------------------------------------------- -/** - *

Generates "muon event": simulate muon spin polarization 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. Use complex polarization - * functions. - * 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, Px(t_i); else calculate spin precession for t_c. - * 3) Determine next ionization time t_i+1; calculate Px(t_i+1) in Muonium. Polarization - * 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. - * - * - *

Calculates Mu0 polarization in x direction during cyclic charge exchange. - * See M. Senba, J.Phys. B23, 1545 (1990), equations (9), (11) - - * \param muonString if eq. "Mu+" begin with Mu+ precession - */ -Double_t PSimulateMuTransition::Event(const TString muonString) -{ - TComplex complexPolX = 1.0; - Double_t muoniumPolX = 1.0; //initial polarization in x direction - Double_t eventTime, eventDiffTime, captureTime, ionizationTime; - - eventTime = 0.; - eventDiffTime = 0.; - - if (fDebugFlag) cout << "Decay time = " << fMuonDecayTime << endl; - - // charge-exchange loop until muon decays - while (1) { - if (muonString == "Mu+"){// Mu+ initial state; get next electron capture time - captureTime = NextEventTime(fCaptureRate); - eventTime += captureTime; - - if (fDebugFlag) cout << "Capture time = " << captureTime << " PolX = " << complexPolX.Re() << endl; - - if (eventTime < fMuonDecayTime) - complexPolX *= GTFunction(captureTime, "Mu+"); - else{ //muon decays; handle precession prior to muon decay - eventDiffTime = fMuonDecayTime - (eventTime - captureTime); - complexPolX *= GTFunction(eventDiffTime, "Mu+"); - break; - } - // now, we have Mu0; get next ionization time - ionizationTime = NextEventTime(fIonizationRate); - eventTime += ionizationTime; - - if (fDebugFlag) cout << "Ioniza. time = " << ionizationTime << " PolX = " << complexPolX.Re() << endl; - - if (eventTime < fMuonDecayTime) - complexPolX *= GTFunction(ionizationTime, "Mu0"); - else{ //muon decays; handle precession prior to muon decay - eventDiffTime = fMuonDecayTime - (eventTime - ionizationTime); - complexPolX *= GTFunction(eventDiffTime, "Mu0"); - break; - } - } - else{// Mu0 as initial state; get next ionization time - ionizationTime = NextEventTime(fIonizationRate); - eventTime += ionizationTime; - - if (fDebugFlag) - cout << "Mu Ioniza. time = " << ionizationTime << " PolX = " << complexPolX.Re() << endl; - - if (eventTime < fMuonDecayTime) - complexPolX *= GTFunction(ionizationTime, "Mu0"); - else{ //muon decays; handle precession prior to muon decay - eventDiffTime = fMuonDecayTime - (eventTime - ionizationTime); - complexPolX *= GTFunction(eventDiffTime, "Mu0"); - break; - } - - // Mu+ state; get next electron capture time - captureTime = NextEventTime(fCaptureRate); - eventTime += captureTime; - - if (fDebugFlag) cout << "Capture time = " << captureTime << " PolX = " << complexPolX.Re() << endl; - - if (eventTime < fMuonDecayTime) - complexPolX *= GTFunction(captureTime, "Mu+"); - else{ //muon decays; handle precession prior to muon decay - eventDiffTime = fMuonDecayTime - (eventTime - captureTime); - complexPolX *= GTFunction(eventDiffTime, "Mu+"); - break; - } - } - } - - muoniumPolX = complexPolX.Re(); - if (fDebugFlag) cout << " Final PolX = " << muoniumPolX << endl; - - return muoniumPolX; -} diff --git a/src/tests/MuTransition/PSimulateMuTransition.h b/src/tests/MuTransition/PSimulateMuTransition.h deleted file mode 100644 index aab0e8d3..00000000 --- a/src/tests/MuTransition/PSimulateMuTransition.h +++ /dev/null @@ -1,119 +0,0 @@ -/*************************************************************************** - - 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 -#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 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) - 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 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); - - 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 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) - 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 - 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 - 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); - virtual TComplex GTFunction(const Double_t &time, const TString chargeState); //!< transverse field polarization function of Mu0 or Mu+ - virtual Double_t GTSpinFlip(const Double_t &time); //!< transverse field polarization function after spin-flip collisions - virtual Double_t Event(const TString muonString); - - ClassDef(PSimulateMuTransition, 0) -}; - -#endif // _PSIMULATEMUTRANSITION_H_ diff --git a/src/tests/MuTransition/PSimulateMuTransitionLinkDef.h b/src/tests/MuTransition/PSimulateMuTransitionLinkDef.h deleted file mode 100644 index d806ffcf..00000000 --- a/src/tests/MuTransition/PSimulateMuTransitionLinkDef.h +++ /dev/null @@ -1,40 +0,0 @@ -/*************************************************************************** - - PSimulateMuTransitionLinkDef.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. * - ***************************************************************************/ - -// changed for root6 -#ifdef __CINT__ - -#pragma link off all globals; -#pragma link off all classes; -#pragma link off all functions; - -#pragma link C++ class PSimulateMuTransition+; - -#endif diff --git a/src/tests/MuTransition/runMuSimulation.C b/src/tests/MuTransition/runMuSimulation.C deleted file mode 100644 index d5fa87bb..00000000 --- a/src/tests/MuTransition/runMuSimulation.C +++ /dev/null @@ -1,270 +0,0 @@ -/*************************************************************************** - - 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. * - ***************************************************************************/ -// -// either do -// root> .L runMuSimulation.C -// root> runMuSimulation() -// -// or -// root> .x runMuSimulation.C -// - -#include "TMusrRunHeader.h" -#include "PSimulateMuTransition.h" - -#define NDECAYHISTS 2 - - -void runMuSimulation() -{ -// load libraries during root startup, defined in rootlogon.C - -// gSystem->Load("$ROOTSYS/lib/libPSimulateMuTransition"); -// gSystem->Load("$ROOTSYS/lib/libTMusrRunHeader.so"); - - char titleStr[256]; - TFolder *histosFolder; - TFolder *decayAnaModule, *scAnaModule; - TFolder *gRunHeader; - TString runTitle; - TString histogramFileName; - TObjArray Slist(0); - TMusrRunPhysicalQuantity prop; - - //prepare to run simulation; here: isotropic Mu with A0 = 100.0 MHz - UInt_t runNo = 7701; - 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.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.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 - - TTimeStamp *timeStampStart = new TTimeStamp(); - cout << endl << "Simulation started on:" << endl; - timeStampStart->Print("l"); - cout << endl; - - 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, Freq23, MuFrac23, ionRate, capRate, spinFlipRate, B); - runTitle = TString("0"); - runTitle += runNo; - runTitle += TString(titleStr); - - cout << runTitle << endl << endl; - - 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->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->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 - - // feed run info header - gRunHeader = gROOT->GetRootFolder()->AddFolder("RunHeader", "MuTransition Simulation Header Info"); - gROOT->GetListOfBrowsables()->Add(gRunHeader, "RunHeader"); -// header = new TLemRunHeader(); - TMusrRunHeader *header = new TMusrRunHeader(true); - - 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", "2016-03-01 06:20:00"); - header->Set("RunInfo/Run Stop Time", "2016-03-01 06:20:11"); - - prop.Set("Run Duration", 11.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"); - - prop.Set("Sample Temperature", MRH_UNDEFINED, T, 0.01, "K"); - header->Set("RunInfo/Sample Temperature", prop); - - 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); - - prop.Set("Implantation Energy", 0, "keV"); - header->Set("RunInfo/Implantation Energy", prop); - - prop.Set("Muon Spin Angle", 0, "degree along x"); - header->Set("RunInfo/Muon Spin Angle", prop); - - header->Set("DetectorInfo/Detector001/Name", "e+ forward"); - header->Set("DetectorInfo/Detector001/Histo Number", 1); - header->Set("DetectorInfo/Detector001/Histo Length", 18001); - header->Set("DetectorInfo/Detector001/Time Zero Bin", 0.001); //doesn't like 0.0 as time zero - header->Set("DetectorInfo/Detector001/First Good Bin", 1); - header->Set("DetectorInfo/Detector001/Last Good Bin", 18000); - - header->Set("DetectorInfo/Detector002/Name", "e+ backward"); - header->Set("DetectorInfo/Detector002/Histo Number", 2); - header->Set("DetectorInfo/Detector002/Histo Length", 18001); - header->Set("DetectorInfo/Detector002/Time Zero Bin", 0.001); - header->Set("DetectorInfo/Detector002/First Good Bin", 1); - header->Set("DetectorInfo/Detector002/Last Good Bin", 18000); - - // simulation parameters - header->Set("Simulation/Mu0 Precession frequency 12", Freq12); - 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/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); - header->Set("Simulation/Decay Asymmetry", Asym); - - header->Set("SampleEnvironmentInfo/Cryo", "no cryostat"); - header->Set("MagneticFieldEnvironmentInfo/Magnet Name", "Field along z"); - header->Set("BeamlineInfo/Name", "Monte-Carlo setup"); - - histosFolder = gROOT->GetRootFolder()->AddFolder("histos", "Histograms"); - gROOT->GetListOfBrowsables()->Add(histosFolder, "histos"); - decayAnaModule = histosFolder->AddFolder("DecayAnaModule", "muSR decay histograms"); - scAnaModule = histosFolder->AddFolder("SCAnaModule", "SlowControl 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]); - - // write file - TFile *fout = new TFile(histogramFileName.Data(), "RECREATE", "Midas MC Histograms"); - if (fout == 0) { - cout << endl << "**ERROR** Couldn't create ROOT file"; - cout << endl << endl; - exit(0); - } - - fout->cd(); - - header->FillFolder(gRunHeader); - gRunHeader->Add(&Slist); - Slist.SetName("RunSummary"); - histosFolder->Write(); - gRunHeader->Write(); - fout->Close(); - cout << "Histograms written to " << histogramFileName.Data() << endl; - - cout << endl << "Simulation stopped on:" << endl; - TTimeStamp *timeStampEnd = new TTimeStamp(); - timeStampEnd->Print("l"); - cout << endl; - -// delete fout; -// delete header; -// delete histo[0]; -// delete histo[1]; -// delete gRunHeader; -} diff --git a/src/tests/MuTransition/testAnalysis.C b/src/tests/MuTransition/testAnalysis.C deleted file mode 100644 index 762d64f7..00000000 --- a/src/tests/MuTransition/testAnalysis.C +++ /dev/null @@ -1,59 +0,0 @@ -/*************************************************************************** - - testAnalysis.C - - Author: Thomas Prokscha - Date: 26-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. * - ***************************************************************************/ -{ - //gROOT->Reset(); - Int_t rebin = 5; - Double_t t0, tbin; - - t0 = 0.; - tbin = 1.; - gROOT->ProcessLine(".x defineFit.C"); //define fit function muonLife - - histos = dynamic_cast (gDirectory->Get("histos")); - hDecay00 = dynamic_cast (histos->FindObjectAny("hDecay00")); - hDecay01 = dynamic_cast (histos->FindObjectAny("hDecay01")); - - hForward = hDecay00->Rebin(rebin, "hForward"); - hBackward = hDecay01->Rebin(rebin, "hBackward"); - - //hForward->Draw(); - hSum = (TH1*) hForward->Clone("hSum"); - hSum->Reset(); - hSum->Add(hForward, hBackward, 1., 1.); - //hSum->Draw(); - - hAsym = (TH1*) hForward->GetAsymmetry(hBackward, 1., 0.); - hAsym->SetName("Asymmetry"); - hAsym->SetTitle("Asymmetry"); - hAsym->GetXaxis()->SetTitle("Time (ns)"); - hAsym->Draw(); - - asymTFExpo->SetLineColor(2); //fit function red -}