Hadr08/src/BiasingOperation.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 BiasingOperation.cc
0027 /// \brief Implementation of the BiasingOperation class
0028 //
0029 //
0030
0031 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0032 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0033
0034 #include "BiasingOperation.hh"
0035
0036 #include "G4BGGNucleonInelasticXS.hh"
0037 #include "G4BGGPionInelasticXS.hh"
0038 #include "G4BiasingProcessInterface.hh"
0039 #include "G4CascadeInterface.hh"
0040 #include "G4CrossSectionDataStore.hh"
0041 #include "G4ExcitedStringDecay.hh"
0042 #include "G4FTFModel.hh"
0043 #include "G4GeneratorPrecompoundInterface.hh"
0044 #include "G4HadronInelasticProcess.hh"
0045 #include "G4HadronicParameters.hh"
0046 #include "G4INCLXXInterface.hh"
0047 #include "G4LundStringFragmentation.hh"
0048 #include "G4NeutronInelasticXS.hh"
0049 #include "G4TheoFSGenerator.hh"
0050 #include "G4VParticleChange.hh"
0051
0052 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0053
0054 BiasingOperation::BiasingOperation(G4String name) : G4VBiasingOperation(name)
0055 {
0056   // Create the inelastic processes for  p , n , pi+ , pi-
0057   fProtonInelasticProcess = new G4HadronInelasticProcess("protonInelastic", G4Proton::Definition());
0058   fNeutronInelasticProcess =
0059     new G4HadronInelasticProcess("neutronInelastic", G4Neutron::Definition());
0060   fPionPlusInelasticProcess =
0061     new G4HadronInelasticProcess("pi+Inelastic", G4PionPlus::Definition());
0062   fPionMinusInelasticProcess =
0063     new G4HadronInelasticProcess("pi-Inelastic", G4PionMinus::Definition());
0064
0065   // Set the energy ranges
0066   const G4double maxBERT = 41.0 * CLHEP::MeV;
0067   const G4double maxINCLXX = 12.0 * CLHEP::GeV;
0068   const G4double minINCLXX = 40.0 * CLHEP::MeV;
0069   const G4double minFTFP = 3.0 * CLHEP::GeV;
0070   const G4double maxFTFP = G4HadronicParameters::Instance()->GetMaxEnergy();
0071
0072   // Create the hadronic models (to replace FTFP_BERT with "FTFP_INCLXX",
0073   // keeping the same energy ranges for the transition between models).
0074   // Notice that it is better to create the models here from scratch,
0075   // instead of reusing the existing ones, because we might pick up the
0076   // existing ones associated to the wrong particles...
0077   // --- FTFP model ---
0078   G4FTFModel* theStringModel = new G4FTFModel;
0079   G4LundStringFragmentation* theLund = new G4LundStringFragmentation;
0080   G4ExcitedStringDecay* theStringDecay = new G4ExcitedStringDecay(theLund);
0081   theStringModel->SetFragmentationModel(theStringDecay);
0082   G4GeneratorPrecompoundInterface* thePrecoInterface = new G4GeneratorPrecompoundInterface;
0083   G4TheoFSGenerator* theHighEnergyModel = new G4TheoFSGenerator("FTFP");
0084   theHighEnergyModel->SetHighEnergyGenerator(theStringModel);
0085   theHighEnergyModel->SetTransport(thePrecoInterface);
0086   theHighEnergyModel->SetMinEnergy(minFTFP);
0087   theHighEnergyModel->SetMaxEnergy(maxFTFP);
0088   // Bertini : create a new model to be used below INCLXX limit
0089   G4CascadeInterface* theBertiniModel = new G4CascadeInterface();
0090   theBertiniModel->SetMinEnergy(0.0);
0091   theBertiniModel->SetMaxEnergy(maxBERT);
0092   // --- Preco ---
0093   // --- INCLXX model ---
0094   // The instance for nucleons:
0095   G4INCLXXInterface* theInclxxModel = new G4INCLXXInterface();
0096   theInclxxModel->SetMinEnergy(minINCLXX);
0097   theInclxxModel->SetMaxEnergy(maxINCLXX);  // Use the same as for FTFP_BERT
0098
0099   // Register the models
0100   fProtonInelasticProcess->RegisterMe(theHighEnergyModel);
0101   fProtonInelasticProcess->RegisterMe(theInclxxModel);
0102   fProtonInelasticProcess->RegisterMe(theBertiniModel);
0103   fNeutronInelasticProcess->RegisterMe(theHighEnergyModel);
0104   fNeutronInelasticProcess->RegisterMe(theInclxxModel);
0105   fNeutronInelasticProcess->RegisterMe(theBertiniModel);
0106   fPionPlusInelasticProcess->RegisterMe(theHighEnergyModel);
0107   fPionPlusInelasticProcess->RegisterMe(theInclxxModel);
0108   fPionPlusInelasticProcess->RegisterMe(theBertiniModel);
0109   fPionMinusInelasticProcess->RegisterMe(theHighEnergyModel);
0110   fPionMinusInelasticProcess->RegisterMe(theInclxxModel);
0111   fPionMinusInelasticProcess->RegisterMe(theBertiniModel);
0112
0113   G4VCrossSectionDataSet* theProtonXSdata = new G4BGGNucleonInelasticXS(G4Proton::Definition());
0114   theProtonXSdata->BuildPhysicsTable(*(G4Proton::Definition()));
0115   fProtonInelasticProcess->AddDataSet(theProtonXSdata);
0116   G4VCrossSectionDataSet* theNeutronXSdata = new G4NeutronInelasticXS;
0117   theNeutronXSdata->BuildPhysicsTable(*(G4Neutron::Definition()));
0118   fNeutronInelasticProcess->AddDataSet(theNeutronXSdata);
0119   G4VCrossSectionDataSet* thePionPlusXSdata = new G4BGGPionInelasticXS(G4PionPlus::Definition());
0120   thePionPlusXSdata->BuildPhysicsTable(*(G4PionPlus::Definition()));
0121   fPionPlusInelasticProcess->AddDataSet(thePionPlusXSdata);
0122   G4VCrossSectionDataSet* thePionMinusXSdata = new G4BGGPionInelasticXS(G4PionMinus::Definition());
0123   thePionMinusXSdata->BuildPhysicsTable(*(G4PionMinus::Definition()));
0124   fPionMinusInelasticProcess->AddDataSet(thePionMinusXSdata);
0125 }
0126
0127 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0128
0129 BiasingOperation::~BiasingOperation() {}
0130
0131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0132
0133 G4VParticleChange* BiasingOperation::ApplyFinalStateBiasing(const G4BiasingProcessInterface*,
0134                                                             const G4Track* track,
0135                                                             const G4Step* step, G4bool&)
0136 {
0137   if (track->GetParticleDefinition() == G4Proton::Definition()) {
0138     auto particle = track->GetDynamicParticle();
0139     auto material = track->GetMaterial();
0140     fProtonInelasticProcess->GetCrossSectionDataStore()->ComputeCrossSection(particle, material);
0141     return fProtonInelasticProcess->PostStepDoIt(*track, *step);
0142   }
0143   else if (track->GetParticleDefinition() == G4Neutron::Definition()) {
0144     auto particle = track->GetDynamicParticle();
0145     auto material = track->GetMaterial();
0146     fNeutronInelasticProcess->GetCrossSectionDataStore()->ComputeCrossSection(particle, material);
0147     return fNeutronInelasticProcess->PostStepDoIt(*track, *step);
0148   }
0149   else if (track->GetParticleDefinition() == G4PionPlus::Definition()) {
0150     auto particle = track->GetDynamicParticle();
0151     auto material = track->GetMaterial();
0152     fPionPlusInelasticProcess->GetCrossSectionDataStore()->ComputeCrossSection(particle, material);
0153     return fPionPlusInelasticProcess->PostStepDoIt(*track, *step);
0154   }
0155   else if (track->GetParticleDefinition() == G4PionMinus::Definition()) {
0156     auto particle = track->GetDynamicParticle();
0157     auto material = track->GetMaterial();
0158     fPionMinusInelasticProcess->GetCrossSectionDataStore()->ComputeCrossSection(particle, material);
0159     return fPionMinusInelasticProcess->PostStepDoIt(*track, *step);
0160   }
0161   else {
0162     G4cerr << "ERROR in BiasingOperation::ApplyFinalStateBiasing : unexpected particle = "
0163            << track->GetParticleDefinition()->GetParticleName() << G4endl;
0164     return 0;
0165   }
0166 }
0167
0168 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......