electromagnetic/TestEm0/DirectAccess.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 /// \file electromagnetic/TestEm0/DirectAccess.cc
0027 /// \brief Main program of the electromagnetic/TestEm0 example
0028 //
0029 //
0030 //
0031 // ------------------------------------------------------------
0032 //
0033 //  To print cross sections per atom and mean free path for simple material
0034 //
0035 #include "G4BetheBlochModel.hh"
0036 #include "G4BetheHeitlerModel.hh"
0037 #include "G4BraggModel.hh"
0038 #include "G4DataVector.hh"
0039 #include "G4Electron.hh"
0040 #include "G4Gamma.hh"
0041 #include "G4KleinNishinaCompton.hh"
0042 #include "G4Material.hh"
0043 #include "G4MollerBhabhaModel.hh"
0044 #include "G4MuBetheBlochModel.hh"
0045 #include "G4MuBremsstrahlungModel.hh"
0046 #include "G4MuPairProductionModel.hh"
0047 #include "G4MuonPlus.hh"
0048 #include "G4NistManager.hh"
0049 #include "G4PEEffectFluoModel.hh"
0050 #include "G4ParticleTable.hh"
0051 #include "G4Positron.hh"
0052 #include "G4Proton.hh"
0053 #include "G4SeltzerBergerModel.hh"
0054 #include "G4SystemOfUnits.hh"
0055 #include "G4UnitsTable.hh"
0056 #include "G4eeToTwoGammaModel.hh"
0057 #include "globals.hh"
0058
0059 int main()
0060 {
0061   G4UnitDefinition::BuildUnitsTable();
0062
0063   G4ParticleDefinition* gamma = G4Gamma::Gamma();
0064   G4ParticleDefinition* posit = G4Positron::Positron();
0065   G4ParticleDefinition* elec = G4Electron::Electron();
0066   G4ParticleDefinition* prot = G4Proton::Proton();
0067   G4ParticleDefinition* muon = G4MuonPlus::MuonPlus();
0068   G4ParticleTable* partTable = G4ParticleTable::GetParticleTable();
0069   partTable->SetReadiness();
0070
0071   G4DataVector cuts;
0072   cuts.push_back(1 * keV);
0073
0074   // define materials
0075   //
0076   G4Material* material = G4NistManager::Instance()->FindOrBuildMaterial("G4_Fe");
0077
0078   G4cout << *(G4Material::GetMaterialTable()) << G4endl;
0079
0080   G4MaterialCutsCouple* couple = new G4MaterialCutsCouple(material);
0081   couple->SetIndex(0);
0082
0083   // work only for simple materials
0084   G4double Z = material->GetZ();
0085   G4double A = material->GetA();
0086
0087   // initialise gamma processes (models)
0088   //
0089   G4VEmModel* phot = new G4PEEffectFluoModel();
0090   G4VEmModel* comp = new G4KleinNishinaCompton();
0091   G4VEmModel* conv = new G4BetheHeitlerModel();
0092   phot->Initialise(gamma, cuts);
0093   comp->Initialise(gamma, cuts);
0094   conv->Initialise(gamma, cuts);
0095
0096   // valid pointer to a couple is needed for this model
0097   phot->SetCurrentCouple(couple);
0098
0099   // compute CrossSection per atom and MeanFreePath
0100   //
0101   G4double Emin = 1.01 * MeV, Emax = 2.01 * MeV, dE = 100 * keV;
0102
0103   G4cout << "\n #### Gamma : CrossSectionPerAtom and MeanFreePath for " << material->GetName()
0104          << G4endl;
0105   G4cout << "\n Energy \t PhotoElec \t Compton \t Conversion \t";
0106   G4cout << "\t PhotoElec \t Compton \t Conversion" << G4endl;
0107
0108   for (G4double Energy = Emin; Energy <= Emax; Energy += dE) {
0109     G4cout << "\n " << G4BestUnit(Energy, "Energy") << "\t"
0110            << G4BestUnit(phot->ComputeCrossSectionPerAtom(gamma, Energy, Z), "Surface") << "\t"
0111            << G4BestUnit(comp->ComputeCrossSectionPerAtom(gamma, Energy, Z), "Surface") << "\t"
0112            << G4BestUnit(conv->ComputeCrossSectionPerAtom(gamma, Energy, Z), "Surface") << "\t \t"
0113            << G4BestUnit(phot->ComputeMeanFreePath(gamma, Energy, material), "Length") << "\t"
0114            << G4BestUnit(comp->ComputeMeanFreePath(gamma, Energy, material), "Length") << "\t"
0115            << G4BestUnit(conv->ComputeMeanFreePath(gamma, Energy, material), "Length");
0116   }
0117
0118   G4cout << G4endl;
0119
0120   // initialise positron annihilation (model)
0121   //
0122   G4VEmModel* anni = new G4eeToTwoGammaModel();
0123   anni->Initialise(posit, cuts);
0124
0125   // compute CrossSection per atom and MeanFreePath
0126   //
0127   Emin = 1.01 * MeV;
0128   Emax = 2.01 * MeV;
0129   dE = 100 * keV;
0130
0131   G4cout << "\n #### e+ annihilation : CrossSectionPerAtom and MeanFreePath"
0132          << " for " << material->GetName() << G4endl;
0133   G4cout << "\n Energy \t e+ annihil \t";
0134   G4cout << "\t e+ annihil" << G4endl;
0135
0136   for (G4double Energy = Emin; Energy <= Emax; Energy += dE) {
0137     G4cout << "\n " << G4BestUnit(Energy, "Energy") << "\t"
0138            << G4BestUnit(anni->ComputeCrossSectionPerAtom(posit, Energy, Z), "Surface") << "\t \t"
0139            << G4BestUnit(anni->ComputeMeanFreePath(posit, Energy, material), "Length");
0140   }
0141
0142   G4cout << G4endl;
0143
0144   // initialise electron processes (models)
0145   //
0146   G4VEmModel* ioni = new G4MollerBhabhaModel();
0147   G4VEmModel* brem = new G4SeltzerBergerModel();
0148   ioni->Initialise(elec, cuts);
0149   brem->Initialise(elec, cuts);
0150
0151   // compute CrossSection per atom and MeanFreePath
0152   //
0153   Emin = 1.01 * MeV;
0154   Emax = 101.01 * MeV;
0155   dE = 10 * MeV;
0156   G4double Ecut = 100 * keV;
0157
0158   G4cout << "\n ####electron: CrossSection, MeanFreePath and StoppingPower"
0159          << " for " << material->GetName() << ";\tEnergy cut = " << G4BestUnit(Ecut, "Energy")
0160          << G4endl;
0161
0162   G4cout << "\n Energy \t ionization \t bremsstra \t";
0163   G4cout << "\t ionization \t bremsstra \t";
0164   G4cout << "\t ionization \t bremsstra" << G4endl;
0165
0166   for (G4double Energy = Emin; Energy <= Emax; Energy += dE) {
0167     G4cout << "\n " << G4BestUnit(Energy, "Energy") << "\t"
0168            << G4BestUnit(ioni->ComputeCrossSectionPerAtom(elec, Energy, Z, A, Ecut), "Surface")
0169            << "\t"
0170            << G4BestUnit(brem->ComputeCrossSectionPerAtom(elec, Energy, Z, A, Ecut), "Surface")
0171            << "\t \t"
0172            << G4BestUnit(ioni->ComputeMeanFreePath(elec, Energy, material, Ecut), "Length") << "\t"
0173            << G4BestUnit(brem->ComputeMeanFreePath(elec, Energy, material, Ecut), "Length")
0174            << "\t \t"
0175            << G4BestUnit(ioni->ComputeDEDXPerVolume(material, elec, Energy, Ecut), "Energy/Length")
0176            << "\t"
0177            << G4BestUnit(brem->ComputeDEDXPerVolume(material, elec, Energy, Ecut), "Energy/Length");
0178   }
0179
0180   G4cout << G4endl;
0181
0182   // initialise proton processes (models)
0183   //
0184   ioni = new G4BetheBlochModel();
0185   ioni->Initialise(prot, cuts);
0186
0187   // compute CrossSection per atom and MeanFreePath
0188   //
0189   Emin = 1.01 * MeV;
0190   Emax = 102.01 * MeV;
0191   dE = 10 * MeV;
0192   Ecut = 100 * keV;
0193
0194   G4cout << "\n #### proton : CrossSection, MeanFreePath and StoppingPower"
0195          << " for " << material->GetName() << ";\tEnergy cut = " << G4BestUnit(Ecut, "Energy")
0196          << G4endl;
0197
0198   G4cout << "\n Energy \t ionization \t";
0199   G4cout << "\t ionization \t";
0200   G4cout << "\t ionization" << G4endl;
0201
0202   for (G4double Energy = Emin; Energy <= Emax; Energy += dE) {
0203     G4cout << "\n " << G4BestUnit(Energy, "Energy") << "\t"
0204            << G4BestUnit(ioni->ComputeCrossSectionPerAtom(prot, Energy, Z, A, Ecut), "Surface")
0205            << "\t \t"
0206            << G4BestUnit(ioni->ComputeMeanFreePath(prot, Energy, material, Ecut), "Length")
0207            << "\t \t"
0208            << G4BestUnit(ioni->ComputeDEDXPerVolume(material, prot, Energy, Ecut), "Energy/Length");
0209   }
0210
0211   G4cout << G4endl;
0212
0213   // low energy : Bragg Model
0214   ioni = new G4BraggModel(prot);
0215   ioni->Initialise(prot, cuts);
0216
0217   // compute CrossSection per atom and MeanFreePath
0218   //
0219   Emin = 1.1 * keV;
0220   Emax = 2.01 * MeV;
0221   dE = 300 * keV;
0222   Ecut = 10 * keV;
0223
0224   G4cout << "\n #### proton : low energy model (Bragg) "
0225          << ";\tEnergy cut = " << G4BestUnit(Ecut, "Energy") << G4endl;
0226
0227   G4cout << "\n Energy \t ionization \t";
0228   G4cout << "\t ionization \t";
0229   G4cout << "\t ionization" << G4endl;
0230
0231   for (G4double Energy = Emin; Energy <= Emax; Energy += dE) {
0232     G4cout << "\n " << G4BestUnit(Energy, "Energy") << "\t"
0233            << G4BestUnit(ioni->ComputeCrossSectionPerAtom(prot, Energy, Z, A, Ecut), "Surface")
0234            << "\t \t"
0235            << G4BestUnit(ioni->ComputeMeanFreePath(prot, Energy, material, Ecut), "Length")
0236            << "\t \t"
0237            << G4BestUnit(ioni->ComputeDEDXPerVolume(material, prot, Energy, Ecut), "Energy/Length");
0238   }
0239
0240   G4cout << G4endl;
0241
0242   // initialise muon processes (models)
0243   //
0244   ioni = new G4MuBetheBlochModel();
0245   brem = new G4MuBremsstrahlungModel();
0246   G4VEmModel* pair = new G4MuPairProductionModel();
0247   ioni->Initialise(muon, cuts);
0248   brem->Initialise(muon, cuts);
0249   pair->Initialise(muon, cuts);
0250
0251   // compute CrossSection per atom and MeanFreePath
0252   //
0253   Emin = 1.01 * GeV;
0254   Emax = 101.01 * GeV;
0255   dE = 10 * GeV;
0256   Ecut = 10 * MeV;
0257
0258   G4cout << "\n ####muon: CrossSection and MeanFreePath for " << material->GetName()
0259          << ";\tEnergy cut = " << G4BestUnit(Ecut, "Energy") << G4endl;
0260
0261   G4cout << "\n Energy \t ionization \t bremsstra \t pair_prod \t";
0262   G4cout << "\t ionization \t bremsstra \t pair_prod" << G4endl;
0263
0264   for (G4double Energy = Emin; Energy <= Emax; Energy += dE) {
0265     G4cout << "\n " << G4BestUnit(Energy, "Energy") << "\t"
0266            << G4BestUnit(ioni->ComputeCrossSectionPerAtom(muon, Energy, Z, A, Ecut), "Surface")
0267            << "\t"
0268            << G4BestUnit(brem->ComputeCrossSectionPerAtom(muon, Energy, Z, A, Ecut), "Surface")
0269            << "\t"
0270            << G4BestUnit(pair->ComputeCrossSectionPerAtom(muon, Energy, Z, A, Ecut), "Surface")
0271            << "\t \t"
0272            << G4BestUnit(ioni->ComputeMeanFreePath(muon, Energy, material, Ecut), "Length") << "\t"
0273            << G4BestUnit(brem->ComputeMeanFreePath(muon, Energy, material, Ecut), "Length") << "\t"
0274            << G4BestUnit(pair->ComputeMeanFreePath(muon, Energy, material, Ecut), "Length");
0275   }
0276
0277   G4cout << G4endl;
0278
0279   G4cout << "\n ####muon: StoppingPower for " << material->GetName()
0280          << ";\tEnergy cut = " << G4BestUnit(Ecut, "Energy") << G4endl;
0281
0282   G4cout << "\n Energy \t ionization \t bremsstra \t pair_prod \t" << G4endl;
0283
0284   for (G4double Energy = Emin; Energy <= Emax; Energy += dE) {
0285     G4cout << "\n " << G4BestUnit(Energy, "Energy") << "\t"
0286            << G4BestUnit(ioni->ComputeDEDXPerVolume(material, muon, Energy, Ecut), "Energy/Length")
0287            << "\t"
0288            << G4BestUnit(brem->ComputeDEDXPerVolume(material, muon, Energy, Ecut), "Energy/Length")
0289            << "\t"
0290            << G4BestUnit(pair->ComputeDEDXPerVolume(material, muon, Energy, Ecut), "Energy/Length");
0291   }
0292
0293   G4cout << G4endl;
0294   return EXIT_SUCCESS;
0295 }
0296
0297 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......