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

Fixed problems with run header

Browse Source
This commit is contained in:
Thomas Prokscha 2016-03-03 17:28:50 +01:00
parent 3ff1bebd9d
commit cd04164943
2 changed files with 144 additions and 47 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

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

@ -43,38 +43,38 @@ void runMuSimulation()
char titleStr[256]; char titleStr[256];
TFolder *histosFolder; TFolder *histosFolder;
TFolder *decayAnaModule; TFolder *decayAnaModule, *scAnaModule;
TFolder *gRunHeader; TFolder *gRunHeader;
TString runTitle; TString runTitle;
TString histogramFileName; TString histogramFileName;
TObjArray Slist(0); TObjArray Slist(0);
TMusrRunPhysicalQuantity prop; TMusrRunPhysicalQuantity prop;
//prepare to run simulation; here: isotropic Mu in Germanium //prepare to run simulation; here: isotropic Mu with A0 = 100.0 MHz
UInt_t runNo = 9900; UInt_t runNo = 7701;
Double_t T = 300.; //temperature Double_t T = 300.; //temperature
Double_t spinFlipRate = 0.001; //if spinFlipRate > 0.001 only spin-flip processes will be simulated Double_t EA = 100; //activation energy (meV)
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 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) Double_t ionRate; //assume Arrhenius behaviour ionRate = preFac*exp(-EA/kT)
Double_t EA = 100.; //activation energy (meV) ionRate = 0.0001; //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 = 106.5; //field in G
Double_t B = 100.; //field in G
Double_t Bvar = 0.; //field variance Double_t Bvar = 0.; //field variance
Double_t Freq12 = 134.858; //Mu freq of the 12 transition Double_t Freq12 = 40.433; //Mu freq of the 12 transition
Double_t Freq34 = 4328.142; //Mu freq of the 34 transition Double_t Freq34 = 59.567; //Mu freq of the 34 transition
Double_t Freq23 = 143.713; //Mu freq of the 23 transition Double_t Freq23 = 256.245; //Mu freq of the 23 transition
Double_t Freq14 = 4606.713; //Mu freq of the 14 transition Double_t Freq14 = 356.245; //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 = 2*0.266; //Mu in states 12 and 34 Double_t MuFrac12 = 2*0.487; //Mu in states 12 and 34
Double_t MuFrac23 = 2*0.234; //Mu in states 23 and 14 Double_t MuFrac23 = 2*0.013; //Mu in states 23 and 14
Int_t Nmuons = 1e6; //number of muons Int_t Nmuons = 5e6; //number of muons
Double_t Asym = 0.27; //muon decay asymmetry Double_t Asym = 0.27; //muon decay asymmetry
histogramFileName = TString("0"); histogramFileName = TString("0");
histogramFileName += runNo; histogramFileName += runNo;
histogramFileName += TString(".root"); histogramFileName += TString(".root");
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); 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 = TString("0");
runTitle += runNo; runTitle += runNo;
runTitle += TString(titleStr); runTitle += TString(titleStr);
@ -103,76 +103,92 @@ void runMuSimulation()
gRunHeader = gROOT->GetRootFolder()->AddFolder("RunHeader", "MuTransition Simulation Header Info"); gRunHeader = gROOT->GetRootFolder()->AddFolder("RunHeader", "MuTransition Simulation Header Info");
gROOT->GetListOfBrowsables()->Add(gRunHeader, "RunHeader"); gROOT->GetListOfBrowsables()->Add(gRunHeader, "RunHeader");
// header = new TLemRunHeader(); // header = new TLemRunHeader();
header = new TMusrRunHeader(true); TMusrRunHeader *header = new TMusrRunHeader(true);
header->FillFolder(gRunHeader);
gRunHeader->Add(&Slist);
Slist.SetName("RunSummary");
header->Set("RunInfo/Generic Validator URL", "http://lmu.web.psi.ch/facilities/software/MusrRoot/validation/MusrRoot.xsd"); 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/Specific Validator URL", "http://lmu.web.psi.ch/facilities/software/MusrRoot/validation/MusrRootLEM.xsd");
header->Set("RunInfo/Generator", "runMuSimulation"); header->Set("RunInfo/Generator", "runMuSimulation");
header->Set("RunInfo/File Name", histogramFileName.Data()); header->Set("RunInfo/File Name", histogramFileName.Data());
header->Set("RunInfo/Run Title", runTitle.Data()); header->Set("RunInfo/Run Title", runTitle.Data());
header->Set("RunInfo/Run Number", (Int_t) runNo); header->Set("RunInfo/Run Number", (Int_t) runNo);
header->Set("RunInfo/Run Start Time", "0"); header->Set("RunInfo/Run Start Time", "2016-03-01 06:20:00");
header->Set("RunInfo/Run Stop Time", "1"); header->Set("RunInfo/Run Stop Time", "2016-03-01 06:20:11");
prop.Set("Run Duration", 1.0, "sec");
prop.Set("Run Duration", 11.0, "sec");
header->Set("RunInfo/Run Duration", prop); header->Set("RunInfo/Run Duration", prop);
header->Set("RunInfo/Laboratory", "PSI"); header->Set("RunInfo/Laboratory", "PSI");
header->Set("RunInfo/Instrument", "MC-Simulation"); header->Set("RunInfo/Instrument", "MC-Simulation");
prop.Set("Muon Beam Momentum", 0.0, "MeV/c"); prop.Set("Muon Beam Momentum", 0.0, "MeV/c");
header->Set("RunInfo/Muon Beam Momentum", prop); header->Set("RunInfo/Muon Beam Momentum", prop);
header->Set("RunInfo/Muon Species", "positive muon and muonium"); header->Set("RunInfo/Muon Species", "positive muon and muonium");
header->Set("RunInfo/Muon Source", "Simulation"); header->Set("RunInfo/Muon Source", "Simulation");
header->Set("RunInfo/Setup", "Monte-Carlo setup"); 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/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 Name", "Monte-Carlo");
header->Set("RunInfo/Sample Temperature", 300);
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"); prop.Set("Sample Magnetic Field", MRH_UNDEFINED, B, Bvar, "G");
header->Set("RunInfo/Sample Magnetic Field", prop); header->Set("RunInfo/Sample Magnetic Field", prop);
header->Set("RunInfo/No of Histos", 2); header->Set("RunInfo/No of Histos", 2);
prop.Set("Time Resolution", 1.0, "ns", "Simulation"); prop.Set("Time Resolution", 1.0, "ns", "Simulation");
header->Set("RunInfo/Time Resolution", prop); header->Set("RunInfo/Time Resolution", prop);
header->Set("DetectorInfo/Detector001/Name", "hDecay001"); prop.Set("Implantation Energy", 0, "keV");
header->Set("DetectorInfo/Detector001/Histo Number", 1); header->Set("RunInfo/Implantation Energy", prop);
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"); prop.Set("Muon Spin Angle", 0, "degree along x");
header->Set("DetectorInfo/Detector002/Histo Number", 1); header->Set("RunInfo/Muon Spin Angle", prop);
header->Set("DetectorInfo/Detector002/Histo Length", 12001);
header->Set("DetectorInfo/Detector002/Time Zero Bin", 0); 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/First Good Bin", 1);
header->Set("DetectorInfo/Detector002/Last Good Bin", 12001); header->Set("DetectorInfo/Detector002/Last Good Bin", 18000);
// simulation parameters // simulation parameters
header->Set("Simulation/Mu Precession frequency 12", Freq12); header->Set("Simulation/Mu0 Precession frequency 12", Freq12);
header->Set("Simulation/Mu Precession frequency 34", Freq34); header->Set("Simulation/Mu0 Precession frequency 34", Freq34);
header->Set("Simulation/Mu Precession frequency 23", Freq23); header->Set("Simulation/Mu0 Precession frequency 23", Freq23);
header->Set("Simulation/Mu Precession frequency 14", Freq14); header->Set("Simulation/Mu0 Precession frequency 14", Freq14);
header->Set("Simulation/Mu Fraction", MuFrac); header->Set("Simulation/Mu0 Fraction", MuFrac);
header->Set("Simulation/Mu Fraction 12", MuFrac12); header->Set("Simulation/Mu0 Fraction 12", MuFrac12);
header->Set("Simulation/Mu Fraction 23", MuFrac23); header->Set("Simulation/Mu0 Fraction 23", MuFrac23);
header->Set("Simulation/Mu+ Capture Rate", capRate); header->Set("Simulation/muon Capture Rate", capRate);
header->Set("Simulation/Mu0 Ionization Rate", ionRate); header->Set("Simulation/Mu0 Ionization Rate", ionRate);
header->Set("Simulation/Mu0 Spin Flip Rate", spinFlipRate); header->Set("Simulation/Mu0 Spin Flip Rate", spinFlipRate);
header->Set("Simulation/Number of Muons", Nmuons); header->Set("Simulation/Number of Muons", Nmuons);
header->Set("Simulation/Decay Asymmetry", Asym); 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"); 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");
scAnaModule = histosFolder->AddFolder("SCAnaModule", "SlowControl histograms");
TH1F *histo[NDECAYHISTS]; TH1F *histo[NDECAYHISTS];
char str[128]; char str[128];
for (UInt_t i=0; i<NDECAYHISTS; i++) { for (UInt_t i=0; i<NDECAYHISTS; i++) {
sprintf(str, "hDecay00%d", (Int_t)i+1); 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, 18001, -0.5, 18000.5);
} }
for (i=0; i<NDECAYHISTS; i++) for (i=0; i<NDECAYHISTS; i++)
@ -191,9 +207,12 @@ void runMuSimulation()
} }
fout->cd(); fout->cd();
header->FillFolder(gRunHeader); header->FillFolder(gRunHeader);
gRunHeader->Write(); gRunHeader->Add(&Slist);
Slist.SetName("RunSummary");
histosFolder->Write(); histosFolder->Write();
gRunHeader->Write();
fout->Close(); fout->Close();
cout << "Histograms written to " << histogramFileName.Data() << endl; cout << "Histograms written to " << histogramFileName.Data() << endl;
delete fout; delete fout;