ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GB01BOptrChangeCrossSection.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file GB01BOptrChangeCrossSection.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
28 //
32 
33 #include "G4ParticleDefinition.hh"
34 #include "G4ParticleTable.hh"
35 #include "G4VProcess.hh"
36 
37 #include "Randomize.hh"
38 
40 
41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
42 
44  G4String name)
45  : G4VBiasingOperator(name),
46  fSetup(true)
47 {
49 
50  if ( fParticleToBias == 0 )
51  {
53  ed << "Particle `" << particleName << "' not found !" << G4endl;
54  G4Exception("GB01BOptrChangeCrossSection(...)",
55  "exGB01.01",
57  ed);
58  }
59 }
60 
61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
62 
64 {
65  for ( std::map< const G4BiasingProcessInterface*, G4BOptnChangeCrossSection* >::iterator
68  it++ ) delete (*it).second;
69 }
70 
71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
72 
74 {
75  // --------------
76  // -- Setup stage:
77  // ---------------
78  // -- Start by collecting processes under biasing, create needed biasing
79  // -- operations and associate these operations to the processes:
80  if ( fSetup )
81  {
82  const G4ProcessManager* processManager = fParticleToBias->GetProcessManager();
83  const G4BiasingProcessSharedData* sharedData =
85  if ( sharedData ) // -- sharedData tested, as is can happen a user attaches an operator to a
86  // -- volume but without defined BiasingProcessInterface processes.
87  {
88  for ( size_t i = 0 ; i < (sharedData->GetPhysicsBiasingProcessInterfaces()).size(); i++ )
89  {
90  const G4BiasingProcessInterface* wrapperProcess =
91  (sharedData->GetPhysicsBiasingProcessInterfaces())[i];
92  G4String operationName = "XSchange-" +
93  wrapperProcess->GetWrappedProcess()->GetProcessName();
94  fChangeCrossSectionOperations[wrapperProcess] =
95  new G4BOptnChangeCrossSection(operationName);
96  }
97  }
98  fSetup = false;
99  }
100 }
101 
102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
103 
107  callingProcess)
108 {
109 
110  // -----------------------------------------------------
111  // -- Check if current particle type is the one to bias:
112  // -----------------------------------------------------
113  if ( track->GetDefinition() != fParticleToBias ) return 0;
114 
115 
116  // ---------------------------------------------------------------------
117  // -- select and setup the biasing operation for current callingProcess:
118  // ---------------------------------------------------------------------
119  // -- Check if the analog cross-section well defined : for example, the conversion
120  // -- process for a gamma below e+e- creation threshold has an DBL_MAX interaction
121  // -- length. Nothing is done in this case (ie, let analog process to deal with the case)
122  G4double analogInteractionLength =
123  callingProcess->GetWrappedProcess()->GetCurrentInteractionLength();
124  if ( analogInteractionLength > DBL_MAX/10. ) return 0;
125 
126  // -- Analog cross-section is well-defined:
127  G4double analogXS = 1./analogInteractionLength;
128 
129  // -- Choose a constant cross-section bias. But at this level, this factor can be made
130  // -- direction dependent, like in the exponential transform MCNP case, or it
131  // -- can be chosen differently, depending on the process, etc.
132  G4double XStransformation = 2.0 ;
133 
134  // -- fetch the operation associated to this callingProcess:
135  G4BOptnChangeCrossSection* operation = fChangeCrossSectionOperations[callingProcess];
136  // -- get the operation that was proposed to the process in the previous step:
137  G4VBiasingOperation* previousOperation = callingProcess->GetPreviousOccurenceBiasingOperation();
138 
139  // -- now setup the operation to be returned to the process: this
140  // -- consists in setting the biased cross-section, and in asking
141  // -- the operation to sample its exponential interaction law.
142  // -- To do this, to first order, the two lines:
143  // operation->SetBiasedCrossSection( XStransformation * analogXS );
144  // operation->Sample();
145  // -- are correct and sufficient.
146  // -- But, to avoid having to shoot a random number each time, we sample
147  // -- only on the first time the operation is proposed, or if the interaction
148  // -- occured. If the interaction did not occur for the process in the previous,
149  // -- we update the number of interaction length instead of resampling.
150  if ( previousOperation == 0 )
151  {
152  operation->SetBiasedCrossSection( XStransformation * analogXS );
153  operation->Sample();
154  }
155  else
156  {
157  if ( previousOperation != operation )
158  {
159  // -- should not happen !
161  ed << " Logic problem in operation handling !" << G4endl;
162  G4Exception("GB01BOptrChangeCrossSection::ProposeOccurenceBiasingOperation(...)",
163  "exGB01.02",
164  JustWarning,
165  ed);
166  return 0;
167  }
168  if ( operation->GetInteractionOccured() )
169  {
170  operation->SetBiasedCrossSection( XStransformation * analogXS );
171  operation->Sample();
172  }
173  else
174  {
175  // -- update the 'interaction length' and underneath 'number of interaction lengths'
176  // -- for past step (this takes into accout the previous step cross-section value)
177  operation->UpdateForStep( callingProcess->GetPreviousStepSize() );
178  // -- update the cross-section value:
179  operation->SetBiasedCrossSection( XStransformation * analogXS );
180  // -- forces recomputation of the 'interaction length' taking into account above
181  // -- new cross-section value [tricky : to be improved]
182  operation->UpdateForStep( 0.0 );
183  }
184  }
185 
186  return operation;
187 
188 }
189 
190 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
191 
195  G4VBiasingOperation* occurenceOperationApplied,
196  G4double,
198  const G4VParticleChange* )
199 {
200  G4BOptnChangeCrossSection* operation = fChangeCrossSectionOperations[callingProcess];
201  if ( operation == occurenceOperationApplied ) operation->SetInteractionOccured();
202 }
203 
204 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......