Added Mu0 spin-flip process

src/tests/MuTransition/PSimulateMuTransition.cpp

@@ -104,6 +104,7 @@ PSimulateMuTransition::PSimulateMuTransition(UInt_t seed)
   fBfield           = 0.01;  // magnetic field (T)
   fCaptureRate      = 0.01;  // Mu+ capture rate (MHz)
   fIonizationRate   = 10.;   // Mu0 ionization rate (MHz)
+  fSpinFlipRate     = 0.001; // Mu0 spin flip rate (MHz)
   fInitialPhase     = 0.;
   fMuonPhase        = fInitialPhase;
   fMuonDecayTime    = 0.;
@@ -142,7 +143,9 @@ void PSimulateMuTransition::PrintSettings() const
   cout << endl << "Mu precession frequency 14 (MHz)  = " << fMuPrecFreq14;
   cout << endl << "B field (T)                           = " << fBfield;
   cout << endl << "Mu+ electron capture rate (MHz)       = " << fCaptureRate;
-  cout << endl << "Mu ionizatioan rate (MHz)             = " << fIonizationRate;
+  cout << endl << "Mu0 ionizatioan rate (MHz)            = " << fIonizationRate;
+  cout << endl << "Mu0 spin-flip rate (MHz)              = " << fSpinFlipRate;
+  cout << endl << "!!! Note: if spin-flip rate > 0.001 only spin-flip process is considered!!!";
   cout << endl << "Decay asymmetry                       = " << fAsymmetry;
   cout << endl << "Muonium fraction                      = " << fMuFraction;
   cout << endl << "Muonium fraction state12              = " << fMuFractionState12;
