diff --git a/src/tests/MuTransition/Makefile.PSimulateMuTransition b/src/tests/MuTransition/Makefile.PSimulateMuTransition new file mode 100644 index 00000000..f6e8fc0e --- /dev/null +++ b/src/tests/MuTransition/Makefile.PSimulateMuTransition @@ -0,0 +1,100 @@ +#--------------------------------------------------- +# Makefile.PSimulateMuTransition +# +# Author: Thomas Prokscha +# Date: 25-Feb-2010 +# +# $Id:$ +# +#--------------------------------------------------- + +#--------------------------------------------------- +# get compilation and library flags from root-config + +ROOTCFLAGS = $(shell $(ROOTSYS)/bin/root-config --cflags) +ROOTLIBS = $(shell $(ROOTSYS)/bin/root-config --libs) +ROOTGLIBS = $(shell $(ROOTSYS)/bin/root-config --glibs) + +#--------------------------------------------------- +# depending on the architecture, choose the compiler, +# linker, and the flags to use +# + +OSTYPE = $(shell uname) + +ifeq ($(OSTYPE),Linux) +OS = LINUX +endif +ifeq ($(OSTYPE),Linux-gnu) +OS = LINUX +endif +ifeq ($(OSTYPE),darwin) +OS = DARWIN +endif + +# -- Linux +ifeq ($(OS),LINUX) +CXX = g++ +CXXFLAGS = -Wall -Wno-trigraphs -fPIC +INCLUDES = -I../include +LD = g++ +LDFLAGS = -g +SOFLAGS = -O -shared +endif + +# -- Darwin +ifeq ($(OS),DARWIN) +CXX = g++ +CXXFLAGS = -Wall -Wno-trigraphs -fPIC +INCLUDES = -I../include +LD = g++ +LDFLAGS = -g +SOFLAGS = -dynamic +endif + +# the output from the root-config script: +CXXFLAGS += $(ROOTCFLAGS) +LDFLAGS += + +# the ROOT libraries (G = graphic) +LIBS = $(ROOTLIBS) -lXMLParser +GLIBS = $(ROOTGLIBS) -lXMLParser + +# some definitions: headers (used to generate *Dict* stuff), sources, objects,... +OBJS = +OBJS += PSimulateMuTransition.o PSimulateMuTransitionDict.o + +SHLIB = libPSimulateMuTransition.so + +# make the shared lib: +# +all: $(SHLIB) + +$(SHLIB): $(OBJS) + @echo "---> Building shared library $(SHLIB) ..." + /bin/rm -f $(SHLIB) + $(LD) $(OBJS) $(SOFLAGS) -o $(SHLIB) $(LIBS) + @echo "done" + +# clean up: remove all object file (and core files) +# semicolon needed to tell make there is no source +# for this target! +# +clean:; @rm -f $(OBJS) *Dict* core* + @echo "---> removing $(OBJS)" + @echo "---> removing *~" + @rm -f *~ +# +$(OBJS): %.o: %.cpp + $(CXX) $(INCLUDES) $(CXXFLAGS) -c $< + +PSimulateMuTransitionDict.cpp: PSimulateMuTransition.h PSimulateMuTransitionLinkDef.h + @echo "Generating dictionary $@..." + rootcint -f $@ -c -p $^ + +install: all + @echo "Installing shared lib: libPSimulateMuTransition.so" +ifeq ($(OS),LINUX) + cp -pv $(SHLIB) $(ROOTSYS)/lib + cp -pv PSimulateMuTransition.h $(ROOTSYS)/include +endif diff --git a/src/tests/MuTransition/PSimulateMuTransition.cpp b/src/tests/MuTransition/PSimulateMuTransition.cpp new file mode 100644 index 00000000..e989d915 --- /dev/null +++ b/src/tests/MuTransition/PSimulateMuTransition.cpp @@ -0,0 +1,225 @@ +/*************************************************************************** + + PSimulateMuTransition.cpp + + Author: Thomas Prokscha + Date: 25-Feb-2010 + + $Id$ + +***************************************************************************/ + +/*************************************************************************** + * Copyright (C) 2010 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 +using namespace std; + +#include + +#include "PSimulateMuTransition.h" + +ClassImp(PSimulateMuTransition) + +//-------------------------------------------------------------------------- +// Constructor +//-------------------------------------------------------------------------- +/** + *

