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

Added MuTransition simulation (Mu+/Mu0 charge-exchange).

Browse Source
This commit is contained in:
nemu 2010-02-26 15:39:07 +00:00
parent 96c78a0f02
commit a5e89b2070
6 changed files with 631 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

100
src/tests/MuTransition/Makefile.PSimulateMuTransition Normal file
View File

@ -0,0 +1,100 @@
#---------------------------------------------------
# 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
# 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) $(ROOTSYS)/lib
cp -pv PSimulateMuTransition.h $(ROOTSYS)/include
endif

225
src/tests/MuTransition/PSimulateMuTransition.cpp Normal file
View File

@ -0,0 +1,225 @@
/***************************************************************************
PSimulateMuTransition.cpp
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 <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
fMuCoupling = 4463.; // vacuum Mu hyperfine coupling constant
fBfield = 0.01; // magnetic field (T)
fCaptureRate = 0.01; // Mu+ capture rate (MHz)
fIonizationRate = 10.; // Mu0 ionization rate (MHz)
fInitialPhase = 0.;
fMuonPhase = fInitialPhase;
fMuonDecayTime = 0.;
fAsymmetry = 0.27;
}
//--------------------------------------------------------------------------
// 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 hyperfine couling (MHz) = " << fMuCoupling;
cout << endl << "B field (T) = " << fBfield;
cout << endl << "Mu+ electron capture rate (MHz) = " << fCaptureRate;
cout << endl << "Mu ionizatioan rate (MHz) = " << fIonizationRate;
cout << endl << "Number of particles to simulate = " << fNmuons;
cout << endl << "Initial muon spin phase (degree) = " << fInitialPhase;
cout << endl << "Decay asymmetry = " << fAsymmetry;
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++){
Event();
//cout << "Final phase: " << fMuonPhase << endl;
// 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);
* \param phase initial muon phase (degree);
*/
Double_t PSimulateMuTransition::PrecessionPhase(const Double_t &time, const Double_t &frequency)
{
return TMath::TwoPi()*frequency*time;
}
//--------------------------------------------------------------------------
// Event (private)
//--------------------------------------------------------------------------
/**
* <p> Generates "one muon event": simulate muon phase under free precession at
* external field and Mu precession
* initial muon state: Mu+
*
* \param muonDecayTime muon decay time (us)
* \param muonPhase muon spin phase in [0,180] degree
*/
void PSimulateMuTransition::Event()
{
Double_t eventTime, eventDiffTime, captureTime, ionizationTime;
Double_t muonPrecessionFreq; // MHz
muonPrecessionFreq = fMuonGyroRatio * fBfield;
fMuonPhase = TMath::TwoPi() * fInitialPhase/360.; // transform to radians
fMuonDecayTime = NextEventTime(fMuonDecayRate);
// charge-exchange loop until muon decay
eventTime = 0.;
eventDiffTime = 0.;
while (1) {
// assume Mu+ as initial state; get next electron capture time
captureTime = NextEventTime(fCaptureRate);
eventTime += captureTime;
if (eventTime < fMuonDecayTime)
fMuonPhase += PrecessionPhase(captureTime, muonPrecessionFreq);
else{ //muon decays; handle precession prior to muon decay
eventDiffTime = fMuonDecayTime - (eventTime - captureTime);
fMuonPhase += PrecessionPhase(eventDiffTime, muonPrecessionFreq);
break;
}
// now, we have Mu0; get next ionization time
ionizationTime = NextEventTime(fIonizationRate);
eventTime += ionizationTime;
if (eventTime < fMuonDecayTime)
fMuonPhase += PrecessionPhase(ionizationTime, fMuCoupling);
else{ //muon decays; handle precession prior to muon decay
eventDiffTime = fMuonDecayTime - (eventTime - ionizationTime);
fMuonPhase += PrecessionPhase(eventDiffTime, fMuCoupling);
break;
}
}
//fMuonPhase = TMath::ACos(TMath::Cos(fMuonPhase))*360./TMath::TwoPi(); //transform back to [0, 180] degree interval
return;
}