@@ -184,12 +187,17 @@ void PSimulateMuTransition::Run(TH1F *histoForward, TH1F *histoBackward)
   for (i = 0; i<fNmuons; i++){
     fMuonPhase         = TMath::TwoPi() * fInitialPhase/360.; // transform to radians
     fMuonDecayTime     = NextEventTime(fMuonDecayRate);
+    
+    if (fSpinFlipRate > 0.001){// consider only Mu0 spin-flip in this case
+      fMuonPhase = TMath::ACos(GTSpinFlip(fMuonDecayTime));
+    }
+    else{
       // initial muon state Mu+ or Mu0?
       if (fRandom->Rndm() <= 1.-fMuFraction) 
         Event("Mu+");
       else
         Event("");
-
+    }
     // fill 50% in "forward", and 50% in "backward" detector to get independent
     // events in "forward" and "backward" histograms. This allows "normal" uSR
     // analysis of the data
@@ -228,7 +236,7 @@ Double_t PSimulateMuTransition::NextEventTime(const Double_t &EventRate)
  * <p>Determines phase of the muon spin
  *
  * \param time duration of precession (us);
- * \param frequency muon spin precession frequency (MHz);
+ * \param chargeState charge state of Mu ("Mu+" or  "Mu0")
  */
 Double_t PSimulateMuTransition::PrecessionPhase(const Double_t &time, const TString chargeState)
 {
@@ -238,9 +246,7 @@ Double_t PSimulateMuTransition::PrecessionPhase(const Double_t &time, const TStr
   if (chargeState == "Mu+")
     muonPhaseX = TMath::TwoPi()*fMuonPrecFreq*time;
   else if (chargeState == "Mu0"){
-    muoniumPolX = 0.5 * 
+    muoniumPolX = GTFunction(time);
-    (fMuFractionState12 * (TMath::Cos(TMath::TwoPi()*fMuPrecFreq12*time) + TMath::Cos(TMath::TwoPi()*fMuPrecFreq34*time)) + 
-     fMuFractionState23 * (TMath::Cos(TMath::TwoPi()*fMuPrecFreq23*time) + TMath::Cos(TMath::TwoPi()*fMuPrecFreq14*time)));
     muonPhaseX = TMath::ACos(muoniumPolX);
   }
   else
@@ -249,6 +255,60 @@ Double_t PSimulateMuTransition::PrecessionPhase(const Double_t &time, const TStr
   return muonPhaseX;
 }
 
+//--------------------------------------------------------------------------
+// Mu0 transverse field polarization function (private)
+//--------------------------------------------------------------------------
+/**
+ * <p>Calculates Mu0 polarization in x direction by superposition of four Mu0 frequencies
+ *
+ * \param time  (us);
+ */
+Double_t PSimulateMuTransition::GTFunction(const Double_t &time)
+{
+  Double_t muoniumPolX = 0;
+
+  muoniumPolX = 0.5 * 
+   (fMuFractionState12 * (TMath::Cos(TMath::TwoPi()*fMuPrecFreq12*time) + TMath::Cos(TMath::TwoPi()*fMuPrecFreq34*time)) + 
+    fMuFractionState23 * (TMath::Cos(TMath::TwoPi()*fMuPrecFreq23*time) + TMath::Cos(TMath::TwoPi()*fMuPrecFreq14*time)));
+  
+  return muoniumPolX;
+}
+
+//--------------------------------------------------------------------------
+// Mu0 transverse field polarization function after n spin-flip collisions (private)
+//--------------------------------------------------------------------------
+/**
+ * <p>Calculates Mu0 polarization in x direction after n spin flip collisions.
+ * See M. Senba, J.Phys. B24, 3531 (1991), equation (17)
+ *
+ * \param time  (us);
+ */
+Double_t PSimulateMuTransition::GTSpinFlip(const Double_t &time)
+{
+  Double_t muoniumPolX = 1.0; //initial polarization in x direction
+  Double_t eventTime =  0;
+  Double_t eventDiffTime = 0;
+  Double_t lastEventTime = 0;
+  
+  eventTime += NextEventTime(fSpinFlipRate);
+  if (eventTime >= time){
+   muoniumPolX = GTFunction(time); 
+  }
+  else{
+   while (eventTime < time){
+     eventDiffTime = eventTime - lastEventTime;
+     muoniumPolX = muoniumPolX * GTFunction(eventDiffTime);
+     lastEventTime = eventTime;
+     eventTime += NextEventTime(fSpinFlipRate);
+   }
+   // calculate for the last collision
+   eventDiffTime = time - lastEventTime;
+   muoniumPolX = muoniumPolX * GTFunction(eventDiffTime);
+  }
+ 
+  return muoniumPolX;
+}
+
 //--------------------------------------------------------------------------
 // Event (private)
 //--------------------------------------------------------------------------

src/tests/MuTransition/PSimulateMuTransition.h

@@ -55,6 +55,7 @@ class PSimulateMuTransition : public TObject
     virtual void SetMuPrecFreq14(Double_t value) { fMuPrecFreq14 = value; } //!< sets Mu transition frequency (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 SetSpinFlipRate(Double_t value){ fSpinFlipRate = value; }     //!< sets Mu0 spin flip rate (MHz)
     virtual void SetDecayAsymmetry(Double_t value){ fAsymmetry = value; }      //!< muon decay asymmetry
     virtual void SetMuFraction(Double_t value){ fMuFraction = value; }         //!< Muonium fraction
     virtual void SetMuFractionState12(Double_t value){ fMuFractionState12 = value; }  
@@ -80,6 +81,7 @@ class PSimulateMuTransition : public TObject
     Double_t  fMuonPrecFreq;    //!< muon precession frequency (MHz)
     Double_t  fCaptureRate;     //!< Mu+ electron capture rate (MHz)
     Double_t  fIonizationRate;  //!< Mu0 ionization rate (MHz)
+    Double_t  fSpinFlipRate;    //!< Mu0 spin-flip rate (MHz)  
     Double_t  fInitialPhase;    //!< initial muon spin phase
     Double_t  fMuonDecayTime;   //!< muon decay time (us)
     Double_t  fMuonPhase;       //!< phase of muon spin
@@ -93,6 +95,8 @@ class PSimulateMuTransition : public TObject
     virtual Double_t NextEventTime(const Double_t &EventRate);
 //     virtual Double_t PrecessionPhase(const Double_t &time, const Double_t &frequency);
     virtual Double_t PrecessionPhase(const Double_t &time, const TString chargeState);
+    virtual Double_t GTFunction(const Double_t &time); //!< transverse field polarization function of Mu0
+    virtual Double_t GTSpinFlip(const Double_t &time); //!< transverse field polarization function after spin-flip collisions
     virtual void     Event(const TString muonString);
  
   ClassDef(PSimulateMuTransition, 0)

src/tests/MuTransition/runMuSimulation.C

@@ -1,3 +1,155 @@
+<<<<<<< HEAD
+/***************************************************************************
+
+  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");
+
+  //prepare to run simulation; here: isotropic Mu in Germanium
+  UInt_t   runNo    = 9903;
+  Double_t T        = 300.; //temperature
+  Double_t capRate  = 1.0;//*sqrt(T/200.); 
+  Double_t spinFlipRate = 0.001;
+  //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];
+  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;
+  }
+
+  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->SetSpinFlipRate(spinFlipRate); // 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]);
+
+  // 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;
+}
+=======
 /***************************************************************************
 
   runMuSimulation.C
@@ -146,3 +298,4 @@ void runMuSimulation()
 
   delete [] histo;
 }
+>>>>>>> 4fec25e423493c58fa21cedec161430f90ffc10d