| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
>> 26 // $Id$
26 // 27 //
27 /// \file GB03BOptnSplitOrKillOnBoundary.cc 28 /// \file GB03BOptnSplitOrKillOnBoundary.cc
28 /// \brief Implementation of the GB03BOptnSpli 29 /// \brief Implementation of the GB03BOptnSplitOrKillOnBoundary class
29 30
30 #include "GB03BOptnSplitOrKillOnBoundary.hh" <<
31 <<
32 #include "Randomize.hh" 31 #include "Randomize.hh"
>> 32 #include "GB03BOptnSplitOrKillOnBoundary.hh"
33 33
34 //....oooOO0OOooo........oooOO0OOooo........oo 34 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
35 35
36 GB03BOptnSplitOrKillOnBoundary::GB03BOptnSplit 36 GB03BOptnSplitOrKillOnBoundary::GB03BOptnSplitOrKillOnBoundary(G4String name)
37 : G4VBiasingOperation(name), fParticleChange << 37 : G4VBiasingOperation(name),
>> 38 fParticleChange(),
>> 39 fParticleChangeForNothing()
38 {} 40 {}
39 41
40 //....oooOO0OOooo........oooOO0OOooo........oo 42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
41 43
42 GB03BOptnSplitOrKillOnBoundary::~GB03BOptnSpli << 44 GB03BOptnSplitOrKillOnBoundary::~GB03BOptnSplitOrKillOnBoundary()
>> 45 {}
43 46
44 //....oooOO0OOooo........oooOO0OOooo........oo 47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
45 48
46 G4double GB03BOptnSplitOrKillOnBoundary::Dista << 49 G4double GB03BOptnSplitOrKillOnBoundary::
47 << 50 DistanceToApplyOperation( const G4Track*,
>> 51 G4double,
>> 52 G4ForceCondition* condition)
48 { 53 {
49 // -- return "infinite" distance for interac 54 // -- return "infinite" distance for interaction, but asks for GenerateBiasingFinalState
50 // -- being called anyway at the end of the 55 // -- being called anyway at the end of the step, by returning the "Forced" condition
51 // -- flag. 56 // -- flag.
52 *condition = Forced; 57 *condition = Forced;
53 return DBL_MAX; 58 return DBL_MAX;
54 } 59 }
55 60
56 //....oooOO0OOooo........oooOO0OOooo........oo 61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
57 62
58 G4VParticleChange* GB03BOptnSplitOrKillOnBound << 63 G4VParticleChange*
59 << 64 GB03BOptnSplitOrKillOnBoundary::GenerateBiasingFinalState( const G4Track* track,
>> 65 const G4Step* step )
60 { 66 {
>> 67
61 // Check if step is limited by the geometry: 68 // Check if step is limited by the geometry: as we attached the biasing operator
62 // to the absorber layer, this volume bounda 69 // to the absorber layer, this volume boundary is the one of the absorber.
63 // (check of current step # of track is inel 70 // (check of current step # of track is inelegant, but is to fix a "feature"
64 // that a cloned track can wrongly be seen i 71 // that a cloned track can wrongly be seen in the wrong volume, because of numerical
65 // precision issue. In this case it makes a 72 // precision issue. In this case it makes a tiny step, which should be disregarded).
66 if ((step->GetPostStepPoint()->GetStepStatus << 73 if ( ( step->GetPostStepPoint()->GetStepStatus() == fGeomBoundary ) &&
67 && (track->GetCurrentStepNumber() != 1)) << 74 ( track->GetCurrentStepNumber() != 1 ) )
68 { << 75 {
69 // -- Before deciding for killing or split << 76
70 // -- the technique or not: << 77 // -- Before deciding for killing or splitting, we make decision on applying
71 G4double trial = G4UniformRand(); // -- N << 78 // -- the technique or not:
72 // -- e << 79 G4double trial = G4UniformRand(); // -- Note: G4UniformRand() is thread-safe
73 if (trial <= fApplyProbability) { << 80 // -- engine for random numbers
74 // -- Get z component of track, to see i << 81 if ( trial <= fApplyProbability )
75 G4double pz = track->GetMomentum().z(); << 82 {
76 << 83 // -- Get z component of track, to see if it moves forward or backward:
77 if (pz > 0.0) { << 84 G4double pz = track->GetMomentum().z();
78 // ----------------------------------- << 85
79 // Here, we are moving "forward". We d << 86 if ( pz > 0.0 )
80 // ----------------------------------- << 87 {
81 << 88 // -------------------------------------------------
82 // Get track weight: << 89 // Here, we are moving "forward". We do "splitting":
83 G4double initialWeight = track->GetWei << 90 // -------------------------------------------------
84 // Define the tracks weight: << 91
85 G4double weightOfTrack = initialWeight << 92 // Get track weight:
86 << 93 G4double initialWeight = track->GetWeight();
87 // The "particle change" is the object << 94 // Define the tracks weight:
88 // the tracking the update of the prim << 95 G4double weightOfTrack = initialWeight/fSplittingFactor;
89 // secondary tracks. << 96
90 fParticleChange.Initialize(*track); << 97 // The "particle change" is the object to be used to communicate to
91 << 98 // the tracking the update of the primary state and/or creation
92 // ask currect track weight to be chan << 99 // secondary tracks.
93 fParticleChange.ProposeParentWeight(we << 100 fParticleChange.Initialize(*track);
94 << 101
95 // Now make clones of this track (this << 102 // ask currect track weight to be changed to new value:
96 // we will then have the primary and N << 103 fParticleChange.ProposeParentWeight( weightOfTrack );
97 // splitting by a factor N: << 104
98 fParticleChange.SetNumberOfSecondaries << 105 // Now make clones of this track (this is the actual splitting):
99 for (G4int iSplit = 1; iSplit < fSplit << 106 // we will then have the primary and N-1 clones of it, hence the
100 G4Track* clone = new G4Track(*track) << 107 // splitting by a factor N:
101 clone->SetWeight(weightOfTrack); << 108 fParticleChange.SetNumberOfSecondaries( fSplittingFactor-1 );
102 fParticleChange.AddSecondary(clone); << 109 for ( G4int iSplit = 1 ; iSplit < fSplittingFactor ; iSplit++ )
103 } << 110 {
104 fParticleChange.SetSecondaryWeightByPr << 111 G4Track* clone = new G4Track( *track );
105 // -- take it as is ;) [though not nec << 112 clone->SetWeight( weightOfTrack );
106 << 113 fParticleChange.AddSecondary( clone );
107 // this new final state is returned to << 114 }
108 return &fParticleChange; << 115 fParticleChange.SetSecondaryWeightByProcess(true); // -- tricky
109 } << 116 // -- take it as is ;) [though not necessary here, put for safety]
110 << 117
111 else { << 118 // this new final state is returned to the tracking;
112 // ----------------------------------- << 119 return &fParticleChange;
113 // Here, we are moving backward. We do << 120
114 // roulette, killing 1/fSplittingFacto << 121 }
115 // ----------------------------------- << 122
116 << 123 else
117 // Get track weight: << 124
118 G4double initialWeight = track->GetWei << 125 {
119 << 126 // --------------------------------------------------------------
120 // The "particle change" is the object << 127 // Here, we are moving backward. We do killing, playing a russian
121 // the tracking the update of the prim << 128 // roulette, killing 1/fSplittingFactor of the tracks in average:
122 // secondary tracks. << 129 // --------------------------------------------------------------
123 fParticleChange.Initialize(*track); << 130
124 << 131 // Get track weight:
125 // Shoot a random number (in ]0,1[ seg << 132 G4double initialWeight = track->GetWeight();
126 G4double random = G4UniformRand(); << 133
127 << 134 // The "particle change" is the object to be used to communicate to
128 // Decide to kill, keeping 1/fSplittin << 135 // the tracking the update of the primary state and/or creation
129 G4double survivingProbability = 1.0 / << 136 // secondary tracks.
130 if (random > survivingProbability) { << 137 fParticleChange.Initialize(*track);
131 // We ask for the the track to be ki << 138
132 fParticleChange.ProposeTrackStatus(f << 139 // Shoot a random number (in ]0,1[ segment):
133 } << 140 G4double random = G4UniformRand();
134 else { << 141
135 // In this case, the track survives. << 142 // Decide to kill with 1/fSplittingFactor probability:
136 // to conserve weight among killed a << 143 G4double killingProbability = 1.0/fSplittingFactor;
137 fParticleChange.ProposeParentWeight( << 144 if ( random < killingProbability )
138 } << 145 {
139 << 146 // We ask for the the track to be killed:
140 // this new final state is returned to << 147 fParticleChange.ProposeTrackStatus(fStopAndKill);
141 return &fParticleChange; << 148 }
142 } << 149 else
143 } // -- end of : if ( trial > probaForApp << 150 {
144 } // -- end of : if ( ( step->GetPostStepPo << 151 // In this case, the track survives. We change its weight
145 // << 152 // to conserve weight among killed and survival tracks:
>> 153 fParticleChange.ProposeParentWeight( initialWeight*fSplittingFactor );
>> 154 }
>> 155
>> 156 // this new final state is returned to the tracking;
>> 157 return &fParticleChange;
>> 158 }
>> 159 } // -- end of : if ( trial > probaForApplying )
>> 160 } // -- end of : if ( ( step->GetPostStepPoint()->GetStepStatus() ==
>> 161 // fGeomBoundary ) ...
>> 162
146 163
147 // Here, the step was not limited by the geo 164 // Here, the step was not limited by the geometry (but certainly by a physics
148 // process). We do nothing: ie we make no ch 165 // process). We do nothing: ie we make no change to the current track.
149 fParticleChangeForNothing.Initialize(*track) 166 fParticleChangeForNothing.Initialize(*track);
150 return &fParticleChangeForNothing; 167 return &fParticleChangeForNothing;
>> 168
151 } 169 }
152 170
153 //....oooOO0OOooo........oooOO0OOooo........oo 171 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
154 172