Constructor. + * + * \param seed for the random number generator + */ +PSimulateMuTransition::PSimulateMuTransition(UInt_t seed) +{ + fValid = true; + + fRandom = new TRandom2(seed); + if (fRandom == 0) { + fValid = false; + } + + fNmuons = 100; // number of muons to simulate + fMuCoupling = 4463.; // vacuum Mu hyperfine coupling constant + fBfield = 0.01; // magnetic field (T) + fCaptureRate = 0.01; // Mu+ capture rate (MHz) + fIonizationRate = 10.; // Mu0 ionization rate (MHz) + fInitialPhase = 0.; + fMuonPhase = fInitialPhase; + fMuonDecayTime = 0.; + fAsymmetry = 0.27; +} + +//-------------------------------------------------------------------------- +// Destructor +//-------------------------------------------------------------------------- +/** + *

Destructor. + * + */ +PSimulateMuTransition::~PSimulateMuTransition() +{ + if (fRandom) { + delete fRandom; + fRandom = 0; + } +} +//-------------------------------------------------------------------------- +// Output of current settings +//-------------------------------------------------------------------------- +/*! + *

Prints the current settings onto std output. + */ +void PSimulateMuTransition::PrintSettings() const +{ + cout << endl << "Mu hyperfine couling (MHz) = " << fMuCoupling; + cout << endl << "B field (T) = " << fBfield; + cout << endl << "Mu+ electron capture rate (MHz) = " << fCaptureRate; + cout << endl << "Mu ionizatioan rate (MHz) = " << fIonizationRate; + cout << endl << "Number of particles to simulate = " << fNmuons; + cout << endl << "Initial muon spin phase (degree) = " << fInitialPhase; + cout << endl << "Decay asymmetry = " << fAsymmetry; + cout << endl << endl; +} + +//-------------------------------------------------------------------------- +// SetSeed (public) +//-------------------------------------------------------------------------- +/** + *

Sets the seed of the random generator. + * + * \param seed for the random number generator + */ +void PSimulateMuTransition::SetSeed(UInt_t seed) +{ + if (!fValid) + return; + + fRandom->SetSeed(seed); +} + +//-------------------------------------------------------------------------- +// Run (public) +//-------------------------------------------------------------------------- +/** + * \param histoForward + */ +void PSimulateMuTransition::Run(TH1F *histoForward, TH1F *histoBackward) +{ + Int_t i; + if (histoForward == 0 || histoBackward == 0) + return; + + for (i = 0; iRndm() <= 0.5) + histoForward->Fill(fMuonDecayTime*1000., 1. + fAsymmetry*TMath::Cos(fMuonPhase)); + else + histoBackward->Fill(fMuonDecayTime*1000., 1. - fAsymmetry*TMath::Cos(fMuonPhase)); + + if ( (i%100000) == 0) cout << "number of events processed: " << i << endl; + } + cout << "number of events processed: " << i << endl; + return; +} + +//-------------------------------------------------------------------------- +// NextEventTime (private) +//-------------------------------------------------------------------------- +/** + *

Determine time of next event, assuming "Poisson" distribution in time + * + * \param EventRate event rate in MHz; returns next event time in micro-seconds + */ +Double_t PSimulateMuTransition::NextEventTime(const Double_t &EventRate) +{ + if (EventRate <= 0.) + return -1.; // signal error + + return -1./EventRate * TMath::Log(fRandom->Rndm()); +} + +//-------------------------------------------------------------------------- +// Phase (private) +//-------------------------------------------------------------------------- +/** + *

Determines phase of the muon spin + * + * \param time duration of precession (us); + * \param frequency muon spin precession frequency (MHz); + * \param phase initial muon phase (degree); + */ +Double_t PSimulateMuTransition::PrecessionPhase(const Double_t &time, const Double_t &frequency) +{ + return TMath::TwoPi()*frequency*time; +} + +//-------------------------------------------------------------------------- +// Event (private) +//-------------------------------------------------------------------------- +/** + *

Generates "one muon event": simulate muon phase under free precession at + * external field and Mu precession + * initial muon state: Mu+ + * + * \param muonDecayTime muon decay time (us) + * \param muonPhase muon spin phase in [0,180] degree + */ +void PSimulateMuTransition::Event() +{ + Double_t eventTime, eventDiffTime, captureTime, ionizationTime; + Double_t muonPrecessionFreq; // MHz + + muonPrecessionFreq = fMuonGyroRatio * fBfield; + fMuonPhase = TMath::TwoPi() * fInitialPhase/360.; // transform to radians + fMuonDecayTime = NextEventTime(fMuonDecayRate); + + // charge-exchange loop until muon decay + eventTime = 0.; + eventDiffTime = 0.; + while (1) { + // assume Mu+ as initial state; get next electron capture time + captureTime = NextEventTime(fCaptureRate); + eventTime += captureTime; + if (eventTime < fMuonDecayTime) + fMuonPhase += PrecessionPhase(captureTime, muonPrecessionFreq); + else{ //muon decays; handle precession prior to muon decay + eventDiffTime = fMuonDecayTime - (eventTime - captureTime); + fMuonPhase += PrecessionPhase(eventDiffTime, muonPrecessionFreq); + break; + } + + // now, we have Mu0; get next ionization time + ionizationTime = NextEventTime(fIonizationRate); + eventTime += ionizationTime; + if (eventTime < fMuonDecayTime) + fMuonPhase += PrecessionPhase(ionizationTime, fMuCoupling); + else{ //muon decays; handle precession prior to muon decay + eventDiffTime = fMuonDecayTime - (eventTime - ionizationTime); + fMuonPhase += PrecessionPhase(eventDiffTime, fMuCoupling); + break; + } + } + //fMuonPhase = TMath::ACos(TMath::Cos(fMuonPhase))*360./TMath::TwoPi(); //transform back to [0, 180] degree interval + return; +} diff --git a/src/tests/MuTransition/PSimulateMuTransition.h b/src/tests/MuTransition/PSimulateMuTransition.h new file mode 100644 index 00000000..172c64af --- /dev/null +++ b/src/tests/MuTransition/PSimulateMuTransition.h @@ -0,0 +1,86 @@ +/*************************************************************************** + + PSimulateMuTransition.h + + Author: Thomas Prokscha + Date: 25-Feb-2010 + + $Id$ + +***************************************************************************/ + +/*************************************************************************** + * Copyright (C) 2010 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. * + ***************************************************************************/ + +#ifndef _PSIMULATEMUTRANSITION_H_ +#define _PSIMULATEMUTRANSITION_H_ + +#include +#include +#include + +// global constants +const Double_t fMuonGyroRatio = 135.54; //!< muon gyromagnetic ratio (MHz/T) +const Double_t fMuonDecayRate = 0.4551; //!< muon decay rate (1/tau_mu, MHz) + +class PSimulateMuTransition : public TObject +{ + public: + PSimulateMuTransition(UInt_t seed = 0); + virtual ~PSimulateMuTransition(); + + virtual void PrintSettings() const; + virtual void SetNmuons(Int_t value) { fNmuons = value; } //!< number of muons + virtual void SetBfield(Double_t value) { fBfield = value; } //!< sets magnetic field (T) + virtual void SetMuCoupling(Double_t value) { fMuCoupling = value; } //!< sets Mu hyperfine coupling (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 SetDecayAsymmetry(Double_t value){ fAsymmetry = value; } //!< muon decay asymmetry + + virtual Bool_t IsValid() { return fValid; } + virtual void SetSeed(UInt_t seed); + + virtual Double_t GetBfield() { return fBfield; } //!< returns the magnetic field (T) + virtual Double_t GetMuCoupling() { return fMuCoupling; } //!< returns the Mu hyperfine coupling (MHz) + virtual Double_t GetCaptureRate() { return fCaptureRate; } //!< returns Mu+ electron capture rate (MHz) + virtual Double_t GetIonizationRate() { return fIonizationRate; } //!< returns Mu0 ionization rate (MHz) + virtual void Run(TH1F *histoForward, TH1F *histoBackward); + + private: + Bool_t fValid; + TRandom2 *fRandom; + + Double_t fBfield; //!< magnetic field (T) + Double_t fMuCoupling; //!< Mu hyperfine coupling constant (MHz) + Double_t fCaptureRate; //!< Mu+ electron capture rate (MHz) + Double_t fIonizationRate; //!< Mu0 ionization rate (MHz) + Double_t fInitialPhase; //!< initial muon spin phase + Double_t fMuonDecayTime; //!< muon decay time (us) + Double_t fMuonPhase; //!< phase of muon spin + Double_t fAsymmetry; //!< muon decay asymmetry + Int_t fNmuons; //!< number of muons to simulate + + virtual Double_t NextEventTime(const Double_t &EventRate); + virtual Double_t PrecessionPhase(const Double_t &time, const Double_t &frequency); + virtual void Event(); + + ClassDef(PSimulateMuTransition, 0) +}; + +#endif // _PSIMULATEMUTRANSITION_H_ diff --git a/src/tests/MuTransition/PSimulateMuTransitionLinkDef.h b/src/tests/MuTransition/PSimulateMuTransitionLinkDef.h new file mode 100644 index 00000000..1538a7ee --- /dev/null +++ b/src/tests/MuTransition/PSimulateMuTransitionLinkDef.h @@ -0,0 +1,39 @@ +/*************************************************************************** + + PSimulateMuTransitionLinkDef.h + + Author: Thomas Prokscha + Date: 25-Feb-2010 + + $Id$ + +***************************************************************************/ + +/*************************************************************************** + * Copyright (C) 2010 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. * + ***************************************************************************/ + +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class PSimulateMuTransition+; + +#endif diff --git a/src/tests/MuTransition/runMuSimulation.C b/src/tests/MuTransition/runMuSimulation.C new file mode 100644 index 00000000..dfca1a18 --- /dev/null +++ b/src/tests/MuTransition/runMuSimulation.C @@ -0,0 +1,122 @@ +/*************************************************************************** + + runMuSimulation.C + + Author: Thomas Prokscha + Date: 25-Feb-2010 + + $Id$ + +***************************************************************************/ + +/*************************************************************************** + * Copyright (C) 2010 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. * + ***************************************************************************/ + +void runMuSimulation() +{ + // load library + gSystem->Load("$ROOTSYS/lib/libPSimulateMuTransition"); + + // generate data + TFolder *histosFolder; + TFolder *decayAnaModule; + TFolder *runInfo; + + histosFolder = gROOT->GetRootFolder()->AddFolder("histos", "Histograms"); + gROOT->GetListOfBrowsables()->Add(histosFolder, "histos"); + decayAnaModule = histosFolder->AddFolder("DecayAnaModule", "muSR decay histograms"); + + // feed run info header + UInt_t runNo = 9010; + TString tstr; + runInfo = gROOT->GetRootFolder()->AddFolder("RunInfo", "LEM RunInfo"); + gROOT->GetListOfBrowsables()->Add(runInfo, "RunInfo"); + header = new TLemRunHeader(); + tstr = TString("0"); + tstr += runNo; + tstr += TString(" - test"); + 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(0.2, 0.001); + header->SetSampleBField(-1.0, 0.1); + header->SetTimeResolution(0.1953125); + header->SetNChannels(12001); + header->SetNHist(2); + 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(); + if (!simulateMuTransition->IsValid()) { + cerr << endl << "**ERROR** while invoking PSimulateTransition" << endl; + return; + } + + //prepare to run simulation + simulateMuTransition->SetMuCoupling(100.); // MHz + simulateMuTransition->SetBfield(0.01); // Tesla + simulateMuTransition->SetCaptureRate(1.0); // MHz + simulateMuTransition->SetIonizationRate(500.0); // MHz + simulateMuTransition->SetNmuons(1e6); + simulateMuTransition->SetDecayAsymmetry(0.27); + + simulateMuTransition->PrintSettings(); + + simulateMuTransition->Run(histo[0], histo[1]); + + for (UInt_t i=0; i<4; i++) + decayAnaModule->Add(histo[i]); + + // write file + tstr = TString("0"); + tstr += runNo; + tstr += TString(".root"); + TFile *fout = new TFile(tstr.Data(), "RECREATE", "Midas Fake Histograms"); + if (fout == 0) { + cout << endl << "**ERROR** Couldn't create ROOT file"; + cout << endl << endl; + exit(0); + } + + fout->cd(); + runInfo->Write(); + histosFolder->Write(); + fout->Close(); + cout << "Histograms written to " << tstr.Data() << endl; + delete fout; + + delete [] histo; +} diff --git a/src/tests/MuTransition/testAnalysis.C b/src/tests/MuTransition/testAnalysis.C new file mode 100644 index 00000000..762d64f7 --- /dev/null +++ b/src/tests/MuTransition/testAnalysis.C @@ -0,0 +1,59 @@ +/*************************************************************************** + + testAnalysis.C + + Author: Thomas Prokscha + Date: 26-Feb-2010 + + $Id$ + +***************************************************************************/ + +/*************************************************************************** + * Copyright (C) 2010 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. * + ***************************************************************************/ +{ + //gROOT->Reset(); + Int_t rebin = 5; + Double_t t0, tbin; + + t0 = 0.; + tbin = 1.; + gROOT->ProcessLine(".x defineFit.C"); //define fit function muonLife + + histos = dynamic_cast (gDirectory->Get("histos")); + hDecay00 = dynamic_cast (histos->FindObjectAny("hDecay00")); + hDecay01 = dynamic_cast (histos->FindObjectAny("hDecay01")); + + hForward = hDecay00->Rebin(rebin, "hForward"); + hBackward = hDecay01->Rebin(rebin, "hBackward"); + + //hForward->Draw(); + hSum = (TH1*) hForward->Clone("hSum"); + hSum->Reset(); + hSum->Add(hForward, hBackward, 1., 1.); + //hSum->Draw(); + + hAsym = (TH1*) hForward->GetAsymmetry(hBackward, 1., 0.); + hAsym->SetName("Asymmetry"); + hAsym->SetTitle("Asymmetry"); + hAsym->GetXaxis()->SetTitle("Time (ns)"); + hAsym->Draw(); + + asymTFExpo->SetLineColor(2); //fit function red +}