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

Merge branch 'master' of ssh://gitorious.psi.ch/nemu/musrfit

Browse Source
This commit is contained in:
suter_a
2016-02-23 13:24:19 +01:00
parent fe48da74a1 7ae1b9ab5e
commit a0f7d5e934
4 changed files with 834 additions and 631 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

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

@ -73,6 +73,7 @@
using namespace std;
#include <TMath.h>
#include <TComplex.h>
#include "PSimulateMuTransition.h"
@ -104,6 +105,7 @@ PSimulateMuTransition::PSimulateMuTransition(UInt_t seed)
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.;
@ -142,7 +144,9 @@ void PSimulateMuTransition::PrintSettings() const
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 << "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;
@ -184,12 +188,17 @@ void PSimulateMuTransition::Run(TH1F *histoForward, TH1F *histoBackward)
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("");
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
@ -228,7 +237,7 @@ Double_t PSimulateMuTransition::NextEventTime(const Double_t &EventRate)
* <p>Determines phase of the muon spin
*
* \param time duration of precession (us);
* \param frequency muon spin precession frequency (MHz);
* \param chargeState charge state of Mu ("Mu+" or "Mu0")
*/
Double_t PSimulateMuTransition::PrecessionPhase(const Double_t &time, const TString chargeState)
{
@ -238,9 +247,7 @@ Double_t PSimulateMuTransition::PrecessionPhase(const Double_t &time, const TStr
if (chargeState == "Mu+")
muonPhaseX = TMath::TwoPi()*fMuonPrecFreq*time;
else if (chargeState == "Mu0"){
muoniumPolX = 0.5 *
(fMuFractionState12 * (TMath::Cos(TMath::TwoPi()*fMuPrecFreq12*time) + TMath::Cos(TMath::TwoPi()*fMuPrecFreq34*time)) +
fMuFractionState23 * (TMath::Cos(TMath::TwoPi()*fMuPrecFreq23*time) + TMath::Cos(TMath::TwoPi()*fMuPrecFreq14*time)));
muoniumPolX = GTFunction(time).Re();
muonPhaseX = TMath::ACos(muoniumPolX);
}
else
@ -249,6 +256,76 @@ Double_t PSimulateMuTransition::PrecessionPhase(const Double_t &time, const TStr
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)
//--------------------------------------------------------------------------

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

@ -34,6 +34,7 @@
#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)
@ -55,6 +56,7 @@ class PSimulateMuTransition : public TObject
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; }
@ -80,6 +82,7 @@ class PSimulateMuTransition : public TObject
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
@ -93,6 +96,8 @@ class PSimulateMuTransition : public TObject
virtual Double_t NextEventTime(const Double_t &EventRate);
// virtual Double_t PrecessionPhase(const Double_t &time, const Double_t &frequency);
virtual Double_t PrecessionPhase(const Double_t &time, const TString chargeState);
virtual 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)

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

