Added Landau energy loss distribution at carbon foil of LEM apparatus

include/musrDetectorConstruction.hh

@@ -41,6 +41,7 @@ class G4VPhysicalVolume;
 class G4Material;
 class musrDetectorMessenger;
 class musrScintSD;
+class musrMuEnergyLossLandau;
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 

include/musrMuEnergyLossLandau.hh

@@ -0,0 +1,95 @@
+/***************************************************************************
+ *  musrSim - the program for the simulation of (mainly) muSR instruments. *
+ *          More info on http://lmu.web.psi.ch/simulation/index.html .     *
+ *          musrSim is based od Geant4 (http://geant4.web.cern.ch/geant4/) *
+ *                                                                         *
+ *  Copyright (C) 2009 by Paul Scherrer Institut, 5232 Villigen PSI,       *
+ *                                                       Switzerland       *
+ *                                                                         *
+ *  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., 675 Mass Ave, Cambridge, MA 02139, USA.              *
+ ***************************************************************************/
+
+//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+//  Muonium Formation according to yield.cc function (through GetYields method).
+//  Id    : musrMuFormation.hh, v 1.4
+//  Author: Taofiq PARAISO, T. Shiroka
+//  Date  : 2007-12
+//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+                              
+#ifndef   musrMuEnergyLossLandau_h
+#define   musrMuEnergyLossLandau_h 1
+
+#include "G4VDiscreteProcess.hh"
+#include "G4ParticleTable.hh"
+
+#include <TMath.h>
+#include <TRandom.h>
+
+
+/*! musrMuEnergyLossLandau class defines the energy loss given by the Landau distribution,
+   in, for example, the thin Carbon foil of the LEM beam line 
+  */
+
+class musrMuEnergyLossLandau : public G4VDiscreteProcess
+{
+ public:
+   
+   musrMuEnergyLossLandau(const G4String& name = "MuEnergyLossLandau", // process description
+		   G4ProcessType aType = fElectromagnetic);
+
+  ~musrMuEnergyLossLandau();
+  
+  //! - Main method. muon energy loss process is executed at the END of a step. */
+  G4VParticleChange* PostStepDoIt(
+			     const G4Track&,
+			     const G4Step&);
+  
+  G4double GetMeanFreePath(const G4Track& aTrack,
+			   G4double previousStepSize,
+			   G4ForceCondition* condition);
+
+
+  //! Condition for process application (step Object).
+  G4bool CheckCondition(const G4Step& aStep);
+  
+  //! Condition for process application (step Pointer).
+  G4bool CheckCondition(const G4Step* aStep);
+  
+  
+  G4String  p_name;
+  G4bool condition;
+  
+  
+  void GetFinalEnergy( const G4Step* aStep);
+  TRandom *random;
+  static double landauMPV;
+  static double landauSigma;
+
+  // model parameters
+  G4ParticleTable* particleTable; 
+  G4ParticleDefinition* particle;
+  G4double rnd;
+  G4DynamicParticle *DP;
+  
+  //! The particle change object.
+  G4VParticleChange fParticleChange; 
+ 
+  void  PrepareSecondary(const G4Track&);
+  G4Track* aSecondary;
+
+  void InitializeSecondaries(const G4Track&);
+};
+
+#endif

src/musrDetectorConstruction.cc

@@ -29,6 +29,7 @@
 
 #include "musrScintSD.hh"
 #include "musrEventAction.hh"
+#include "musrMuEnergyLossLandau.hh"
 
 #include "G4GeometryManager.hh"
 #include "G4PhysicalVolumeStore.hh"
@@ -74,6 +75,8 @@
 #include "musrErrorMessage.hh"
 #include "musrSteppingAction.hh"
 
+using namespace std;
+
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 musrDetectorConstruction::musrDetectorConstruction(G4String steeringFileName)
@@ -1538,6 +1541,14 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	;  // processes are interpreded later in  musrPhysicsList.cc
       }
       
+      else if (strcmp(tmpString1,"SetLandauMPV")==0){
+	sscanf(&line[0],"%*s %*s %lf",&musrMuEnergyLossLandau::landauMPV); 
+	cout << "landauMPV = " << musrMuEnergyLossLandau::landauMPV << endl;
+      }
+      else if (strcmp(tmpString1,"SetLandauSigma")==0){
+	sscanf(&line[0],"%*s %*s %lf",&musrMuEnergyLossLandau::landauSigma);
+	cout << "landauSigma = " << musrMuEnergyLossLandau::landauSigma << endl;
+      }
 
       else ReportGeometryProblem(line);
       

src/musrMuEnergyLossLandau.cc

