LMU/musrfit
5
0
Fork 0
Code Issues 5 Pull Requests Actions Packages Projects Releases 4 Activity

merge into master

Browse Source
This commit is contained in:
suter_a
2016-03-04 13:23:12 +01:00
parent 18564964fa ec5a26d96d
commit b157d01506
316 changed files with 35986 additions and 6332 deletions
Download Patch File Download Diff File Expand all files Collapse all files

11
src/tests/GbG_LF/CMakeLists.txt Normal file
View File

@@ -0,0 +1,11 @@
cmake_minimum_required(VERSION 2.8)
project(GbG_LF)
set(DependOnLibs gsl)
set(DependOnLibs gslcblas)
add_executable(GbG_LF GbG_LF.cpp)
target_link_libraries(GbG_LF gslcblas gsl)

116
src/tests/GbG_LF/GbG_LF.cpp Normal file
View File

@@ -0,0 +1,116 @@
#include <iostream>
#include <fstream>
#include <vector>
using namespace std;
#include <gsl/gsl_integration.h>
#define GAMMA_BAR_MUON 1.35538817e-2
#define TWO_PI 6.28318530717958647692528676656
typedef struct {
double wExt;
double s0;
double s1;
} fun_params;
void syntax()
{
cout << endl;
cout << "usage: GbG_LF Bext sigma0 sigmaGbG tmax fln" << endl;
cout << " Bext in (G)" << endl;
cout << " sigma0 in (1/us)" << endl;
cout << " sigmaGbG in (1/us)" << endl;
cout << " tmax in (us)" << endl;
cout << " fln : output file name" << endl;
cout << endl << endl;
}
double pz_GbG_2(double x, void *param)
{
fun_params *par = (fun_params*) param;
double s02 = par->s0 * par->s0;
double s12 = par->s1 * par->s1;
double x2 = x*x;
double aa = 1.0+x2*s12;
return 2.0*(s02*s02+3.0*s12*s12*aa*aa+6.0*s02*s12*aa)/(pow(par->wExt, 3.0)*pow(aa, 4.5))*exp(-0.5*s02*x2/aa);
}
int main(int argc, char *argv[])
{
if (argc != 6) {
syntax();
return -1;
}
double Bext = atof(argv[1]);
double sigma0 = atof(argv[2]);
double sigma1 = atof(argv[3]);
double tmax = atof(argv[4]);
// calculate needed time step size
double dt = 1.0 / (GAMMA_BAR_MUON*Bext) / 10.0;
if (10.0/dt < 500.0)
dt = 0.02;
cout << "dt = " << dt << " (us)" << endl;
double wExt = TWO_PI * GAMMA_BAR_MUON * Bext;
double s0 = sigma0;
double s1 = sigma1;
double wExt2 = wExt*wExt;
double s02 = s0*s0;
double s12 = s1*s1;
gsl_integration_workspace *w = gsl_integration_workspace_alloc(65536);
gsl_function fun;
fun.function = &pz_GbG_2;
fun_params param = { wExt, s0, s1 };
fun.params = &param;
gsl_integration_qawo_table *tab = gsl_integration_qawo_table_alloc(wExt, tmax, GSL_INTEG_SINE, 16);
double result = 0.0, err = 0.0;
double t = 0.0, t2 = 0.0, dval=0.0, aa=0.0;
vector<double> time;
vector<double> pol, pol1, pol2;
do {
time.push_back(t);
t2 = t*t;
dval = 0.0;
// P_z^LF (GbG, 1st part)
aa = 1.0+t2*s12;
dval = 1.0 - 2.0*(s02+s12)/wExt2 + 2.0*(s02+s12*aa)/(wExt2*pow(aa,2.5))*cos(wExt*t)*exp(-0.5*s02*t2/aa);
pol1.push_back(dval);
// P_z^LF (GbG, 2nd part)
gsl_integration_qawo_table_set_length(tab, t);
gsl_integration_qawo(&fun, 0.0, 1.0e-9, 1.0e-9, 1024, w, tab, &result, &err);
pol2.push_back(result);
dval += result;
pol.push_back(dval);
t += dt;
} while (t <= tmax);
gsl_integration_workspace_free(w);
gsl_integration_qawo_table_free(tab);
ofstream fout(argv[5], ofstream::out);
fout << "% Bext (G) : " << Bext << endl;
fout << "% sigma0 (1/us) : " << sigma0 << endl;
fout << "% sigma1 (1/us) : " << sigma1 << endl;
fout << "%----" << endl;
fout << "% time (us), pol1, pol2, pol" << endl;
for (unsigned int i=0; i<time.size(); i++) {
fout << time[i] << ", " << pol1[i] << ", " << pol2[i] << ", " << pol[i] << endl;
}
fout.close();
return 0;
}

