diff --git a/run/macros/makeDecayHists.C b/run/macros/makeDecayHists.C new file mode 100644 index 0000000..4096db0 --- /dev/null +++ b/run/macros/makeDecayHists.C @@ -0,0 +1,44 @@ +// makeDecayHists.C +// +// Simple root macro to generate positron decay histograms from musrSim data + +// #include "NewSpec.h" +// +// void makeDecayHists(const char* musrSimFileName) +// { +// gROOT->ProcessLine(".L NewSpec.C"); +// TFile *fout = new TFile(musrSimFileName, "UPDATE"); +// if (fout == 0) { +// cout << endl << "**ERROR** Couldn't open file" << musrSimFileName; +// cout << endl << endl; +// return; +// } +// fout->cd(); +// +// // TFile::Open(musrSimFileName, "update"); +// NewSpec spec; +// spec.TBLRCoinIO(1,0.5,0); //all muon decays +// +// // hT->Write();hB->Write();hL->Write();hR->Write();hAsyLR->Write();hAsyTB->Write(); +// // hEdeposited->Write();hTof->Write();hDetz->Write();hBeamSpot->Write(); +// +// fout->Close(); +// return; +// } +{ +// gROOT->ProcessLine(".L NewSpec.C"); + TH1F *hT, *hB, *hL, *hR, *hAsyLR, *hAsyTB, *hEdeposited, *hTof, *hDetz, *hBeamSpot; + TFile::Open("data/musr_13047.root", "update"); + NewSpec spec; + spec.TBLRCoinIO(1,0.5,1); //0: all muon decays, 1: muons decaying outside sample plate + hT = (TH1F*)gDirectory->FindObjectAny("hT"); hB = (TH1F*)gDirectory->FindObjectAny("hB"); + hL = (TH1F*)gDirectory->FindObjectAny("hL"); hR = (TH1F*)gDirectory->FindObjectAny("hR"); + hAsyLR = (TH1F*)gDirectory->FindObjectAny("hAsyLR"); hAsyTB = (TH1F*)gDirectory->FindObjectAny("hAsyTB"); + hEdeposited = (TH1F*)gDirectory->FindObjectAny("hEdeposited"); + hTof = (TH1F*)gDirectory->FindObjectAny("hTof"); + hDetz = (TH1F*)gDirectory->FindObjectAny("hDetz"); + hBeamSpot = (TH1F*)gDirectory->FindObjectAny("hBeamSpot"); + + hT->Write();hB->Write();hL->Write();hR->Write();hAsyLR->Write();hAsyTB->Write(); + hEdeposited->Write();hTof->Write();hDetz->Write();hBeamSpot->Write(); +} \ No newline at end of file diff --git a/run/macros/writeMusrRoot.C b/run/macros/writeMusrRoot.C new file mode 100644 index 0000000..fd41acb --- /dev/null +++ b/run/macros/writeMusrRoot.C @@ -0,0 +1,192 @@ +/*************************************************************************** + + writeMusrRoot.C + + Author: Thomas Prokscha + Date: 29-Nov-2017 + + +***************************************************************************/ + +/*************************************************************************** + * Copyright (C) 2017 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 "TMusrRunHeader.h" + +void writeMusrRoot(Int_t runNo, Double_t B, Double_t ModHV, Double_t SampleHV, Double_t ImpEnergy, Double_t TimeZero) +{ +// load libraries during root startup, defined in rootlogon.C + +// gSystem->Load("$ROOTSYS/lib/libPSimulateMuTransition"); +// gSystem->Load("$ROOTSYS/lib/libTMusrRunHeader.so"); + + TH1F *hDecay001, *hDecay002, *hDecay003, *hDecay004; + char titleStr[256]; + TFolder *histosFolder; + TFolder *decayAnaModule, *scAnaModule; + TFolder *gRunHeader; + TString runTitle; + TString histogramFileName; + TObjArray Slist(0); + TMusrRunPhysicalQuantity prop; + +// Int_t runNo = 13031; +// Double_t B = 1500.; //field in G +// Double_t ModHV = 12.0; //kV +// Double_t SampleHV = 10.0; //kV +// Double_t ImpEnergy = 1.1; //keV +// Double_t TimeZero = 550.; + + Int_t Nmuons = 2000000; //number of muons started in simulation + Double_t SpinRot = -10.0; //degree + Double_t Bvar = 0.1; //field variation in G + +// histogramFileName = TString("0"); + histogramFileName += runNo; + histogramFileName += TString("_his.root"); + sprintf(titleStr,"musrSim, WEW=%5.1f G, Mod=%5.1f kV, Sample=%5.1f kV, E=%5.1f keV, SR=%5.1f deg.", B, ModHV, SampleHV, ImpEnergy, SpinRot); +// runTitle = TString("0"); + runTitle += runNo; + runTitle += TString(", "); + runTitle += TString(titleStr); + + // feed run info header + gRunHeader = gROOT->GetRootFolder()->AddFolder("RunHeader", "MuTransition Simulation Header Info"); + gROOT->GetListOfBrowsables()->Add(gRunHeader, "RunHeader"); + 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", "musrSim LEM"); + 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", "LEM setup"); + header->Set("RunInfo/Comment", "musrSim/Geant4.10.3 beam transport simulation"); + header->Set("RunInfo/Sample Name", "Sample Plate"); + + prop.Set("Sample Temperature", MRH_UNDEFINED, 300, 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", 4); + + prop.Set("Time Resolution", 1.0, "ns", "Simulation"); + header->Set("RunInfo/Time Resolution", prop); + + prop.Set("Implantation Energy", ImpEnergy, "keV"); + header->Set("RunInfo/Implantation Energy", prop); + + prop.Set("Muon Spin Angle", SpinRot, "degree along x"); + header->Set("RunInfo/Muon Spin Angle", prop); + + header->Set("DetectorInfo/Detector001/Name", "e+ Left"); + header->Set("DetectorInfo/Detector001/Histo Number", 1); + header->Set("DetectorInfo/Detector001/Histo Length", 13000); + header->Set("DetectorInfo/Detector001/Time Zero Bin", TimeZero); + header->Set("DetectorInfo/Detector001/First Good Bin", TimeZero+50); + header->Set("DetectorInfo/Detector001/Last Good Bin", 13000); + + header->Set("DetectorInfo/Detector002/Name", "e+ Top"); + header->Set("DetectorInfo/Detector002/Histo Number", 2); + header->Set("DetectorInfo/Detector002/Histo Length", 13000); + header->Set("DetectorInfo/Detector002/Time Zero Bin", TimeZero); + header->Set("DetectorInfo/Detector002/First Good Bin", TimeZero+50); + header->Set("DetectorInfo/Detector002/Last Good Bin", 13000); + + header->Set("DetectorInfo/Detector003/Name", "e+ Right"); + header->Set("DetectorInfo/Detector003/Histo Number", 3); + header->Set("DetectorInfo/Detector003/Histo Length", 13000); + header->Set("DetectorInfo/Detector003/Time Zero Bin", TimeZero); + header->Set("DetectorInfo/Detector003/First Good Bin", TimeZero+50); + header->Set("DetectorInfo/Detector003/Last Good Bin", 13000); + + header->Set("DetectorInfo/Detector004/Name", "e+ Bottom"); + header->Set("DetectorInfo/Detector004/Histo Number", 4); + header->Set("DetectorInfo/Detector004/Histo Length", 13000); + header->Set("DetectorInfo/Detector004/Time Zero Bin", TimeZero); + header->Set("DetectorInfo/Detector004/First Good Bin", TimeZero+50); + header->Set("DetectorInfo/Detector004/Last Good Bin", 13000); + + // simulation parameters + header->Set("Simulation/Number of Muons", Nmuons); + + header->Set("SampleEnvironmentInfo/Cryo", "no cryostat"); + header->Set("MagneticFieldEnvironmentInfo/Magnet Name", "Field along z"); + header->Set("BeamlineInfo/Name", "musrSim LEM 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"); + + gDirectory->ls(); + if (gDirectory->FindObjectAny("hL;1")) cout << "Histo found" << endl; + hDecay001 = (TH1F*)gDirectory->FindObjectAny("hL;1"); hDecay002 = (TH1F*)gDirectory->FindObjectAny("hT;1"); + hDecay003 = (TH1F*)gDirectory->FindObjectAny("hR;1"); hDecay004 = (TH1F*)gDirectory->FindObjectAny("hB;1"); + + hDecay001->SetNameTitle("hDecay001", "hDecay001"); + hDecay002->SetNameTitle("hDecay002", "hDecay002"); + hDecay003->SetNameTitle("hDecay003", "hDecay003"); + hDecay004->SetNameTitle("hDecay004", "hDecay004"); + + decayAnaModule->Add(hDecay001);decayAnaModule->Add(hDecay002);decayAnaModule->Add(hDecay003);decayAnaModule->Add(hDecay004); + + // write file + TFile *fout = new TFile(histogramFileName.Data(), "RECREATE", "Midas musrSim 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 header; +// delete histo[0]; +// delete histo[1]; +// delete gRunHeader; +}