Geant4_10
G4BOptrForceCollision.cc
Go to the documentation of this file.
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 //
26 #include "G4BOptrForceCollision.hh"
28 
32 #include "G4BOptnCloning.hh"
33 
34 #include "G4Step.hh"
35 #include "G4StepPoint.hh"
36 #include "G4VProcess.hh"
37 
38 #include "G4ParticleDefinition.hh"
39 #include "G4ParticleTable.hh"
40 
41 #include "G4SystemOfUnits.hh"
42 
44  : G4VBiasingOperator(name),
45  fFirstProcess(0), fLastProcess(0),
46  fSetup(true)
47 {
48  fSharedForceInteractionOperation = new G4BOptnForceCommonTruncatedExp("SharedForceInteraction");
49  fCloningOperation = new G4BOptnCloning("Cloning");
50  fParticle = G4ParticleTable::GetParticleTable()->FindParticle(particleName);
51 
52  if ( fParticle == 0 )
53  {
55  ed << " Particle `" << particleName << "' not found !" << G4endl;
56  G4Exception(" G4BOptrForceCollision::G4BOptrForceCollision(...)",
57  "BIAS.GEN.07",
59  ed);
60  }
61 }
62 
64  : G4VBiasingOperator(name),
65  fFirstProcess(0), fLastProcess(0),
66  fSetup(true)
67 {
68  fSharedForceInteractionOperation = new G4BOptnForceCommonTruncatedExp("SharedForceInteraction");
69  fCloningOperation = new G4BOptnCloning("Cloning");
70  fParticle = particle;
71 }
72 
74 {
75  for ( std::map< const G4BiasingProcessInterface*, G4BOptnForceFreeFlight* >::iterator it = fFreeFlightOperations.begin() ;
76  it != fFreeFlightOperations.end() ;
77  it++ ) delete (*it).second;
78  delete fSharedForceInteractionOperation;
79  delete fCloningOperation;
80 }
81 
82 
83 G4VBiasingOperation* G4BOptrForceCollision::ProposeOccurenceBiasingOperation(const G4Track* track, const G4BiasingProcessInterface* callingProcess)
84 {
85  if ( track->GetDefinition() != fParticle ) return 0;
86 
87  // -- start by remembering processes under biasing, and create needed biasing operations:
88  // -- ( Might consider moving this in a less called method. )
89  if ( fSetup )
90  {
91  if ( ( fFirstProcess == 0 ) && ( callingProcess->GetIsFirstPostStepGPILInterface() ) ) fFirstProcess = callingProcess;
92  if ( fLastProcess == 0 )
93  {
94  fProcesses.push_back(callingProcess);
95  G4String operationName = "FreeFlight-"+callingProcess->GetWrappedProcess()->GetProcessName();
96  fFreeFlightOperations[callingProcess] = new G4BOptnForceFreeFlight(operationName);
97  if ( callingProcess->GetIsLastPostStepGPILInterface() )
98  {
99  fLastProcess = callingProcess;
100  fSetup = false;
101  }
102  }
103  }
104 
105 
106  // -- Send force interaction operation to the callingProcess:
107  // ----------------------------------------------------------
108  // -- at this level, a copy of the track entering the volume was
109  // -- generated (borned) earlier. This copy will make the forced
110  // -- interaction in the volume.
111  if ( GetBirthOperation( track ) == fCloningOperation )
112  {
113  // -- forced interaction already occured, abort.
114  // -- [Note that if ones would redo a forced interaction, it
115  // -- would require an other cloning before, if one wants to conserve
116  // -- the weight.]
117  if ( fSharedForceInteractionOperation->GetInteractionOccured() ) return 0;
118 
119  if ( callingProcess == fFirstProcess )
120  {
121  // -- first step of cloned track, initialize the forced interaction operation:
122  if ( track->GetCurrentStepNumber() == 1 ) fSharedForceInteractionOperation->Initialize( track );
123  else
124  {
125  if ( fSharedForceInteractionOperation->GetInitialMomentum() != track->GetMomentum() )
126  {
127  // -- means that some other physics process, not under control of the forced interaction operation,
128  // -- has occured, need to re-initialize the operation as distance to boundary has changed.
129  // -- [ Note the re-initialization is only possible for a Markovian law. ]
130  fSharedForceInteractionOperation->Initialize( track );
131  }
132  else
133  {
134  // -- means that some other non-physics process (biasing or not, like step limit), has occured,
135  // -- but track conserves its momentum direction, only need to reduced the maximum distance for
136  // -- forced interaction.
137  // -- [ Note the update is only possible for a Markovian law. ]
138  fSharedForceInteractionOperation->UpdateForStep( track->GetStep() );
139  }
140  }
141  }
142 
143  // -- Sanity check : it may happen in limit cases that distance to out is zero,
144  // -- weight would be infinite in this case. Abort forced interaction.
145  if ( fSharedForceInteractionOperation->GetMaximumDistance() < DBL_MIN ) return 0;
146 
147  // -- conditions to apply forced interaction are met, set up physics:
148  // -- Collects well-defined cross-sections...
149  G4double currentInteractionLength = callingProcess->GetWrappedProcess()->GetCurrentInteractionLength();
150  G4VBiasingOperation* operationToReturn = 0;
151  if ( currentInteractionLength < DBL_MAX/10. )
152  {
153  fSharedForceInteractionOperation->AddCrossSection( callingProcess->GetWrappedProcess(), 1.0/currentInteractionLength );
154  operationToReturn = fSharedForceInteractionOperation;
155  }
156  // -- ... if current process is the last one, cross-section collection is finished.
157  // -- Operation is ready to sample the force interaction law and to randomly select
158  // -- the process to be applied.
159  if ( ( callingProcess == fLastProcess ) &&
160  ( fSharedForceInteractionOperation->GetCommonTruncatedExpLaw()->GetNumberOfSharing() > 0 ) )
161  fSharedForceInteractionOperation->Sample();
162 
163  // -- back to current process (last one or not), if meaningful cross section, bias it returning the force operation:
164  return operationToReturn;
165 
166  } // -- end of " if ( GetBirthOperation( track ) == fCloningOperation ) "
167 
168 
169 
170  // -- Send force free flight to the callingProcess:
171  // ------------------------------------------------
172  // -- At this point a track cloning happened in the previous step,
173  // -- we request the continuing/current track to be forced flight.
174  // -- Note this track will fly with 0.0 weight during its forced flight:
175  // -- it is to forbid the double counting with the force interaction track.
176  // -- Its weight is restored at the end of its free flight, this weight
177  // -- being its initial weight * the weight for the free flight travel,
178  // -- this last one being per process. The initial weight is common, and is
179  // -- arbitrary asked to the first operation to take care of it.
180  if ( fPreviousOperationApplied == fCloningOperation )
181  {
182  G4BOptnForceFreeFlight* operation = fFreeFlightOperations[callingProcess];
183  if ( callingProcess->GetWrappedProcess()->GetCurrentInteractionLength() < DBL_MAX/10. )
184  {
185  // -- the initial track weight will be restored only by the first DoIt free flight:
186  operation->ResetInitialTrackWeight(fInitialTrackWeight);
187  return operation;
188  }
189  }
190  // -- If forced flight was already requested for the calling process, this flight
191  // -- may have been interupted by some other non-physics process : request
192  // -- process to continue with same forced flight:
193  // -- ** ?? ** Incorrect here in case of flight interrupted by a *physics process*. Ie case
194  // -- ** ?? ** of a subset of physics processes forced is not treated correctly here ** ?? **
195  else if ( callingProcess->GetPreviousOccurenceBiasingOperation() == fFreeFlightOperations[callingProcess] )
196  return fFreeFlightOperations[callingProcess];
197 
198  // -- other cases here: particle appearing in the volume by some
199  // -- previous interaction : we decide to not bias these.
200  return 0;
201 
202 }
203 
204 
205 G4VBiasingOperation* G4BOptrForceCollision::ProposeNonPhysicsBiasingOperation(const G4Track* track,
207 {
208  if ( track->GetDefinition() != fParticle ) return 0;
209 
210  // -- bias only tracks entering the volume.
211  // -- A "cloning" is done:
212  // -- - the primary will be forced flight under a zero weight up to volume exit,
213  // -- where the weight will be restored with proper weight for free flight
214  // -- - the clone will be forced to interact in the volume.
215  if ( track->GetStep()->GetPreStepPoint()->GetStepStatus() == fGeomBoundary )
216  {
217  fSharedForceInteractionOperation->SetInteractionOccured( false );
218  fInitialTrackWeight = track->GetWeight();
219  fCloningOperation->SetCloneWeights(0.0, fInitialTrackWeight);
220  return fCloningOperation;
221  }
222 
223  // -- for all other cases: does nothing
224  return 0;
225 }
226 
227 
229 {
230  fPreviousOperationApplied = 0;
231 }
232 
233 
235 {
236  fPreviousOperationApplied = 0;
237  // -- clean-up track for what happens with this operator:
238  ForgetTrack ( track );
239 }
240 
241 
243  G4VBiasingOperation* operationApplied, const G4VParticleChange* particleChangeProduced )
244 {
245  fPreviousOperationApplied = operationApplied;
246  if ( operationApplied == fCloningOperation )
247  RememberSecondaries( callingProcess, operationApplied, particleChangeProduced );
248 }
virtual void OperationApplied(const G4BiasingProcessInterface *callingProcess, G4BiasingAppliedCase biasingCase, G4VBiasingOperation *operationApplied, const G4VParticleChange *particleChangeProduced)
G4ParticleDefinition * GetDefinition() const
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
virtual void StartTracking(const G4Track *track)
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4StepStatus GetStepStatus() const
void AddCrossSection(const G4VProcess *, G4double)
const XML_Char * name
Definition: expat.h:151
const G4ThreeVector & GetInitialMomentum() const
const G4Step * GetStep() const
G4BOptrForceCollision(G4String particleToForce, G4String name="ForceCollision")
G4BiasingAppliedCase
G4VProcess * GetWrappedProcess() const
G4bool GetIsFirstPostStepGPILInterface(G4bool physOnly=true) const
G4StepPoint * GetPreStepPoint() const
void ResetInitialTrackWeight(G4double w)
G4int GetCurrentStepNumber() const
G4bool GetIsLastPostStepGPILInterface(G4bool physOnly=true) const
G4double GetCurrentInteractionLength() const
Definition: G4VProcess.hh:462
const G4String & GetProcessName() const
Definition: G4VProcess.hh:408
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
void ForgetTrack(const G4Track *track)
void SetCloneWeights(G4double clone1Weight, G4double clone2Weight)
G4ThreeVector GetMomentum() const
static G4ParticleTable * GetParticleTable()
void RememberSecondaries(const G4BiasingProcessInterface *callingProcess, const G4VBiasingOperation *operationApplied, const G4VParticleChange *particleChangeProduced)
#define DBL_MIN
Definition: templates.hh:75
G4double GetWeight() const
#define G4endl
Definition: G4ios.hh:61
G4ILawCommonTruncatedExp * GetCommonTruncatedExpLaw()
double G4double
Definition: G4Types.hh:76
#define DBL_MAX
Definition: templates.hh:83
const G4VBiasingOperation * GetBirthOperation(const G4Track *)
G4VBiasingOperation * GetPreviousOccurenceBiasingOperation() const
virtual void ExitBiasing(const G4Track *, const G4BiasingProcessInterface *)