GB01/src/GB01BOptrChangeCrossSection.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 GB01/src/GB01BOptrChangeCrossSection.cc
0027 /// \brief Implementation of the GB01BOptrChangeCrossSection class
0028 //
0029 #include "GB01BOptrChangeCrossSection.hh"
0030
0031 #include "G4BOptnChangeCrossSection.hh"
0032 #include "G4BiasingProcessInterface.hh"
0033 #include "G4InteractionLawPhysical.hh"
0034 #include "G4ParticleDefinition.hh"
0035 #include "G4ParticleTable.hh"
0036 #include "G4VProcess.hh"
0037 #include "Randomize.hh"
0038
0039 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0040
0041 GB01BOptrChangeCrossSection::GB01BOptrChangeCrossSection(G4String particleName, G4String name)
0042   : G4VBiasingOperator(name), fSetup(true)
0043 {
0044   fParticleToBias = G4ParticleTable::GetParticleTable()->FindParticle(particleName);
0045
0046   if (fParticleToBias == 0) {
0047     G4ExceptionDescription ed;
0048     ed << "Particle `" << particleName << "' not found !" << G4endl;
0049     G4Exception("GB01BOptrChangeCrossSection(...)", "exGB01.01", JustWarning, ed);
0050   }
0051 }
0052
0053 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0054
0055 GB01BOptrChangeCrossSection::~GB01BOptrChangeCrossSection()
0056 {
0057   for (std::map<const G4BiasingProcessInterface*, G4BOptnChangeCrossSection*>::iterator it =
0058          fChangeCrossSectionOperations.begin();
0059        it != fChangeCrossSectionOperations.end(); it++)
0060     delete (*it).second;
0061 }
0062
0063 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0064
0065 void GB01BOptrChangeCrossSection::StartRun()
0066 {
0067   // --------------
0068   // -- Setup stage:
0069   // ---------------
0070   // -- Start by collecting processes under biasing, create needed biasing
0071   // -- operations and associate these operations to the processes:
0072   if (fSetup) {
0073     const G4ProcessManager* processManager = fParticleToBias->GetProcessManager();
0074     const G4BiasingProcessSharedData* sharedData =
0075       G4BiasingProcessInterface::GetSharedData(processManager);
0076     if (sharedData)  // -- sharedData tested, as is can happen a user attaches an operator to a
0077                      // -- volume but without defined BiasingProcessInterface processes.
0078     {
0079       for (size_t i = 0; i < (sharedData->GetPhysicsBiasingProcessInterfaces()).size(); i++) {
0080         const G4BiasingProcessInterface* wrapperProcess =
0081           (sharedData->GetPhysicsBiasingProcessInterfaces())[i];
0082         G4String operationName =
0083           "XSchange-" + wrapperProcess->GetWrappedProcess()->GetProcessName();
0084         fChangeCrossSectionOperations[wrapperProcess] =
0085           new G4BOptnChangeCrossSection(operationName);
0086       }
0087     }
0088     fSetup = false;
0089   }
0090 }
0091
0092 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0093
0094 G4VBiasingOperation* GB01BOptrChangeCrossSection::ProposeOccurenceBiasingOperation(
0095   const G4Track* track, const G4BiasingProcessInterface* callingProcess)
0096 {
0097   // -----------------------------------------------------
0098   // -- Check if current particle type is the one to bias:
0099   // -----------------------------------------------------
0100   if (track->GetDefinition() != fParticleToBias) return 0;
0101
0102   // ---------------------------------------------------------------------
0103   // -- select and setup the biasing operation for current callingProcess:
0104   // ---------------------------------------------------------------------
0105   // -- Check if the analog cross-section well defined : for example, the conversion
0106   // -- process for a gamma below e+e- creation threshold has an DBL_MAX interaction
0107   // -- length. Nothing is done in this case (ie, let analog process to deal with the case)
0108   G4double analogInteractionLength =
0109     callingProcess->GetWrappedProcess()->GetCurrentInteractionLength();
0110   if (analogInteractionLength > DBL_MAX / 10.) return 0;
0111
0112   // -- Analog cross-section is well-defined:
0113   G4double analogXS = 1. / analogInteractionLength;
0114
0115   // -- Choose a constant cross-section bias. But at this level, this factor can be made
0116   // -- direction dependent, like in the exponential transform MCNP case, or it
0117   // -- can be chosen differently, depending on the process, etc.
0118   G4double XStransformation = 2.0;
0119
0120   // -- fetch the operation associated to this callingProcess:
0121   G4BOptnChangeCrossSection* operation = fChangeCrossSectionOperations[callingProcess];
0122   // -- get the operation that was proposed to the process in the previous step:
0123   G4VBiasingOperation* previousOperation = callingProcess->GetPreviousOccurenceBiasingOperation();
0124
0125   // -- now setup the operation to be returned to the process: this
0126   // -- consists in setting the biased cross-section, and in asking
0127   // -- the operation to sample its exponential interaction law.
0128   // -- To do this, to first order, the two lines:
0129   //        operation->SetBiasedCrossSection( XStransformation * analogXS );
0130   //        operation->Sample();
0131   // -- are correct and sufficient.
0132   // -- But, to avoid having to shoot a random number each time, we sample
0133   // -- only on the first time the operation is proposed, or if the interaction
0134   // -- occured. If the interaction did not occur for the process in the previous,
0135   // -- we update the number of interaction length instead of resampling.
0136   if (previousOperation == 0) {
0137     operation->SetBiasedCrossSection(XStransformation * analogXS);
0138     operation->Sample();
0139   }
0140   else {
0141     if (previousOperation != operation) {
0142       // -- should not happen !
0143       G4ExceptionDescription ed;
0144       ed << " Logic problem in operation handling !" << G4endl;
0145       G4Exception("GB01BOptrChangeCrossSection::ProposeOccurenceBiasingOperation(...)", "exGB01.02",
0146                   JustWarning, ed);
0147       return 0;
0148     }
0149     if (operation->GetInteractionOccured()) {
0150       operation->SetBiasedCrossSection(XStransformation * analogXS);
0151       operation->Sample();
0152     }
0153     else {
0154       // -- update the 'interaction length' and underneath 'number of interaction lengths'
0155       // -- for past step  (this takes into accout the previous step cross-section value)
0156       operation->UpdateForStep(callingProcess->GetPreviousStepSize());
0157       // -- update the cross-section value:
0158       operation->SetBiasedCrossSection(XStransformation * analogXS);
0159       // -- forces recomputation of the 'interaction length' taking into account above
0160       // -- new cross-section value [tricky : to be improved]
0161       operation->UpdateForStep(0.0);
0162     }
0163   }
0164
0165   return operation;
0166 }
0167
0168 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0169
0170 void GB01BOptrChangeCrossSection::OperationApplied(const G4BiasingProcessInterface* callingProcess,
0171                                                    G4BiasingAppliedCase,
0172                                                    G4VBiasingOperation* occurenceOperationApplied,
0173                                                    G4double, G4VBiasingOperation*,
0174                                                    const G4VParticleChange*)
0175 {
0176   G4BOptnChangeCrossSection* operation = fChangeCrossSectionOperations[callingProcess];
0177   if (operation == occurenceOperationApplied) operation->SetInteractionOccured();
0178 }
0179
0180 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......