musrsim/run/macros/writeMusrRoot.C

193 lines
8.4 KiB
C

/***************************************************************************
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;
}