| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
26 // G4WeightWindowAlgorithm implementation <<
27 // 23 //
28 // Author: Michael Dressel (CERN), 2003 << 24 // $Id: G4WeightWindowAlgorithm.cc,v 1.5 2002/10/16 14:29:07 dressel Exp $
>> 25 // GEANT4 tag $Name: geant4-05-00 $
>> 26 //
>> 27 // ----------------------------------------------------------------------
>> 28 // GEANT 4 class source file
>> 29 //
>> 30 // G4WeightWindowAlgorithm.cc
>> 31 //
29 // ------------------------------------------- 32 // ----------------------------------------------------------------------
30 #include "G4WeightWindowAlgorithm.hh" 33 #include "G4WeightWindowAlgorithm.hh"
31 #include "Randomize.hh" 34 #include "Randomize.hh"
32 35
33 G4WeightWindowAlgorithm:: << 36
34 G4WeightWindowAlgorithm(G4double upperLimitFac << 37 G4WeightWindowAlgorithm::G4WeightWindowAlgorithm() :
35 G4double survivalFacto << 38 fUpper(1),
36 G4int maxNumberOfSplit << 39 fLower(1)
37 : fUpperLimitFactor(upperLimitFactor), << 40 {}
38 fSurvivalFactor(survivalFactor), << 41
39 fMaxNumberOfSplits(maxNumberOfSplits) << 42 G4WeightWindowAlgorithm::~G4WeightWindowAlgorithm()
40 { << 43 {}
>> 44
>> 45
>> 46 void G4WeightWindowAlgorithm::SetUpperLimit(G4double Upper){
>> 47 fUpper = Upper;
>> 48 }
>> 49
>> 50 void G4WeightWindowAlgorithm::SetLowerLimit(G4double Lower){
>> 51 fLower = Lower;
41 } 52 }
42 53
43 G4WeightWindowAlgorithm::~G4WeightWindowAlgori <<
44 54
45 G4Nsplit_Weight 55 G4Nsplit_Weight
46 G4WeightWindowAlgorithm::Calculate(G4double in << 56 G4WeightWindowAlgorithm::Calculate(G4double init_w,
47 G4double lo << 57 G4double importance) const {
48 { << 58
49 G4double survivalWeight = lowerWeightBound * << 59
50 G4double upperWeight = lowerWeightBound * fU << 60 G4Nsplit_Weight nw = {1, init_w};
51 << 61 G4double iw = importance * init_w;
52 // initialize return value for case weight i << 62 if (iw>fUpper) {
53 // << 63 // f is the factor by which the weight is greater
54 G4Nsplit_Weight nw; << 64 // than allowed by fUpper
55 nw.fN = 1; << 65 // it is almost the number of coppies to be produced
56 nw.fW = init_w; << 66 G4double f = iw / fUpper;
57 << 67 // calculate new weight
58 if (init_w > upperWeight) << 68 nw.fW/=f;
59 { << 69
60 // splitting << 70 // calculate the number of coppies
61 // << 71 nw.fN = static_cast<G4int>(f);
62 G4double temp_wi_ws = init_w/upperWeight; << 72 // correct the number of coppies in case f is not an integer
63 auto split_i = static_cast<G4int>(temp_wi << 73 if (static_cast<G4double>(nw.fN) != f) {
64 if(split_i != temp_wi_ws) { ++split_i; } << 74 // probabillity p for splitting into nw.fN+1 particles
65 G4double wi_ws = init_w/split_i; << 75 G4double p = f - nw.fN;
66 << 76 // get a random number out of [0,1)
67 // values in case integer mode or in csae << 77 G4double r = G4UniformRand();
68 // mode and the lower number of splits wil << 78 if (r<p) {
69 // << 79 nw.fN++;
70 nw.fN = split_i; << 80 }
71 nw.fW = wi_ws; << 81 }
72 <<
73 //TB if (wi_ws <= fMaxNumberOfSplits) { <<
74 //TB if (wi_ws > int_wi_ws) { <<
75 //TB // double mode <<
76 //TB G4double p2 = wi_ws - int_wi_ws; <<
77 //TB G4double r = G4UniformRand(); <<
78 //TB if (r<p2) { <<
79 //TB nw.fN = int_wi_ws + 1; <<
80 //TB } <<
81 //TB } <<
82 //TB } <<
83 //TB else { <<
84 //TB // fMaxNumberOfSplits < wi_ws <<
85 //TB nw.fN = fMaxNumberOfSplits; <<
86 //TB nw.fW = init_w/fMaxNumberOfSplits; <<
87 //TB } <<
88 82
89 } 83 }
90 else if (init_w < lowerWeightBound) // Russ << 84 else if (iw < fLower) {
91 { << 85 // play russian roulett
92 G4double wi_ws = init_w/survivalWeight; << 86 G4double p = iw / fLower; // survival prob.
93 G4double p = std::max(wi_ws,1./fMaxNumberO <<
94 G4double r = G4UniformRand(); 87 G4double r = G4UniformRand();
95 if (r<p) << 88 if (r>=p) {
96 { << 89 // kill track
97 nw.fW = init_w/p; <<
98 nw.fN = 1; <<
99 } <<
100 else <<
101 { <<
102 nw.fW = 0; <<
103 nw.fN = 0; 90 nw.fN = 0;
104 } << 91 nw.fW = 0;
105 } << 92 }
106 // else do nothing << 93 else {
107 << 94 nw.fW*=1/p; // to be consistant with the survival prob.
>> 95 }
>> 96 }
108 return nw; 97 return nw;
109 } 98 }
110 99
111 100