@ -28,83 +28,62 @@
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
#include "/apps/cern/root-git/include/TMusrRunHeader.h"
#define NDECAYHISTS 2
void runMuSimulation()
{
// load library
gSystem->Load("$ROOTSYS/lib/libPSimulateMuTransition");
// generate data
// 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;
TFolder *runInfo;
histosFolder = gROOT->GetRootFolder()->AddFolder("histos", "Histograms");
gROOT->GetListOfBrowsables()->Add(histosFolder, "histos");
decayAnaModule = histosFolder->AddFolder("DecayAnaModule", "muSR decay histograms");
TFolder *gRunHeader;
TString runTitle;
TString histogramFileName;
TObjArray Slist(0);
TMusrRunPhysicalQuantity prop;
//prepare to run simulation; here: isotropic Mu in Germanium
UInt_t runNo = 9903;
Double_t T = 300.; //temperature
Double_t capRate = 1.0;//*sqrt(T/200.);
//assume that capture rate varies as sqrt(T), capRate = sigma*v*p , v ~ sqrt(T)
UInt_t runNo = 9900;
Double_t T = 300.; //temperature
Double_t spinFlipRate = 0.001; //if spinFlipRate > 0.001 only spin-flip processes will be simulated
Double_t capRate = 1.0;//*sqrt(T/200.); //assume that capture rate varies as sqrt(T), capRate = sigma*v*p , v ~ sqrt(T)
Double_t ionRate; //assume Arrhenius behaviour ionRate = preFac*exp(-EA/kT)
Double_t EA = 100.; //activation energy (meV)
ionRate = 2.9e7 * exp(-EA/(0.08625*T)); // Ge: 2.9*10^7MHz "attempt" frequency; 1K = 0.08625 meV
Double_t B = 100.; //field in G
Double_t Freq12 = 4463; //Mu freq of the 12 transition
Double_t Freq34 = 4463; //Mu freq of the 34 transition
Double_t Freq23 = 4463; //Mu freq of the 23 transition
Double_t Freq14 = 4463; //Mu freq of the 14 transition
Double_t Bvar = 0.; //field variance
Double_t Freq12 = 134.858; //Mu freq of the 12 transition
Double_t Freq34 = 4328.142; //Mu freq of the 34 transition
Double_t Freq23 = 143.713; //Mu freq of the 23 transition
Double_t Freq14 = 4606.713; //Mu freq of the 14 transition
Double_t MuFrac = 1.0; //total Mu fraction
Double_t MuFrac12 = 0.5; //Mu in states 12 and 34
Double_t MuFrac23 = 0.5; //Mu in states 23 and 14
Int_t Nmuons = 1e7; //number of muons
Double_t MuFrac12 = 2*0.266; //Mu in states 12 and 34
Double_t MuFrac23 = 2*0.234; //Mu in states 23 and 14
Int_t Nmuons = 1e6; //number of muons
Double_t Asym = 0.27; //muon decay asymmetry
// feed run info header
TString tstr;
runInfo = gROOT->GetRootFolder()->AddFolder("RunInfo", "LEM RunInfo");
gROOT->GetListOfBrowsables()->Add(runInfo, "RunInfo");
header = new TLemRunHeader();
tstr = TString("0");
tstr += runNo;
tstr += TString(" - Mu-frac 1.0, Mu12 -4463MHz (0.5), Mu34 -4463MHz(0.5), T=300K/EA=100meV, Cap. 1.0MHz, 10mT");
histogramFileName = TString("0");
histogramFileName += runNo;
histogramFileName += TString(".root");
header->SetRunTitle(tstr.Data());
header->SetLemSetup("trivial");
header->SetRunNumber(runNo);
header->SetStartTime(0);
header->SetStopTime(1);
header->SetModeratorHV(32.0, 0.01);
header->SetSampleHV(0.0, 0.01);
header->SetImpEnergy(31.8);
header->SetSampleTemperature(T, 0.001);
header->SetSampleBField(B, 0.1);
header->SetTimeResolution(1.);
header->SetNChannels(12001);
header->SetNHist(2);
header->SetOffsetPPCHistograms(20);
header->SetCuts("none");
header->SetModerator("none");
Double_t tt0[2] = {0., 0.};
header->SetTimeZero(tt0);
runInfo->Add(header); //add header to RunInfo folder
TH1F *histo[4];
char str[128];
for (UInt_t i=0; i<2; i++) {
sprintf(str, "hDecay0%d", (Int_t)i);
histo[i] = new TH1F(str, str, 12001, -0.5, 12000.5);
sprintf(str, "hDecay2%d", (Int_t)i);
histo[i+2] = new TH1F(str, str, 12001, -0.5, 12000.5);
}
sprintf(titleStr,"- Mu-frac %3.2f, Mu12 %6.2f MHz(%3.2f), Mu23 %6.2f MHz(%3.2f), ionRate %8.2f MHz, capRate %6.2f MHz, SF rate %6.2f MHz, %5.0f G", MuFrac, Freq12, MuFrac12/2, Freq23, MuFrac23/2, ionRate, capRate, spinFlipRate, B);
runTitle = TString("0");
runTitle += runNo;
runTitle += TString(titleStr);
PSimulateMuTransition *simulateMuTransition = new PSimulateMuTransition();
if (!simulateMuTransition->IsValid()) {
if (!simulateMuTransition->IsValid()){
cerr << endl << "**ERROR** while invoking PSimulateTransition" << endl;
return;
}
simulateMuTransition->SetMuPrecFreq12(Freq12); // MHz
simulateMuTransition->SetMuPrecFreq34(Freq34); // MHz
simulateMuTransition->SetMuPrecFreq23(Freq23); // MHz
@ -115,22 +94,96 @@ void runMuSimulation()
simulateMuTransition->SetBfield(B/10000.); // Tesla
simulateMuTransition->SetCaptureRate(capRate); // MHz
simulateMuTransition->SetIonizationRate(ionRate); // MHz
simulateMuTransition->SetSpinFlipRate(spinFlipRate); // MHz
simulateMuTransition->SetNmuons(Nmuons);
simulateMuTransition->SetDecayAsymmetry(Asym);
simulateMuTransition->SetDebugFlag(kFALSE); // to print time and phase during charge-changing cycle
simulateMuTransition->PrintSettings();
// feed run info header
gRunHeader = gROOT->GetRootFolder()->AddFolder("RunHeader", "MuTransition Simulation Header Info");
gROOT->GetListOfBrowsables()->Add(gRunHeader, "RunHeader");
// header = new TLemRunHeader();
header = new TMusrRunHeader(true);
header->FillFolder(gRunHeader);
gRunHeader->Add(&Slist);
Slist.SetName("RunSummary");
simulateMuTransition->Run(histo[0], histo[1]);
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");
for (UInt_t i=0; i<4; i++)
header->Set("RunInfo/File Name", histogramFileName.Data());
header->Set("RunInfo/Run Title", runTitle.Data());
header->Set("RunInfo/Run Number", (Int_t) runNo);
header->Set("RunInfo/Run Start Time", "0");
header->Set("RunInfo/Run Stop Time", "1");
prop.Set("Run Duration", 1.0, "sec");
header->Set("RunInfo/Run Duration", prop);
header->Set("RunInfo/Laboratory", "PSI");
header->Set("RunInfo/Instrument", "MC-Simulation");
prop.Set("Muon Beam Momentum", 0.0, "MeV/c");
header->Set("RunInfo/Muon Beam Momentum", prop);
header->Set("RunInfo/Muon Species", "positive muon and muonium");
header->Set("RunInfo/Muon Source", "Simulation");
header->Set("RunInfo/Setup", "Monte-Carlo setup");
header->Set("RunInfo/Comment", "Testing effect of charge-exchange or Mu0 spin flip processes on uSR signal");
header->Set("RunInfo/Sample Name", "Monte-Carlo");
header->Set("RunInfo/Sample Temperature", 300);
prop.Set("Sample Magnetic Field", MRH_UNDEFINED, B, Bvar, "G");
header->Set("RunInfo/Sample Magnetic Field", prop);
header->Set("RunInfo/No of Histos", 2);
prop.Set("Time Resolution", 1.0, "ns", "Simulation");
header->Set("RunInfo/Time Resolution", prop);
header->Set("DetectorInfo/Detector001/Name", "hDecay001");
header->Set("DetectorInfo/Detector001/Histo Number", 1);
header->Set("DetectorInfo/Detector001/Histo Length", 12001);
header->Set("DetectorInfo/Detector001/Time Zero Bin", 0);
header->Set("DetectorInfo/Detector001/First Good Bin", 1);
header->Set("DetectorInfo/Detector001/Last Good Bin", 12001);
header->Set("DetectorInfo/Detector002/Name", "hDecay002");
header->Set("DetectorInfo/Detector002/Histo Number", 1);
header->Set("DetectorInfo/Detector002/Histo Length", 12001);
header->Set("DetectorInfo/Detector002/Time Zero Bin", 0);
header->Set("DetectorInfo/Detector002/First Good Bin", 1);
header->Set("DetectorInfo/Detector002/Last Good Bin", 12001);
// simulation parameters
header->Set("Simulation/Mu Precession frequency 12", Freq12);
header->Set("Simulation/Mu Precession frequency 34", Freq34);
header->Set("Simulation/Mu Precession frequency 23", Freq23);
header->Set("Simulation/Mu Precession frequency 14", Freq14);
header->Set("Simulation/Mu Fraction", MuFrac);
header->Set("Simulation/Mu Fraction 12", MuFrac12);
header->Set("Simulation/Mu Fraction 23", MuFrac23);
header->Set("Simulation/Mu+ Capture Rate", capRate);
header->Set("Simulation/Mu0 Ionization Rate", ionRate);
header->Set("Simulation/Mu0 Spin Flip Rate", spinFlipRate);
header->Set("Simulation/Number of Muons", Nmuons);
header->Set("Simulation/Decay Asymmetry", Asym);
histosFolder = gROOT->GetRootFolder()->AddFolder("histos", "Histograms");
gROOT->GetListOfBrowsables()->Add(histosFolder, "histos");
decayAnaModule = histosFolder->AddFolder("DecayAnaModule", "muSR decay histograms");
TH1F *histo[NDECAYHISTS];
char str[128];
for (UInt_t i=0; i<NDECAYHISTS; i++) {
sprintf(str, "hDecay00%d", (Int_t)i+1);
histo[i] = new TH1F(str, str, 12001, -0.5, 12000.5);
}
for (i=0; i<NDECAYHISTS; i++)
decayAnaModule->Add(histo[i]);
// run simulation
simulateMuTransition->PrintSettings();
simulateMuTransition->Run(histo[0], histo[1]);
// write file
tstr = TString("0");
tstr += runNo;
tstr += TString(".root");
TFile *fout = new TFile(tstr.Data(), "RECREATE", "Midas Fake Histograms");
TFile *fout = new TFile(histogramFileName.Data(), "RECREATE", "Midas MC Histograms");
if (fout == 0) {
cout << endl << "**ERROR** Couldn't create ROOT file";
cout << endl << endl;
@ -138,10 +191,11 @@ void runMuSimulation()
}
fout->cd();
runInfo->Write();
header->FillFolder(gRunHeader);
gRunHeader->Write();
histosFolder->Write();
fout->Close();
cout << "Histograms written to " << tstr.Data() << endl;
cout << "Histograms written to " << histogramFileName.Data() << endl;
delete fout;
delete [] histo;