78
src/tests/MuTransition/7700.msr Normal file
View File

@@ -0,0 +1,78 @@
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

67
src/tests/MuTransition/9900.msr Normal file
View File

@@ -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

822
src/tests/MuTransition/PSimulateMuTransition.cpp
View File

@@ -1,363 +1,459 @@
/***************************************************************************
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. Up to 3 Mu0 states are considered at the moment
(analogous to MuBC in Si, B||(100)), a non-precessing signal, and two precessing
states ("nu_12" and "nu_34").
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 "3 Mu states"
(as MuBC in Si, B || (100), for example): a non-precessing,
and two precessing signals (nu_12, nu_34).
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, also determine which Mu0 state has been formed
at the capture event. Calculate muon spin precession.
4) get the next electron capture time, continue until t_d is reached.
***************************************************************************/
/***************************************************************************
* 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 <iostream>
using namespace std;
#include <TMath.h>
#include "PSimulateMuTransition.h"
ClassImp(PSimulateMuTransition)
//--------------------------------------------------------------------------
// Constructor
//--------------------------------------------------------------------------
/**
* <p> 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
fBfield = 0.01; // magnetic field (T)
fCaptureRate = 0.01; // Mu+ capture rate (MHz)
fIonizationRate = 10.; // Mu0 ionization rate (MHz)
fInitialPhase = 0.;
fMuonPhase = fInitialPhase;
fMuonDecayTime = 0.;
fAsymmetry = 0.27;
fMuFraction = 0.;
fMuFractionState1 = 0.;
fMuFractionState2 = 0.;
fDebugFlag = kFALSE;
}
//--------------------------------------------------------------------------
// Destructor
//--------------------------------------------------------------------------
/**
* <p> Destructor.
*
*/
PSimulateMuTransition::~PSimulateMuTransition()
{
if (fRandom) {
delete fRandom;
fRandom = 0;
}
}
//--------------------------------------------------------------------------
// Output of current settings
//--------------------------------------------------------------------------
/*!
* <p>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 state1 = " << fMuFractionState1;
cout << endl << "Muonium fraction state2 = " << fMuFractionState2;
cout << endl << "Number of particles to simulate = " << fNmuons;
cout << endl << "Initial muon spin phase (degree) = " << fInitialPhase;
cout << endl << "Debug flag = " << fDebugFlag;
cout << endl << endl;
}
//--------------------------------------------------------------------------
// SetSeed (public)
//--------------------------------------------------------------------------
/**
* <p>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<fNmuons; i++){
fMuonPhase = TMath::TwoPi() * fInitialPhase/360.; // transform to radians
fMuonDecayTime = NextEventTime(fMuonDecayRate);
// 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)
//--------------------------------------------------------------------------
/**
* <p>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)
//--------------------------------------------------------------------------
/**
* <p>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 Double_t &frequency)
{
return TMath::TwoPi()*frequency*time;
}
//--------------------------------------------------------------------------
// Event (private)
//--------------------------------------------------------------------------
/**
* <p> Generates "muon event": simulate muon spin phase under charge-exchange
* with "3 Mu states" (as MuBC in Si, B || (100), for example): a non-precessing,
* and two precessing signals (nu_12, nu_34).
* 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, also determine which Mu0 state has been formed
* at the capture event. Calculate muon spin precession.
* 4) get the next electron capture time, continue until t_d is reached.
*
* <p> For isotropic muonium, TF:
* nu_12 and nu_34 with equal probabilities, probability for both states fMuFractionState1
* ni_23 and nu_14 with equal probabilities, probability for both states fMuFractionState2
*
* \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;
muonPrecessionFreq = 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, muonPrecessionFreq);
else{ //muon decays; handle precession prior to muon decay
eventDiffTime = fMuonDecayTime - (eventTime - captureTime);
fMuonPhase += PrecessionPhase(eventDiffTime, muonPrecessionFreq);
break;
}
// now, we have Mu0; get next ionization time
ionizationTime = NextEventTime(fIonizationRate);
eventTime += ionizationTime;
// determine Mu state
rndm = fRandom->Rndm();
frac1 = 1. - fMuFractionState1 - fMuFractionState2; // non-precessing Mu states
frac2 = 1. - fMuFractionState2;
if ( rndm < frac1 )
muoniumPrecessionFreq = 0.;
else if (rndm >= frac1 && rndm <= frac2){
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 << " Freq = " << muoniumPrecessionFreq
<< " Phase = " << fMuonPhase << endl;
if (eventTime < fMuonDecayTime)
fMuonPhase += PrecessionPhase(ionizationTime, muoniumPrecessionFreq);
else{ //muon decays; handle precession prior to muon decay
eventDiffTime = fMuonDecayTime - (eventTime - ionizationTime);
fMuonPhase += PrecessionPhase(eventDiffTime, muoniumPrecessionFreq);
break;
}
}
else{
// Mu0 as initial state; get next ionization time
ionizationTime = NextEventTime(fIonizationRate);
eventTime += ionizationTime;
// determine Mu state
rndm = fRandom->Rndm();
frac1 = 1. - fMuFractionState1 - fMuFractionState2; // non-precessing Mu states
frac2 = 1. - fMuFractionState2;
if ( rndm < frac1 )
muoniumPrecessionFreq = 0.;
else if (rndm >= frac1 && rndm <= frac2){
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 << " Freq = " << muoniumPrecessionFreq
<< " Phase = " << fMuonPhase << endl;
if (eventTime < fMuonDecayTime)
fMuonPhase += PrecessionPhase(ionizationTime, muoniumPrecessionFreq);
else{ //muon decays; handle precession prior to muon decay
eventDiffTime = fMuonDecayTime - (eventTime - ionizationTime);
fMuonPhase += PrecessionPhase(eventDiffTime, muoniumPrecessionFreq);
break;
}
// Mu+ state; get next electron capture time
captureTime = NextEventTime(fCaptureRate);
eventTime += captureTime;
if (fDebugFlag) cout << "Capture time = " << captureTime << " Phase = " << fMuonPhase << endl;
if (eventTime < fMuonDecayTime)
fMuonPhase += PrecessionPhase(captureTime, muonPrecessionFreq);
else{ //muon decays; handle precession prior to muon decay
eventDiffTime = fMuonDecayTime - (eventTime - captureTime);
fMuonPhase += PrecessionPhase(eventDiffTime, muonPrecessionFreq);
break;
}
}
}
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 <iostream>
using namespace std;
#include <TMath.h>
#include <TComplex.h>
#include "PSimulateMuTransition.h"
ClassImp(PSimulateMuTransition)
//--------------------------------------------------------------------------
// Constructor
//--------------------------------------------------------------------------
/**
* <p> 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
//--------------------------------------------------------------------------
/**
* <p> Destructor.
*
*/
PSimulateMuTransition::~PSimulateMuTransition()
{
if (fRandom) {
delete fRandom;
fRandom = 0;
}
}
//--------------------------------------------------------------------------
// Output of current settings
//--------------------------------------------------------------------------
/*!
* <p>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)
//--------------------------------------------------------------------------
/**
* <p>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<fNmuons; i++){
fMuonPhase = TMath::TwoPi() * fInitialPhase/360.; // transform to radians
fMuonDecayTime = NextEventTime(fMuonDecayRate);
if (fSpinFlipRate > 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)
//--------------------------------------------------------------------------
/**
* <p>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)
//--------------------------------------------------------------------------
/**
* <p>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();
muonPhaseX = TMath::ACos(muoniumPolX);
}
else
muonPhaseX = 0.;
return muonPhaseX;
}
//--------------------------------------------------------------------------
// Mu0 transverse field polarization function (private)
//--------------------------------------------------------------------------
/**
* <p>Calculates Mu0 polarization in x direction by superposition of four Mu0 frequencies
*
* \param time (us);
*/
TComplex PSimulateMuTransition::GTFunction(const Double_t &time)
{
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;
}
//--------------------------------------------------------------------------
// Mu0 transverse field polarization function after n spin-flip collisions (private)
//--------------------------------------------------------------------------
/**
* <p>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).Re();
}
else{
while (eventTime < time){
eventDiffTime = eventTime - lastEventTime;
complexPolX = complexPolX * GTFunction(eventDiffTime);
lastEventTime = eventTime;
eventTime += NextEventTime(fSpinFlipRate);
}
// calculate for the last collision
eventDiffTime = time - lastEventTime;
complexPolX = complexPolX * GTFunction(eventDiffTime);
muoniumPolX = complexPolX.Re();
}
return muoniumPolX;
}
//--------------------------------------------------------------------------
// Event (private)
//--------------------------------------------------------------------------
/**
* <p> 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.
*
* <p> 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;
}

205
src/tests/MuTransition/PSimulateMuTransition.h
View File

@@ -1,99 +1,106 @@
/***************************************************************************
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 <TObject.h>
#include <TH1F.h>
#include <TRandom2.h>
// 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 SetMuFractionState1(Double_t value){ fMuFractionState1 = value; }
virtual void SetMuFractionState2(Double_t value){ fMuFractionState2 = value; }
virtual Bool_t IsValid() { return fValid; }
virtual void SetSeed(UInt_t seed);
virtual Double_t GetBfield() { return fBfield; } //!< returns the magnetic field (T)
virtual Double_t 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 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 fMuFractionState1; //!< fraction of Mu in state 1
Double_t fMuFractionState2; //!< fraction of Mu in state 2
Int_t fNmuons; //!< number of muons to simulate
Bool_t fDebugFlag; //!< debug flag
virtual Double_t NextEventTime(const Double_t &EventRate);
virtual Double_t PrecessionPhase(const Double_t &time, const Double_t &frequency);
virtual void 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 <TObject.h>
#include <TH1F.h>
#include <TRandom2.h>
#include <TComplex.h>
// 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 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);
ClassDef(PSimulateMuTransition, 0)
};
#endif // _PSIMULATEMUTRANSITION_H_

369
src/tests/MuTransition/runMuSimulation.C
View File

@@ -1,148 +1,221 @@
/***************************************************************************
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->SetMuFractionState1(MuFrac12); // Mu in states 12, 34
simulateMuTransition->SetMuFractionState2(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;
}
/***************************************************************************
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. *
***************************************************************************/
#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");
}
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.0001;//*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)
ionRate = 0.0001; //2.9e7 * exp(-EA/(0.08625*T)); // Ge: 2.9*10^7MHz "attempt" frequency; 1K = 0.08625 meV
Double_t B = 106.5; //field in G
Double_t Bvar = 0.; //field variance
Double_t Freq12 = 40.433; //Mu freq of the 12 transition
Double_t Freq34 = 59.567; //Mu freq of the 34 transition
Double_t Freq23 = 256.245; //Mu freq of the 23 transition
Double_t Freq14 = 356.245; //Mu freq of the 14 transition
Double_t MuFrac = 1.0; //total Mu fraction
Double_t MuFrac12 = 2*0.487; //Mu in states 12 and 34
Double_t MuFrac23 = 2*0.013; //Mu in states 23 and 14
Int_t Nmuons = 5e6; //number of muons
Double_t Asym = 0.27; //muon decay asymmetry
histogramFileName = TString("0");
histogramFileName += runNo;
histogramFileName += TString(".root");
sprintf(titleStr,"- complexMuPol, A0 100MHz, Mu-frac %3.2f, Mu12 %6.2f MHz(%3.2f), Mu23 %6.2f MHz(%3.2f), ionRate %8.3f MHz, capRate %6.3f MHz, SF rate %6.3f MHz, %5.1f G", MuFrac, Freq12, MuFrac12/2, Freq23, MuFrac23/2, ionRate, capRate, spinFlipRate, B);
runTitle = TString("0");
runTitle += runNo;
runTitle += TString(titleStr);
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
// 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 Fraction", MuFrac);
header->Set("Simulation/Mu0 Fraction 12", MuFrac12);
header->Set("Simulation/Mu0 Fraction 23", MuFrac23);
header->Set("Simulation/muon 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; i<NDECAYHISTS; i++) {
sprintf(str, "hDecay00%d", (Int_t)i+1);
histo[i] = new TH1F(str, str, 18001, -0.5, 18000.5);
}
for (i=0; i<NDECAYHISTS; i++)
decayAnaModule->Add(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;
delete fout;
delete [] histo;
}