@@ -0,0 +1,132 @@
+/***************************************************************************
+ *  musrSim - the program for the simulation of (mainly) muSR instruments. *
+ *          More info on http://lmu.web.psi.ch/simulation/index.html .     *
+ *          musrSim is based od Geant4 (http://geant4.web.cern.ch/geant4/) *
+ *                                                                         *
+ *  Copyright (C) 2009 by Paul Scherrer Institut, 5232 Villigen PSI,       *
+ *                                                       Switzerland       *
+ *                                                                         *
+ *  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., 675 Mass Ave, Cambridge, MA 02139, USA.              *
+ ***************************************************************************/
+
+//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+//  Muonium Formation according to yield.cc function (through GetYields method).
+//  Id    : musrMuFormation.cc, v 1.4
+//  Author: Taofiq PARAISO, T. Shiroka
+//  Date  : 2007-12
+//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+
+#include "musrMuEnergyLossLandau.hh"
+
+using namespace std;
+
+musrMuEnergyLossLandau::musrMuEnergyLossLandau(const G4String& name, G4ProcessType aType)
+               : G4VDiscreteProcess(name, aType)
+{
+  random = new TRandom();
+}
+
+musrMuEnergyLossLandau::~musrMuEnergyLossLandau(){}
+
+double musrMuEnergyLossLandau::landauMPV   = 0.1;
+double musrMuEnergyLossLandau::landauSigma = 0.1;
+
+
+G4VParticleChange* musrMuEnergyLossLandau::PostStepDoIt(const G4Track& trackData,
+                                                 const G4Step&  aStep)
+{ // Initialize ParticleChange  (by setting all its members equal to
+  //                             the corresponding members in G4Track)
+  fParticleChange.Initialize(trackData);
+
+  G4Track theNewTrack;
+  if(CheckCondition(aStep))
+    {
+      GetFinalEnergy(&aStep);
+      G4Step theStep;
+      PrepareSecondary( trackData);      
+      fParticleChange.AddSecondary(aSecondary);
+      fParticleChange.ProposeTrackStatus(fStopAndKill) ;
+    }  
+  else
+    {
+      fParticleChange.ProposeTrackStatus(trackData.GetTrackStatus()) ;
+    }
+  return &fParticleChange;
+}
+
+
+G4bool musrMuEnergyLossLandau::CheckCondition(const G4Step& aStep)
+{ // Decide when to call the MuEnergyLossLandau - i.e. for muons going through the C foil.
+  condition=false;
+  p_name = aStep.GetTrack()->GetDefinition()->GetParticleName(); // particle name  
+  std::string logVolName = aStep.GetTrack()->GetVolume()->GetLogicalVolume()->GetName();
+  if(p_name == "mu+" && ((logVolName=="log_coulombCFoil")||(logVolName=="log_CFoil")))
+    {
+      condition=true;
+    }
+  return condition;
+}
+
+
+G4double musrMuEnergyLossLandau::GetMeanFreePath(const G4Track&, 
+                                                G4double,
+                                                G4ForceCondition* condition1)
+{
+  *condition1 = Forced;
+   return DBL_MAX;
+}
+
+
+void musrMuEnergyLossLandau::GetFinalEnergy(const G4Step* aStep)
+{    // Particle generation
+     particleTable=G4ParticleTable::GetParticleTable();     
+     G4double Eloss, Efinal;
+     G4double E = aStep->GetTrack()->GetDynamicParticle()->GetKineticEnergy()/CLHEP::keV;
+     
+     // Positive muon
+     if(p_name=="mu+"){
+      do{
+        Eloss = random->Landau(landauMPV,landauSigma);       
+      } while (Eloss < 0.);
+     
+     Efinal = E - Eloss; 
+//      cout << "E, Eloss, Efinal = " << E << ", " << Eloss << ", "  << Efinal << endl;
+ 	     
+     particle = particleTable->FindParticle(p_name) ;
+
+     // Set the new dynamic particle DP
+        DP = new G4DynamicParticle(particle, 
+	         aStep->GetTrack()->GetDynamicParticle()->GetMomentumDirection(),
+	         Efinal/1000.);
+	 
+/*	 IMPORTANT : COPY THOSE DATA TO GET THE SAME PARTICLE PROPERTIES!!!
+	 SHOULD BE KEPT WHEN BUILDING A PARTICLE CHANGE  */
+	 DP->SetProperTime(  aStep->GetTrack()->GetDynamicParticle()->GetProperTime());
+	 DP->SetPolarization(aStep->GetTrack()->GetDynamicParticle()->GetPolarization().x(),
+	                     aStep->GetTrack()->GetDynamicParticle()->GetPolarization().y(),
+			     aStep->GetTrack()->GetDynamicParticle()->GetPolarization().z());
+	 DP->SetPreAssignedDecayProperTime(
+	 aStep->GetTrack()->GetDynamicParticle()->GetPreAssignedDecayProperTime());
+     }
+}
+
+
+void musrMuEnergyLossLandau::PrepareSecondary(const G4Track& track)
+{
+  if(p_name=="mu+")
+    {
+     aSecondary = new G4Track(DP,track.GetGlobalTime(),track.GetPosition());
+    }
+}

src/musrPhysicsList.cc

@@ -263,6 +263,9 @@ void musrPhysicsList::ConstructProcess()
 // For Muonium formation in the Carbon foil
 #include "musrMuFormation.hh" // includes the yield function Y = Y(E).
 
+// For muon energy loss in thin Carbon foil
+#include "musrMuEnergyLossLandau.hh"
+
 // For a simple Muonium "scattering" when Mu hits solid materials
 #include "musrMuScatter.hh"
 
@@ -381,6 +384,7 @@ void musrPhysicsList::ConstructEM()
           //  else if (stringProcessName=="G4hIonisation")        pManager->AddProcess(new G4hIonisation,nr1,nr2,nr3);
 	//  else if (stringProcessName=="G4hLowEnergyIonisation")      pManager->AddProcess(new G4hLowEnergyIonisation,nr1,nr2,nr3);
 	  else if (stringProcessName=="musrMuFormation")              pManager->AddProcess(new musrMuFormation,nr1,nr2,nr3);
+	  else if (stringProcessName=="musrMuEnergyLossLandau")       pManager->AddProcess(new musrMuEnergyLossLandau,nr1,nr2,nr3);
 	  else if (stringProcessName=="G4Scintillation")              {
         	      G4Scintillation* myScintillation = new G4Scintillation();
 		      //		      if (musrParameters::finiteRiseTimeInScintillator) myScintillation->SetFiniteRiseTime(true);