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

Use Musr-Root run header

Browse Source
This commit is contained in:
Thomas Prokscha 2016-02-17 20:42:13 +01:00
parent 96ad763940
commit 1c5c94a86b
2 changed files with 187 additions and 217 deletions
Show Stats 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

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

@ -28,84 +28,62 @@
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/ ***************************************************************************/
#include "/apps/cern/root-git/include/TMusrRunHeader.h"
#define NDECAYHISTS 2
void runMuSimulation() void runMuSimulation()
{ {
// load library // load library
gSystem->Load("$ROOTSYS/lib/libPSimulateMuTransition"); 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 *histosFolder;
TFolder *decayAnaModule; TFolder *decayAnaModule;
TFolder *runInfo; TFolder *gRunHeader;
TString runTitle;
histosFolder = gROOT->GetRootFolder()->AddFolder("histos", "Histograms"); TString histogramFileName;
gROOT->GetListOfBrowsables()->Add(histosFolder, "histos"); TObjArray Slist(0);
decayAnaModule = histosFolder->AddFolder("DecayAnaModule", "muSR decay histograms"); TMusrRunPhysicalQuantity prop;
//prepare to run simulation; here: isotropic Mu in Germanium //prepare to run simulation; here: isotropic Mu in Germanium
UInt_t runNo = 9903; UInt_t runNo = 9900;
Double_t T = 300.; //temperature Double_t T = 300.; //temperature
Double_t capRate = 1.0;//*sqrt(T/200.); Double_t spinFlipRate = 0.001; //if spinFlipRate > 0.001 only spin-flip processes will be simulated
Double_t spinFlipRate = 0.001; Double_t capRate = 1.0;//*sqrt(T/200.); //assume that capture rate varies as sqrt(T), capRate = sigma*v*p , v ~ sqrt(T)
//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 ionRate; //assume Arrhenius behaviour ionRate = preFac*exp(-EA/kT)
Double_t EA = 100.; //activation energy (meV) 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 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 B = 100.; //field in G
Double_t Freq12 = 4463; //Mu freq of the 12 transition Double_t Bvar = 0.; //field variance
Double_t Freq34 = 4463; //Mu freq of the 34 transition Double_t Freq12 = 134.858; //Mu freq of the 12 transition
Double_t Freq23 = 4463; //Mu freq of the 23 transition Double_t Freq34 = 4328.142; //Mu freq of the 34 transition
Double_t Freq14 = 4463; //Mu freq of the 14 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 MuFrac = 1.0; //total Mu fraction
Double_t MuFrac12 = 0.5; //Mu in states 12 and 34 Double_t MuFrac12 = 2*0.266; //Mu in states 12 and 34
Double_t MuFrac23 = 0.5; //Mu in states 23 and 14 Double_t MuFrac23 = 2*0.234; //Mu in states 23 and 14
Int_t Nmuons = 1e7; //number of muons Int_t Nmuons = 1e6; //number of muons
Double_t Asym = 0.27; //muon decay asymmetry Double_t Asym = 0.27; //muon decay asymmetry
// feed run info header histogramFileName = TString("0");
TString tstr; histogramFileName += runNo;
runInfo = gROOT->GetRootFolder()->AddFolder("RunInfo", "LEM RunInfo"); histogramFileName += TString(".root");
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()); 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);
header->SetLemSetup("trivial"); runTitle = TString("0");
header->SetRunNumber(runNo); runTitle += runNo;
header->SetStartTime(0); runTitle += TString(titleStr);
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(); PSimulateMuTransition *simulateMuTransition = new PSimulateMuTransition();
if (!simulateMuTransition->IsValid()){ if (!simulateMuTransition->IsValid()){
cerr << endl << "**ERROR** while invoking PSimulateTransition" << endl; cerr << endl << "**ERROR** while invoking PSimulateTransition" << endl;
return; return;
} }
simulateMuTransition->SetMuPrecFreq12(Freq12); // MHz simulateMuTransition->SetMuPrecFreq12(Freq12); // MHz
simulateMuTransition->SetMuPrecFreq34(Freq34); // MHz simulateMuTransition->SetMuPrecFreq34(Freq34); // MHz
simulateMuTransition->SetMuPrecFreq23(Freq23); // MHz simulateMuTransition->SetMuPrecFreq23(Freq23); // MHz
@ -121,167 +99,91 @@ void runMuSimulation()
simulateMuTransition->SetDecayAsymmetry(Asym); simulateMuTransition->SetDecayAsymmetry(Asym);
simulateMuTransition->SetDebugFlag(kFALSE); // to print time and phase during charge-changing cycle 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++)
decayAnaModule->Add(histo[i]);
// write file header->Set("RunInfo/File Name", histogramFileName.Data());
tstr = TString("0"); header->Set("RunInfo/Run Title", runTitle.Data());
tstr += runNo; header->Set("RunInfo/Run Number", (Int_t) runNo);
tstr += TString(".root"); header->Set("RunInfo/Run Start Time", "0");
TFile *fout = new TFile(tstr.Data(), "RECREATE", "Midas Fake Histograms"); header->Set("RunInfo/Run Stop Time", "1");
if (fout == 0) { prop.Set("Run Duration", 1.0, "sec");
cout << endl << "**ERROR** Couldn't create ROOT file"; header->Set("RunInfo/Run Duration", prop);
cout << endl << endl; header->Set("RunInfo/Laboratory", "PSI");
exit(0); 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);
fout->cd(); header->Set("DetectorInfo/Detector001/Name", "hDecay001");
runInfo->Write(); header->Set("DetectorInfo/Detector001/Histo Number", 1);
histosFolder->Write(); header->Set("DetectorInfo/Detector001/Histo Length", 12001);
fout->Close(); header->Set("DetectorInfo/Detector001/Time Zero Bin", 0);
cout << "Histograms written to " << tstr.Data() << endl; header->Set("DetectorInfo/Detector001/First Good Bin", 1);
delete fout; header->Set("DetectorInfo/Detector001/Last Good Bin", 12001);
delete [] histo; 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);
runMuSimulation.C // simulation parameters
header->Set("Simulation/Mu Precession frequency 12", Freq12);
Author: Thomas Prokscha header->Set("Simulation/Mu Precession frequency 34", Freq34);
Date: 25-Feb-2010 header->Set("Simulation/Mu Precession frequency 23", Freq23);
header->Set("Simulation/Mu Precession frequency 14", Freq14);
$Id$ 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);
* Copyright (C) 2010 by Thomas Prokscha, Paul Scherrer Institut * header->Set("Simulation/Mu0 Spin Flip Rate", spinFlipRate);
* * header->Set("Simulation/Number of Muons", Nmuons);
* This program is free software; you can redistribute it and/or modify * header->Set("Simulation/Decay Asymmetry", Asym);
* 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"); histosFolder = gROOT->GetRootFolder()->AddFolder("histos", "Histograms");
gROOT->GetListOfBrowsables()->Add(histosFolder, "histos"); gROOT->GetListOfBrowsables()->Add(histosFolder, "histos");
decayAnaModule = histosFolder->AddFolder("DecayAnaModule", "muSR decay histograms"); decayAnaModule = histosFolder->AddFolder("DecayAnaModule", "muSR decay histograms");
//prepare to run simulation; here: isotropic Mu in Germanium TH1F *histo[NDECAYHISTS];
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]; char str[128];
for (UInt_t i=0; i<2; i++) { for (UInt_t i=0; i<NDECAYHISTS; i++) {
sprintf(str, "hDecay0%d", (Int_t)i); sprintf(str, "hDecay00%d", (Int_t)i+1);
histo[i] = new TH1F(str, str, 12001, -0.5, 12000.5); 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(); for (i=0; i<NDECAYHISTS; i++)
if (!simulateMuTransition->IsValid()) {
cerr << endl << "**ERROR** while invoking PSimulateTransition" << endl;
return;
}
simulateMuTransition->SetMuPrecFreq12(Freq12); // MHz
simulateMuTransition->SetMuPrecFreq34(Freq34); // MHz
simulateMuTransition->SetMuPrecFreq23(Freq23); // MHz
simulateMuTransition->SetMuPrecFreq14(Freq14); // MHz
simulateMuTransition->SetMuFraction(MuFrac); // initial Mu fraction
simulateMuTransition->SetMuFractionState12(MuFrac12); // Mu in states 12, 34
simulateMuTransition->SetMuFractionState23(MuFrac23); // Mu in states 23, 14
simulateMuTransition->SetBfield(B/10000.); // Tesla
simulateMuTransition->SetCaptureRate(capRate); // MHz
simulateMuTransition->SetIonizationRate(ionRate); // MHz
simulateMuTransition->SetNmuons(Nmuons);
simulateMuTransition->SetDecayAsymmetry(Asym);
simulateMuTransition->SetDebugFlag(kFALSE); // to print time and phase during charge-changing cycle
simulateMuTransition->PrintSettings();
simulateMuTransition->Run(histo[0], histo[1]);
for (UInt_t i=0; i<4; i++)
decayAnaModule->Add(histo[i]); decayAnaModule->Add(histo[i]);
// run simulation
simulateMuTransition->PrintSettings();
simulateMuTransition->Run(histo[0], histo[1]);
// write file // write file
tstr = TString("0"); TFile *fout = new TFile(histogramFileName.Data(), "RECREATE", "Midas MC Histograms");
tstr += runNo;
tstr += TString(".root");
TFile *fout = new TFile(tstr.Data(), "RECREATE", "Midas Fake Histograms");
if (fout == 0) { if (fout == 0) {
cout << endl << "**ERROR** Couldn't create ROOT file"; cout << endl << "**ERROR** Couldn't create ROOT file";
cout << endl << endl; cout << endl << endl;
@ -289,10 +191,11 @@ void runMuSimulation()
} }
fout->cd(); fout->cd();
runInfo->Write(); header->FillFolder(gRunHeader);
gRunHeader->Write();
histosFolder->Write(); histosFolder->Write();
fout->Close(); fout->Close();
cout << "Histograms written to " << tstr.Data() << endl; cout << "Histograms written to " << histogramFileName.Data() << endl;
delete fout; delete fout;
delete [] histo; delete [] histo;