184
src/tests/analyticFakeData/analyticFakeData.C
View File

@@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2012 by Andreas Suter *
* Copyright (C) 2007-2015 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *
@@ -27,6 +27,10 @@
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
#include <iostream>
#include <vector>
using namespace std;
void analyticFakeData(const TString type, UInt_t runNo)
{
// load library
@@ -40,6 +44,10 @@ void analyticFakeData(const TString type, UInt_t runNo)
TFolder *runHeader;
TFolder *runInfo;
UInt_t offset = 0;
const UInt_t noOfHistos = 4;
const UInt_t noOfChannels = 426600;
const Double_t timeResolution = 0.025; // ns
TRandom3 rand;
histosFolder = gROOT->GetRootFolder()->AddFolder("histos", "Histograms");
gROOT->GetListOfBrowsables()->Add(histosFolder, "histos");
@@ -49,8 +57,12 @@ void analyticFakeData(const TString type, UInt_t runNo)
fileName.Form("0%d.root", (Int_t)runNo);
Double_t t0[16] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
cout << ">> define t0" << endl;
vector<Double_t> t0;
for (UInt_t i=0; i<noOfHistos; i++)
t0.push_back(5000.0);
cout << ">> define header" << endl;
if (!type.CompareTo("TLemRunHeader")) {
// feed run info header
runInfo = gROOT->GetRootFolder()->AddFolder("RunInfo", "LEM RunInfo");
@@ -67,8 +79,8 @@ void analyticFakeData(const TString type, UInt_t runNo)
header->SetImpEnergy(31.8);
header->SetSampleTemperature(0.2, 0.001);
header->SetSampleBField(-1.0, 0.1);
header->SetTimeResolution(0.025);
header->SetNChannels(66601);
header->SetTimeResolution(timeResolution);
header->SetNChannels(noOfChannels);
header->SetNHist(8);
header->SetCuts("none");
header->SetModerator("none");
@@ -107,17 +119,16 @@ void analyticFakeData(const TString type, UInt_t runNo)
prop.Set("Sample Magnetic Field", 40.0, "T");
header->Set("RunInfo/Sample Magnetic Field", prop);
prop.Set("Sample Temperature", 1.0, "mK");
header->Set("RunInfo/No of Histos", 8);
prop.Set("Time Resolution", 0.025, "ns");
header->Set("RunInfo/No of Histos", (Int_t)noOfHistos);
prop.Set("Time Resolution", timeResolution, "ns");
header->Set("RunInfo/Time Resolution", prop);
vector<Int_t> ivec;
ivec.push_back(0);
// ivec.push_back(20);
header->Set("RunInfo/RedGreen Offsets", ivec);
offset = 1;
// 2nd write all the required DetectorInfo entries
for (UInt_t i=0; i<16; i++) {
for (UInt_t i=0; i<4; i++) {
tstr.Form("DetectorInfo/Detector%03d/", i+offset);
label = tstr + TString("Name");
tstr1.Form("Detector%03d", (Int_t)(i+offset));
@@ -125,32 +136,14 @@ void analyticFakeData(const TString type, UInt_t runNo)
label = tstr + TString("Histo No");
header->Set(label, (Int_t)(i+offset));
label = tstr + TString("Histo Length");
header->Set(label, 426600);
header->Set(label, (Int_t)noOfChannels);
label = tstr + TString("Time Zero Bin");
header->Set(label, t0[i]);
label = tstr + TString("First Good Bin");
header->Set(label, (Int_t)t0[i]);
label = tstr + TString("Last Good Bin");
header->Set(label, 426600);
header->Set(label, (Int_t)noOfChannels);
}
/*
for (UInt_t i=0; i<8; i++) {
tstr.Form("DetectorInfo/Detector%03d/", 20+i+offset);
label = tstr + TString("Name");
tstr1.Form("Detector%03d", (Int_t)(20+i+offset));
header->Set(label, tstr1);
label = tstr + TString("Histo No");
header->Set(label, (Int_t)(20+i+offset));
label = tstr + TString("Histo Length");
header->Set(label, 66601);
label = tstr + TString("Time Zero Bin");
header->Set(label, t0[i]);
label = tstr + TString("First Good Bin");
header->Set(label, (Int_t)t0[i]);
label = tstr + TString("Last Good Bin");
header->Set(label, 66600);
}
*/
// 3rd write required SampleEnvironmentInfo entries
header->Set("SampleEnvironmentInfo/Cryo", "virtual");
@@ -162,124 +155,138 @@ void analyticFakeData(const TString type, UInt_t runNo)
}
// create histos
Double_t n0[16] = {25.0, 24.8, 25.1, 25.0, 24.8, 24.9, 25.3, 25.0, 25.0, 24.8, 25.1, 25.0, 24.8, 24.9, 25.3, 25.0};
Double_t bkg[16] = {0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05};
cout << ">> define n0" << endl;
vector<Double_t> n0;
for (UInt_t i=0; i<noOfHistos; i++)
n0.push_back(25.0 + 2.0*rand.Rndm());
cout << ">> define bkg" << endl;
vector<Double_t> bkg;
for (UInt_t i=0; i<noOfHistos; i++)
bkg.push_back(0.05 + 0.02*rand.Rndm());
const Double_t tau = 2197.019; // ns
// asymmetry related stuff
const Double_t timeResolution = 0.025; // ns
Double_t a0[16] = {0.25, 0.26, 0.25, 0.26, 0.25, 0.26, 0.25, 0.26, 0.25, 0.26, 0.25, 0.26, 0.25, 0.26, 0.25, 0.26};
// Double_t a1[8] = {0.06, 0.07, 0.06, 0.07, 0.06, 0.07, 0.06, 0.07};
Double_t phase[16] = {5.0*TMath::Pi()/180.0, 27.5*TMath::Pi()/180.0, 50.0*TMath::Pi()/180.0, 72.5*TMath::Pi()/180.0,
95.0*TMath::Pi()/180.0, 117.5*TMath::Pi()/180.0, 140.0*TMath::Pi()/180.0, 162.5*TMath::Pi()/180.0,
185.0*TMath::Pi()/180, 207.5*TMath::Pi()/180.0, 230.0*TMath::Pi()/180.0, 252.5*TMath::Pi()/180.0,
275*TMath::Pi()/180.0, 297.5*TMath::Pi()/180.0, 320.0*TMath::Pi()/180.0, 342.5*TMath::Pi()/180.0};
cout << ">> define a0" << endl;
vector<Double_t> a0;
for (UInt_t i=0; i<noOfHistos; i++)
a0.push_back(0.25 + 0.01*rand.Rndm());
cout << ">> define phase" << endl;
vector<Double_t> phase;
for (UInt_t i=0; i<noOfHistos; i++)
phase.push_back((5.0 + 2.0*rand.Rndm())*TMath::Pi()/180.0 + TMath::TwoPi()/noOfHistos * (Double_t)i);
const Double_t gamma = 0.0000135538817; // gamma/(2pi)
Double_t bb0 = 40000.0; // field in Gauss
// Double_t bb1 = 40015.0; // field in Gauss
Double_t rate0 = 0.25/1000.0; // in 1/ns
// Double_t rate1 = 0.15/1000.0; // in 1/ns
Double_t bb0 = 5000.0; // field in Gauss
Double_t rate0 = 7.0/1000.0; // in 1/ns
Double_t frac0 = 0.5;
Double_t bb1 = bb0 + 200.0; // field in Gauss
Double_t rate1 = 0.75/1000.0; // in 1/ns
Double_t frac1 = 0.2;
Double_t bb2 = bb0 + 600.0; // field in Gauss
Double_t rate2 = 0.25/1000.0; // in 1/ns
// fake function parameters header info: only for test purposes
cout << ">> write fake header for TMusrRoot" << endl;
if (type.CompareTo("TLemRunHeader")) {
TDoubleVector dvec;
header->Set("FakeFct/Def", "N0 exp(-t/tau_mu) [1 + A exp(-1/2 (t sigma)^2) cos(gamma_mu B t + phi)] + bkg");
for (UInt_t i=0; i<16; i++)
header->Set("FakeFct/Def", "N0 exp(-t/tau_mu) [1 + sum_{k=0}^2 frac_k A_0 exp(-1/2 (t sigma_k)^2) cos(gamma_mu B_k t + phi)] + bkg");
for (UInt_t i=0; i<noOfHistos; i++)
dvec.push_back(n0[i]);
header->Set("FakeFct/N0", dvec);
dvec.clear();
for (UInt_t i=0; i<16; i++)
for (UInt_t i=0; i<noOfHistos; i++)
dvec.push_back(bkg[i]);
header->Set("FakeFct/bkg", dvec);
dvec.clear();
for (UInt_t i=0; i<16; i++)
for (UInt_t i=0; i<noOfHistos; i++)
dvec.push_back(a0[i]);
header->Set("FakeFct/A", dvec);
dvec.clear();
for (UInt_t i=0; i<16; i++)
for (UInt_t i=0; i<noOfHistos; i++)
dvec.push_back(phase[i]);
header->Set("FakeFct/phase", dvec);
prop.Set("B", bb0, "G");
header->Set("FakeFct/B", prop);
prop.Set("lambda", rate0, "1/usec");
header->Set("FakeFct/lambda", prop);
prop.Set("B0", bb0, "G");
header->Set("FakeFct/B0", prop);
prop.Set("rate0", 1.0e3*rate0, "1/usec");
header->Set("FakeFct/rate0", prop);
header->Set("FakeFct/frac0", frac0);
prop.Set("B1", bb1, "G");
header->Set("FakeFct/B1", prop);
prop.Set("rate1", 1.0e3*rate1, "1/usec");
header->Set("FakeFct/rate1", prop);
header->Set("FakeFct/frac1", frac1);
prop.Set("B2", bb2, "G");
header->Set("FakeFct/B2", prop);
prop.Set("rate2", 1.0e3*rate2, "1/usec");
header->Set("FakeFct/rate2", prop);
}
TH1F *histo[16];
cout << ">> create histo objects" << endl;
vector<TH1F*> histo;
histo.resize(noOfHistos);
char str[128];
for (UInt_t i=0; i<16; i++) {
for (UInt_t i=0; i<noOfHistos; i++) {
if (!type.CompareTo("TLemRunHeader"))
sprintf(str, "hDecay%02d", (Int_t)(i+offset));
else
sprintf(str, "hDecay%03d", (Int_t)(i+offset));
histo[i] = new TH1F(str, str, 426601, -0.5, 426600.5);
/*
if (!type.CompareTo("TLemRunHeader"))
sprintf(str, "hDecay%02d", (Int_t)(20+i+offset));
else
sprintf(str, "hDecay%03d", (Int_t)(20+i+offset));
histo[i+8] = new TH1F(str, str, 426601, -0.5, 426600.5);
*/
histo[i] = new TH1F(str, str, noOfChannels+1, -0.5, noOfChannels+0.5);
}
Double_t time;
Double_t dval, dval1;
for (UInt_t i=0; i<16; i++) {
for (UInt_t j=0; j<426600; j++) {
Double_t dval;
cout << ">> fill histos" << endl;
for (UInt_t i=0; i<noOfHistos; i++) {
for (UInt_t j=0; j<noOfChannels; j++) {
if (j<t0[i]) {
histo[i]->SetBinContent(j+1, bkg[i]);
} else {
time = (Double_t)(j-t0[i])*timeResolution;
/*
dval = (Double_t)n0[i]*TMath::Exp(-time/tau)*(1.0+a0[i]*TMath::Exp(-time*rate0)*TMath::Cos(TMath::TwoPi()*gamma*bb0*time+phase[i])+
a1[i]*TMath::Exp(-time*rate1)*TMath::Cos(TMath::TwoPi()*gamma*bb1*time+phase[i]))+(Double_t)bkg[i];
*/
dval = (Double_t)n0[i]*TMath::Exp(-time/tau)*(1.0+a0[i]*TMath::Exp(-0.5*TMath::Power(time*rate0,2))*TMath::Cos(TMath::TwoPi()*gamma*bb0*time+phase[i]))+(Double_t)bkg[i];
dval = (Double_t)n0[i]*TMath::Exp(-time/tau)*(1.0+
frac0*a0[i]*TMath::Exp(-0.5*TMath::Power(time*rate0,2))*TMath::Cos(TMath::TwoPi()*gamma*bb0*time+phase[i]) +
frac1*a0[i]*TMath::Exp(-0.5*TMath::Power(time*rate1,2))*TMath::Cos(TMath::TwoPi()*gamma*bb1*time+phase[i]) +
(1.0-frac0-frac1)*a0[i]*TMath::Exp(-0.5*TMath::Power(time*rate2,2))*TMath::Cos(TMath::TwoPi()*gamma*bb2*time+phase[i]))+(Double_t)bkg[i];
histo[i]->SetBinContent(j+1, dval);
}
}
}
/*
cout << ">> add prompt peak" << endl;
// add a promp peak
Double_t ampl = 0.0;
Double_t ampl = 50.0;
Double_t promptLambda = 500.0/1000.0; // Lorentzain in 1/ns
if (ampl != 0.0) {
for (UInt_t i=0; i<8; i++) {
for (UInt_t j=1; j<426601; j++) {
for (UInt_t i=0; i<noOfHistos; i++) {
for (UInt_t j=1; j<noOfChannels+1; j++) {
dval = histo[i]->GetBinContent(j);
dval1 = dval*0.9; // simulate a PPC
time = (Double_t)(j-t0[i])*timeResolution;
dval += ampl*TMath::Exp(-fabs(time)*promptLambda);
dval1 += ampl*TMath::Exp(-fabs(time)*promptLambda);
histo[i]->SetBinContent(j, dval);
histo[i+8]->SetBinContent(j, dval1);
}
}
}
*/
// add Poisson noise
cout << ">> add Poisson noise" << endl;
PAddPoissonNoise *addNoise = new PAddPoissonNoise();
if (!addNoise->IsValid()) {
cerr << endl << "**ERROR** while invoking PAddPoissonNoise" << endl;
return;
}
for (UInt_t i=0; i<16; i++) {
for (UInt_t i=0; i<noOfHistos; i++) {
addNoise->AddNoise(histo[i]);
}
delete addNoise;
addNoise = 0;
/*
for (UInt_t i=0; i<8; i++) {
for (UInt_t j=1; j<histo[i]->GetEntries(); j++) {
histo[i+8]->SetBinContent(j, histo[i]->GetBinContent(j));
}
}
*/
for (UInt_t i=0; i<16; i++)
cout << ">> add histos to decayAnaModule" << endl;
for (UInt_t i=0; i<noOfHistos; i++)
decayAnaModule->Add(histo[i]);
// write file
cout << ">> write file" << endl;
TFile *fout = new TFile(fileName.Data(), "RECREATE", "Midas Fake Histograms");
if (fout == 0) {
cout << endl << "**ERROR** Couldn't create ROOT file";
@@ -299,5 +306,10 @@ void analyticFakeData(const TString type, UInt_t runNo)
fout->Close();
delete fout;
delete [] histo;
/*
cout << ">> cleanup" << endl;
for (UInt_t i=0; i<noOfHistos; i+1)
delete histo[i];
*/
cout << ">> done!" << endl;
}

99
src/tests/iotests/iotests.C Normal file
View File

@@ -0,0 +1,99 @@
// in order to illustrate the c++ io stream formatting
// usage: 1) start root
// 2) root 5.34.x -> .L iotests.C // root 6.x.y -> .L iotests.C+
// 3) iotests()
#include <iostream>
using namespace std;
void dump_ioflags(const ios_base::fmtflags flags)
{
cout << endl << ">> flags decoded:" << endl;
if (flags & ios_base::boolalpha)
cout << ">> read/write bool elements as alphabetic strings (true and false)." << endl;
if (flags & ios_base::showbase)
cout << ">> write integral values preceded by their corresponding numeric base prefix." << endl;
if (flags & ios_base::showpoint)
cout << ">> write floating-point values including always the decimal point." << endl;
if (flags & ios_base::showpos)
cout << ">> write non-negative numerical values preceded by a plus sign (+)." << endl;
if (flags & ios_base::skipws)
cout << ">> skip leading whitespaces on certain input operations." << endl;
if (flags & ios_base::unitbuf)
cout << ">> flush output after each inserting operation." << endl;
if (flags & ios_base::uppercase)
cout << ">> write uppercase letters replacing lowercase letters in certain insertion operations." << endl;
if (flags & ios_base::dec)
cout << ">> read/write integral values using decimal base format." << endl;
if (flags & ios_base::hex)
cout << ">> read/write integral values using hexadecimal base format." << endl;
if (flags & ios_base::oct)
cout << ">> read/write integral values using octal base format." << endl;
if (flags & ios_base::fixed)
cout << ">> write floating point values in fixed-point notation." << endl;
if (flags & ios_base::scientific)
cout << ">> write floating-point values in scientific notation." << endl;
if (flags & ios_base::internal)
cout << ">> the output is padded to the field width by inserting fill characters at a specified internal point." << endl;
if (flags & ios_base::left)
cout << ">> the output is padded to the field width appending fill characters at the end." << endl;
if (flags & ios_base::right)
cout << ">> the output is padded to the field width by inserting fill characters at the beginning." << endl;
cout << "---" << endl;
}
void iotests()
{
Double_t dval = 13.24;
// default settings printout
cout << "default flags = " << cout.flags() << ":" << endl;
dump_ioflags(cout.flags());
cout << "dval = " << dval << endl;
cout << "++++" << endl << endl;
// printout of the 'floatfield' flags
cout << "ios_base::floatfield :" << ios_base::floatfield << endl;
dump_ioflags(ios_base::floatfield);
cout << "set floatfield flags." << endl;
cout.setf(ios_base::floatfield);
cout << "dval = " << dval << endl;
cout << "++++" << endl << endl;
// unset 'floatfield' flags in order to see if the default is recovered
cout.unsetf(ios::floatfield);
cout << "cout.unsetf(ios::floatfield) -> cout.flags()=" << cout.flags() << endl;
dump_ioflags(cout.flags());
cout << "dval = " << dval << endl;
cout << "++++" << endl << endl;
// set 'fixed' flag
cout.setf(ios_base::fixed);
cout << "cout.setf(ios_base::fixed) -> cout.flags()=" << cout.flags() << endl;
dump_ioflags(cout.flags());
cout << "dval = " << dval << endl;
cout << "++++" << endl << endl;
// unset 'fixed' flag and set 'scientific' instead
cout.unsetf(ios_base::fixed);
cout << "cout.unsetf(ios_base::fixed)" << endl;
cout.setf(ios_base::scientific);
cout << "cout.setf(ios_base::scientific) -> cout.flags()=" << cout.flags() << endl;
dump_ioflags(cout.flags());
cout << "dval = " << dval << endl;
cout << "++++" << endl << endl;
// precision playground according to http://www.cplusplus.com/reference/ios/ios_base/precision/
cout << "precision playground according to http://www.cplusplus.com/reference/ios/ios_base/precision/" << endl;
Double_t f = 3.14159;
cout.unsetf ( ios::floatfield ); // floatfield not set
cout << "cout.unsetf ( ios::floatfield ) -> cout.flags()=" << cout.flags() << endl;
dump_ioflags(cout.flags());
cout.precision(5);
cout << "set cout.precision(5)" << endl;
cout << f << endl;
cout.precision(10);
cout << "set cout.precision(10)" << endl;
cout << f << endl;
cout.setf( ios::fixed, ios::floatfield ); // floatfield set to fixed
cout << "cout.setf( ios::fixed, ios::floatfield ) -> cout.flags()=" << cout.flags() << endl;
dump_ioflags(cout.flags());
cout << f << endl;
cout << "++++" << endl << endl;
}