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 DirectAccess.cc
0027 /// \brief Main program of the electromagnetic/TestEm0 example
0028
0029 // ------------------------------------------------------------
0030 //
0031 //  To print cross sections per atom and mean free path for simple material
0032 //
0033 #include "G4BetheBlochModel.hh"
0034 #include "G4BetheHeitlerModel.hh"
0035 #include "G4BraggModel.hh"
0036 #include "G4DataVector.hh"
0037 #include "G4Electron.hh"
0038 #include "G4Gamma.hh"
0039 #include "G4KleinNishinaCompton.hh"
0040 #include "G4Material.hh"
0041 #include "G4MollerBhabhaModel.hh"
0042 #include "G4MuBetheBlochModel.hh"
0043 #include "G4MuBremsstrahlungModel.hh"
0044 #include "G4MuPairProductionModel.hh"
0045 #include "G4MuonPlus.hh"
0046 #include "G4NistManager.hh"
0047 #include "G4PEEffectFluoModel.hh"
0048 #include "G4ParticleTable.hh"
0049 #include "G4Positron.hh"
0050 #include "G4Proton.hh"
0051 #include "G4SeltzerBergerModel.hh"
0052 #include "G4SystemOfUnits.hh"
0053 #include "G4UnitsTable.hh"
0054 #include "G4eeToTwoGammaModel.hh"
0055 #include "globals.hh"
0056
0057 int main()
0058 {
0059   G4UnitDefinition::BuildUnitsTable();
0060
0061   G4ParticleDefinition* gamma = G4Gamma::Gamma();
0062   G4ParticleDefinition* posit = G4Positron::Positron();
0063   G4ParticleDefinition* elec = G4Electron::Electron();
0064   G4ParticleDefinition* prot = G4Proton::Proton();
0065   G4ParticleDefinition* muon = G4MuonPlus::MuonPlus();
0066   G4ParticleTable* partTable = G4ParticleTable::GetParticleTable();
0067   partTable->SetReadiness();
0068
0069   G4DataVector cuts;
0070   cuts.push_back(1 * keV);
0071
0072   // define materials
0073   //
0074   G4Material* material = G4NistManager::Instance()->FindOrBuildMaterial("G4_Fe");
0075
0076   G4cout << *(G4Material::GetMaterialTable()) << G4endl;
0077
0078   G4MaterialCutsCouple* couple = new G4MaterialCutsCouple(material);
0079   couple->SetIndex(0);
0080
0081   // work only for simple materials
0082   G4double Z = material->GetZ();
0083   G4double A = material->GetA();
0084
0085   // initialise gamma processes (models)
0086   //
0087   G4VEmModel* phot = new G4PEEffectFluoModel();
0088   G4VEmModel* comp = new G4KleinNishinaCompton();
0089   G4VEmModel* conv = new G4BetheHeitlerModel();
0090   phot->Initialise(gamma, cuts);
0091   comp->Initialise(gamma, cuts);
0092   conv->Initialise(gamma, cuts);
0093
0094   // valid pointer to a couple is needed for this model
0095   phot->SetCurrentCouple(couple);
0096
0097   // compute CrossSection per atom and MeanFreePath
0098   //
0099   G4double Emin = 1.01 * MeV, Emax = 2.01 * MeV, dE = 100 * keV;
0100
0101   G4cout << "\n #### Gamma : CrossSectionPerAtom and MeanFreePath for " << material->GetName()
0102          << G4endl;
0103   G4cout << "\n Energy \t PhotoElec \t Compton \t Conversion \t";
0104   G4cout << "\t PhotoElec \t Compton \t Conversion" << G4endl;
0105
0106   for (G4double Energy = Emin; Energy <= Emax; Energy += dE) {
0107     G4cout << "\n " << G4BestUnit(Energy, "Energy") << "\t"
0108            << G4BestUnit(phot->ComputeCrossSectionPerAtom(gamma, Energy, Z), "Surface") << "\t"
0109            << G4BestUnit(comp->ComputeCrossSectionPerAtom(gamma, Energy, Z), "Surface") << "\t"
0110            << G4BestUnit(conv->ComputeCrossSectionPerAtom(gamma, Energy, Z), "Surface") << "\t \t"
0111            << G4BestUnit(phot->ComputeMeanFreePath(gamma, Energy, material), "Length") << "\t"
0112            << G4BestUnit(comp->ComputeMeanFreePath(gamma, Energy, material), "Length") << "\t"
0113            << G4BestUnit(conv->ComputeMeanFreePath(gamma, Energy, material), "Length");
0114   }
0115
0116   G4cout << G4endl;
0117
0118   // initialise positron annihilation (model)
0119   //
0120   G4VEmModel* anni = new G4eeToTwoGammaModel();
0121   anni->Initialise(posit, cuts);
0122
0123   // compute CrossSection per atom and MeanFreePath
0124   //
0125   Emin = 1.01 * MeV;
0126   Emax = 2.01 * MeV;
0127   dE = 100 * keV;
0128
0129   G4cout << "\n #### e+ annihilation : CrossSectionPerAtom and MeanFreePath"
0130          << " for " << material->GetName() << G4endl;
0131   G4cout << "\n Energy \t e+ annihil \t";
0132   G4cout << "\t e+ annihil" << G4endl;
0133
0134   for (G4double Energy = Emin; Energy <= Emax; Energy += dE) {
0135     G4cout << "\n " << G4BestUnit(Energy, "Energy") << "\t"
0136            << G4BestUnit(anni->ComputeCrossSectionPerAtom(posit, Energy, Z), "Surface") << "\t \t"
0137            << G4BestUnit(anni->ComputeMeanFreePath(posit, Energy, material), "Length");
0138   }
0139
0140   G4cout << G4endl;
0141
0142   // initialise electron processes (models)
0143   //
0144   G4VEmModel* ioni = new G4MollerBhabhaModel();
0145   G4VEmModel* brem = new G4SeltzerBergerModel();
0146   ioni->Initialise(elec, cuts);
0147   brem->Initialise(elec, cuts);
0148
0149   // compute CrossSection per atom and MeanFreePath
0150   //
0151   Emin = 1.01 * MeV;
0152   Emax = 101.01 * MeV;
0153   dE = 10 * MeV;
0154   G4double Ecut = 100 * keV;
0155
0156   G4cout << "\n ####electron: CrossSection, MeanFreePath and StoppingPower"
0157          << " for " << material->GetName() << ";\tEnergy cut = " << G4BestUnit(Ecut, "Energy")
0158          << G4endl;
0159
0160   G4cout << "\n Energy \t ionization \t bremsstra \t";
0161   G4cout << "\t ionization \t bremsstra \t";
0162   G4cout << "\t ionization \t bremsstra" << G4endl;
0163
0164   for (G4double Energy = Emin; Energy <= Emax; Energy += dE) {
0165     G4cout << "\n " << G4BestUnit(Energy, "Energy") << "\t"
0166            << G4BestUnit(ioni->ComputeCrossSectionPerAtom(elec, Energy, Z, A, Ecut), "Surface")
0167            << "\t"
0168            << G4BestUnit(brem->ComputeCrossSectionPerAtom(elec, Energy, Z, A, Ecut), "Surface")
0169            << "\t \t"
0170            << G4BestUnit(ioni->ComputeMeanFreePath(elec, Energy, material, Ecut), "Length") << "\t"
0171            << G4BestUnit(brem->ComputeMeanFreePath(elec, Energy, material, Ecut), "Length")
0172            << "\t \t"
0173            << G4BestUnit(ioni->ComputeDEDXPerVolume(material, elec, Energy, Ecut), "Energy/Length")
0174            << "\t"
0175            << G4BestUnit(brem->ComputeDEDXPerVolume(material, elec, Energy, Ecut), "Energy/Length");
0176   }
0177
0178   G4cout << G4endl;
0179
0180   // initialise proton processes (models)
0181   //
0182   ioni = new G4BetheBlochModel();
0183   ioni->Initialise(prot, cuts);
0184
0185   // compute CrossSection per atom and MeanFreePath
0186   //
0187   Emin = 1.01 * MeV;
0188   Emax = 102.01 * MeV;
0189   dE = 10 * MeV;
0190   Ecut = 100 * keV;
0191
0192   G4cout << "\n #### proton : CrossSection, MeanFreePath and StoppingPower"
0193          << " for " << material->GetName() << ";\tEnergy cut = " << G4BestUnit(Ecut, "Energy")
0194          << G4endl;
0195
0196   G4cout << "\n Energy \t ionization \t";
0197   G4cout << "\t ionization \t";
0198   G4cout << "\t ionization" << G4endl;
0199
0200   for (G4double Energy = Emin; Energy <= Emax; Energy += dE) {
0201     G4cout << "\n " << G4BestUnit(Energy, "Energy") << "\t"
0202            << G4BestUnit(ioni->ComputeCrossSectionPerAtom(prot, Energy, Z, A, Ecut), "Surface")
0203            << "\t \t"
0204            << G4BestUnit(ioni->ComputeMeanFreePath(prot, Energy, material, Ecut), "Length")
0205            << "\t \t"
0206            << G4BestUnit(ioni->ComputeDEDXPerVolume(material, prot, Energy, Ecut), "Energy/Length");
0207   }
0208
0209   G4cout << G4endl;
0210
0211   // low energy : Bragg Model
0212   ioni = new G4BraggModel(prot);
0213   ioni->Initialise(prot, cuts);
0214
0215   // compute CrossSection per atom and MeanFreePath
0216   //
0217   Emin = 1.1 * keV;
0218   Emax = 2.01 * MeV;
0219   dE = 300 * keV;
0220   Ecut = 10 * keV;
0221
0222   G4cout << "\n #### proton : low energy model (Bragg) "
0223          << ";\tEnergy cut = " << G4BestUnit(Ecut, "Energy") << G4endl;
0224
0225   G4cout << "\n Energy \t ionization \t";
0226   G4cout << "\t ionization \t";
0227   G4cout << "\t ionization" << G4endl;
0228
0229   for (G4double Energy = Emin; Energy <= Emax; Energy += dE) {
0230     G4cout << "\n " << G4BestUnit(Energy, "Energy") << "\t"
0231            << G4BestUnit(ioni->ComputeCrossSectionPerAtom(prot, Energy, Z, A, Ecut), "Surface")
0232            << "\t \t"
0233            << G4BestUnit(ioni->ComputeMeanFreePath(prot, Energy, material, Ecut), "Length")
0234            << "\t \t"
0235            << G4BestUnit(ioni->ComputeDEDXPerVolume(material, prot, Energy, Ecut), "Energy/Length");
0236   }
0237
0238   G4cout << G4endl;
0239
0240   // initialise muon processes (models)
0241   //
0242   ioni = new G4MuBetheBlochModel();
0243   brem = new G4MuBremsstrahlungModel();
0244   G4VEmModel* pair = new G4MuPairProductionModel();
0245   ioni->Initialise(muon, cuts);
0246   brem->Initialise(muon, cuts);
0247   pair->Initialise(muon, cuts);
0248
0249   // compute CrossSection per atom and MeanFreePath
0250   //
0251   Emin = 1.01 * GeV;
0252   Emax = 101.01 * GeV;
0253   dE = 10 * GeV;
0254   Ecut = 10 * MeV;
0255
0256   G4cout << "\n ####muon: CrossSection and MeanFreePath for " << material->GetName()
0257          << ";\tEnergy cut = " << G4BestUnit(Ecut, "Energy") << G4endl;
0258
0259   G4cout << "\n Energy \t ionization \t bremsstra \t pair_prod \t";
0260   G4cout << "\t ionization \t bremsstra \t pair_prod" << G4endl;
0261
0262   for (G4double Energy = Emin; Energy <= Emax; Energy += dE) {
0263     G4cout << "\n " << G4BestUnit(Energy, "Energy") << "\t"
0264            << G4BestUnit(ioni->ComputeCrossSectionPerAtom(muon, Energy, Z, A, Ecut), "Surface")
0265            << "\t"
0266            << G4BestUnit(brem->ComputeCrossSectionPerAtom(muon, Energy, Z, A, Ecut), "Surface")
0267            << "\t"
0268            << G4BestUnit(pair->ComputeCrossSectionPerAtom(muon, Energy, Z, A, Ecut), "Surface")
0269            << "\t \t"
0270            << G4BestUnit(ioni->ComputeMeanFreePath(muon, Energy, material, Ecut), "Length") << "\t"
0271            << G4BestUnit(brem->ComputeMeanFreePath(muon, Energy, material, Ecut), "Length") << "\t"
0272            << G4BestUnit(pair->ComputeMeanFreePath(muon, Energy, material, Ecut), "Length");
0273   }
0274
0275   G4cout << G4endl;
0276
0277   G4cout << "\n ####muon: StoppingPower for " << material->GetName()
0278          << ";\tEnergy cut = " << G4BestUnit(Ecut, "Energy") << G4endl;
0279
0280   G4cout << "\n Energy \t ionization \t bremsstra \t pair_prod \t" << G4endl;
0281
0282   for (G4double Energy = Emin; Energy <= Emax; Energy += dE) {
0283     G4cout << "\n " << G4BestUnit(Energy, "Energy") << "\t"
0284            << G4BestUnit(ioni->ComputeDEDXPerVolume(material, muon, Energy, Ecut), "Energy/Length")
0285            << "\t"
0286            << G4BestUnit(brem->ComputeDEDXPerVolume(material, muon, Energy, Ecut), "Energy/Length")
0287            << "\t"
0288            << G4BestUnit(pair->ComputeDEDXPerVolume(material, muon, Energy, Ecut), "Energy/Length");
0289   }
0290
0291   G4cout << G4endl;
0292   return EXIT_SUCCESS;
0293 }
0294
0295 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......