86
src/tests/MuTransition/PSimulateMuTransition.h Normal file
View File

@ -0,0 +1,86 @@
/***************************************************************************
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 SetBfield(Double_t value) { fBfield = value; } //!< sets magnetic field (T)
virtual void SetMuCoupling(Double_t value) { fMuCoupling = value; } //!< sets Mu hyperfine coupling (MHz)
virtual void SetCaptureRate(Double_t value){ fCaptureRate = value; } //!< sets Mu+ electron capture rate (MHz)
virtual void SetIonizationRate(Double_t value){ fIonizationRate = value; } //!< sets Mu0 ionization rate (MHz)
virtual void SetDecayAsymmetry(Double_t value){ fAsymmetry = value; } //!< muon decay asymmetry
virtual Bool_t IsValid() { return fValid; }
virtual void SetSeed(UInt_t seed);
virtual Double_t GetBfield() { return fBfield; } //!< returns the magnetic field (T)
virtual Double_t GetMuCoupling() { return fMuCoupling; } //!< returns the Mu hyperfine coupling (MHz)
virtual Double_t GetCaptureRate() { return fCaptureRate; } //!< returns Mu+ electron capture rate (MHz)
virtual Double_t GetIonizationRate() { return fIonizationRate; } //!< returns Mu0 ionization rate (MHz)
virtual void Run(TH1F *histoForward, TH1F *histoBackward);
private:
Bool_t fValid;
TRandom2 *fRandom;
Double_t fBfield; //!< magnetic field (T)
Double_t fMuCoupling; //!< Mu hyperfine coupling constant (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
Int_t fNmuons; //!< number of muons to simulate
virtual Double_t NextEventTime(const Double_t &EventRate);
virtual Double_t PrecessionPhase(const Double_t &time, const Double_t &frequency);
virtual void Event();
ClassDef(PSimulateMuTransition, 0)
};
#endif // _PSIMULATEMUTRANSITION_H_

39
src/tests/MuTransition/PSimulateMuTransitionLinkDef.h Normal file
View File

@ -0,0 +1,39 @@
/***************************************************************************
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. *
***************************************************************************/
#ifdef __CINT__
#pragma link off all globals;
#pragma link off all classes;
#pragma link off all functions;
#pragma link C++ class PSimulateMuTransition+;
#endif

122
src/tests/MuTransition/runMuSimulation.C Normal file
View File

@ -0,0 +1,122 @@
/***************************************************************************
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");
// feed run info header
UInt_t runNo = 9010;
TString tstr;
runInfo = gROOT->GetRootFolder()->AddFolder("RunInfo", "LEM RunInfo");
gROOT->GetListOfBrowsables()->Add(runInfo, "RunInfo");
header = new TLemRunHeader();
tstr = TString("0");
tstr += runNo;
tstr += TString(" - test");
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(0.2, 0.001);
header->SetSampleBField(-1.0, 0.1);
header->SetTimeResolution(0.1953125);
header->SetNChannels(12001);
header->SetNHist(2);
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;
}
//prepare to run simulation
simulateMuTransition->SetMuCoupling(100.); // MHz
simulateMuTransition->SetBfield(0.01); // Tesla
simulateMuTransition->SetCaptureRate(1.0); // MHz
simulateMuTransition->SetIonizationRate(500.0); // MHz
simulateMuTransition->SetNmuons(1e6);
simulateMuTransition->SetDecayAsymmetry(0.27);
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;
}

59
src/tests/MuTransition/testAnalysis.C Normal file
View File

@ -0,0 +1,59 @@
/***************************************************************************
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 <TFolder*>(gDirectory->Get("histos"));
hDecay00 = dynamic_cast <TH1F*>(histos->FindObjectAny("hDecay00"));
hDecay01 = dynamic_cast <TH1F*>(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
}