| 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 // G4GeneralParticleSource class implementatio << 26 ///////////////////////////////////////////////////////////////////////////////
>> 27 //
>> 28 // MODULE: G4GeneralParticleSource.cc
>> 29 //
>> 30 // Version: 2.0
>> 31 // Date: 5/02/04
>> 32 // Author: Fan Lei
>> 33 // Organisation: QinetiQ ltd.
>> 34 // Customer: ESA/ESTEC
>> 35 //
>> 36 // Documentation avaialable at http://reat.space.qinetiq.com/gps
>> 37 // These include:
>> 38 // User Requirement Document (URD)
>> 39 // Software Specification Documents (SSD)
>> 40 // Software User Manual (SUM): on-line version available
>> 41 // Technical Note (TN) on the physics and algorithms
>> 42 //
>> 43 ///////////////////////////////////////////////////////////////////////////////
>> 44 //
>> 45 // CHANGE HISTORY
>> 46 // --------------
>> 47 //
>> 48 // Version 2.0, 05/02/2004, Fan Lei, Created.
>> 49 // After changes to version 1.1 as in Geant4 v6.0
>> 50 // - Mutilple particle source definition
>> 51 // - Re-structured commands
>> 52 // - Split the task into smaller classes
>> 53 //
>> 54 // - old commonds have been retained for backward compatibility, will be
>> 55 // removed in the future.
>> 56 //
>> 57 //
>> 58 ///////////////////////////////////////////////////////////////////////////////
27 // 59 //
28 // Author: Fan Lei, QinetiQ ltd - 05/02/2004 <<
29 // Customer: ESA/ESTEC <<
30 // Version: 2.0 <<
31 // ------------------------------------------- <<
32 <<
33 #include "G4Event.hh" 60 #include "G4Event.hh"
34 #include "Randomize.hh" 61 #include "Randomize.hh"
35 #include "G4GeneralParticleSource.hh" 62 #include "G4GeneralParticleSource.hh"
36 #include "G4SingleParticleSource.hh" <<
37 #include "G4UnitsTable.hh" <<
38 <<
39 #include "G4GeneralParticleSourceData.hh" <<
40 <<
41 #include "G4Threading.hh" <<
42 #include "G4AutoLock.hh" <<
43 <<
44 namespace <<
45 { <<
46 G4Mutex messangerInit = G4MUTEX_INITIALIZER; <<
47 } <<
48 63
49 G4GeneralParticleSource::G4GeneralParticleSour 64 G4GeneralParticleSource::G4GeneralParticleSource()
>> 65 : multiple_vertex(false), flat_sampling(false), weight_change(1.)
50 { 66 {
51 GPSData = G4GeneralParticleSourceData::Insta << 67 sourceVector.clear();
52 // currentSource = GPSData->GetCurrentSource << 68 sourceIntensity.clear();
53 // currentSourceIdx = G4int(GPSData->GetSour << 69 sourceProbability.clear();
54 << 70 currentSource = new G4SingleParticleSource();
55 // Messenger is special, only a worker shoul << 71 sourceVector.push_back(currentSource);
56 // Singleton pattern << 72 sourceIntensity.push_back(1.);
57 // << 73 currentSourceIdx = G4int(sourceVector.size() - 1);
58 theMessenger = G4GeneralParticleSourceMessen << 74 theMessenger = new G4GeneralParticleSourceMessenger(this);
59 << 75 theMessenger->SetParticleGun(currentSource);
60 // Some initialization should be done only o << 76 IntensityNormalization();
61 // <<
62 G4AutoLock l(&messangerInit); <<
63 static G4bool onlyOnce = false; <<
64 if ( !onlyOnce ) <<
65 { <<
66 theMessenger->SetParticleGun(GPSData->GetC <<
67 IntensityNormalization(); <<
68 onlyOnce = true; <<
69 } <<
70 } 77 }
71 << 78
72 G4GeneralParticleSource::~G4GeneralParticleSou 79 G4GeneralParticleSource::~G4GeneralParticleSource()
73 { 80 {
74 theMessenger->Destroy(); << 81 delete theMessenger;
75 } 82 }
76 83
77 void G4GeneralParticleSource::AddaSource(G4dou 84 void G4GeneralParticleSource::AddaSource(G4double aV)
78 { 85 {
79 GPSData->Lock(); << 86 currentSource = new G4SingleParticleSource();
80 << 87 theMessenger->SetParticleGun(currentSource);
81 GPSData->AddASource(aV); << 88 sourceVector.push_back(currentSource);
82 theMessenger->SetParticleGun(GPSData->GetCur << 89 sourceIntensity.push_back(aV);
83 << 90 currentSourceIdx = G4int(sourceVector.size() - 1);
84 // TODO: But do we really normalize here aft <<
85 IntensityNormalization(); 91 IntensityNormalization();
86 <<
87 GPSData->Unlock(); <<
88 } 92 }
89 93
90 void G4GeneralParticleSource::IntensityNormali 94 void G4GeneralParticleSource::IntensityNormalization()
91 { 95 {
92 GPSData->IntensityNormalise(); << 96 G4double total = 0.;
93 normalised=GPSData->Normalised(); << 97 size_t i = 0 ;
94 } << 98 for (i = 0; i < sourceIntensity.size(); i++)
>> 99 total += sourceIntensity[i] ;
>> 100 //
>> 101 sourceProbability.clear();
>> 102 std::vector <G4double> sourceNormalizedIntensity;
>> 103 sourceNormalizedIntensity.clear();
>> 104
>> 105 sourceNormalizedIntensity.push_back(sourceIntensity[0]/total);
>> 106 sourceProbability.push_back(sourceNormalizedIntensity[0]);
>> 107
>> 108 for ( i = 1 ; i < sourceIntensity.size(); i++) {
>> 109 sourceNormalizedIntensity.push_back(sourceIntensity[i]/total);
>> 110 sourceProbability.push_back(sourceNormalizedIntensity[i] + sourceProbability[i-1]);
>> 111 }
>> 112
>> 113 // set source weights here based on sampling scheme (analog/flat) and intensities
>> 114 for ( i = 0 ; i < sourceIntensity.size(); i++) {
>> 115 if (!flat_sampling) {
>> 116 sourceVector[i]->GetBiasRndm()->SetIntensityWeight(1.);
>> 117 } else {
>> 118 sourceVector[i]->GetBiasRndm()->SetIntensityWeight(sourceNormalizedIntensity[i]*sourceIntensity.size());
>> 119 }
>> 120 }
>> 121
>> 122 normalised = true;
>> 123 }
95 124
96 void G4GeneralParticleSource::ListSource() 125 void G4GeneralParticleSource::ListSource()
97 { 126 {
98 G4cout << "The number of particle sources is << 127 G4cout << " The number of particle sources is " << sourceIntensity.size() << G4endl;
99 << GPSData->GetIntensityVectorSize() << 128 for (size_t i = 0 ; i < sourceIntensity.size(); i++)
100 G4cout << " Multiple Vertex sources: " << GP << 129 G4cout << " source " << i << " intensity is " << sourceIntensity[i] << G4endl;
101 G4cout << " Flat Sampling flag: " << GPSData <<
102 const G4int currentIdx = GPSData->GetCurrent <<
103 for(G4int i=0; i<GPSData->GetIntensityVector <<
104 { <<
105 G4cout << "\tsource " << i << " with inten <<
106 << GPSData->GetIntensity(i) << G4en <<
107 const G4SingleParticleSource* thisSrc = GP <<
108 G4cout << " \t\tNum Particles: "<<thisSrc- <<
109 << "; Particle type: " <<
110 << thisSrc->GetParticleDefinition() <<
111 G4cout << " \t\tEnergy: " <<
112 << G4BestUnit(thisSrc->GetEneDist()->GetM <<
113 G4cout << " \t\tDirection: " <<
114 << thisSrc->GetAngDist()->GetDirect <<
115 G4cout << G4BestUnit(thisSrc->GetPosDist() <<
116 << G4endl; <<
117 G4cout << " \t\tAngular Distribution: " <<
118 << thisSrc->GetAngDist()->GetDistTy <<
119 G4cout << " \t\tEnergy Distribution: " <<
120 << thisSrc->GetEneDist()->GetEnergy <<
121 G4cout << " \t\tPosition Distribution Type <<
122 << thisSrc->GetPosDist()->GetPosDis <<
123 G4cout << "; Position Shape: " <<
124 << thisSrc->GetPosDist()->GetPosDis <<
125 } <<
126 <<
127 // Set back previous source <<
128 GPSData->GetCurrentSource(currentIdx); <<
129 } 130 }
130 131
131 void G4GeneralParticleSource::SetCurrentSource 132 void G4GeneralParticleSource::SetCurrentSourceto(G4int aV)
132 { 133 {
133 G4int id = aV; << 134 size_t id = size_t (aV) ;
134 if ( id < GPSData->GetIntensityVectorSize() << 135 if ( id <= sourceIntensity.size() ) {
135 { << 136 currentSourceIdx = aV;
136 // currentSourceIdx = aV; << 137 currentSource = sourceVector[id];
137 // currentSource = GPSData->GetCurrentSour << 138 theMessenger->SetParticleGun(currentSource);
138 theMessenger->SetParticleGun(GPSData->GetC << 139 //
139 } << 140 } else {
140 else << 141 G4cout << " source index is invalid " << G4endl;
141 { << 142 G4cout << " it shall be <= " << sourceIntensity.size() << G4endl;
142 G4ExceptionDescription msg; <<
143 msg << "Trying to set source to index " << <<
144 << GPSData->GetIntensityVectorSize() < <<
145 G4Exception("G4GeneralParticleSoruce::SetC <<
146 FatalException, msg); <<
147 } 143 }
148 } 144 }
149 145
150 void G4GeneralParticleSource::SetCurrentSource 146 void G4GeneralParticleSource::SetCurrentSourceIntensity(G4double aV)
151 { 147 {
152 GPSData->Lock(); << 148 sourceIntensity[currentSourceIdx] = aV;
153 GPSData->SetCurrentSourceIntensity(aV); << 149 normalised = false;
154 GPSData->Unlock(); <<
155 normalised = GPSData->Normalised(); <<
156 } 150 }
157 151
158 void G4GeneralParticleSource::ClearAll() 152 void G4GeneralParticleSource::ClearAll()
159 { 153 {
160 GPSData->ClearSources(); << 154 currentSourceIdx = -1;
161 normalised=GPSData->Normalised(); << 155 currentSource = 0;
>> 156 sourceVector.clear();
>> 157 sourceIntensity.clear();
>> 158 sourceProbability.clear();
162 } 159 }
163 160
164 void G4GeneralParticleSource::DeleteaSource(G4 161 void G4GeneralParticleSource::DeleteaSource(G4int aV)
165 { 162 {
166 G4int id = aV; << 163 size_t id = size_t (aV) ;
167 if ( id <= GPSData->GetIntensityVectorSize() << 164 if ( id <= sourceIntensity.size() ) {
168 { << 165 sourceVector.erase(sourceVector.begin()+aV);
169 GPSData->DeleteASource(aV); << 166 sourceIntensity.erase(sourceIntensity.begin()+aV);
170 normalised=GPSData->Normalised(); << 167 normalised = false ;
171 } << 168 if (currentSourceIdx == aV ) {
172 else << 169 if ( sourceIntensity.size() > 0 ) {
173 { << 170 currentSource = sourceVector[0];
>> 171 currentSourceIdx = 1;
>> 172 } else {
>> 173 currentSource = 0;
>> 174 currentSourceIdx = -1;
>> 175 }
>> 176 }
>> 177 } else {
174 G4cout << " source index is invalid " << G 178 G4cout << " source index is invalid " << G4endl;
175 G4cout << " it shall be <= " << 179 G4cout << " it shall be <= " << sourceIntensity.size() << G4endl;
176 << GPSData->GetIntensityVectorSize( <<
177 } 180 }
178 } << 181 }
179 182
180 void G4GeneralParticleSource::GeneratePrimaryV 183 void G4GeneralParticleSource::GeneratePrimaryVertex(G4Event* evt)
181 { 184 {
182 if (!GPSData->GetMultipleVertex()) << 185 if (!multiple_vertex){
183 { << 186 if (sourceIntensity.size() > 1) {
184 G4SingleParticleSource* currentSource = GP << 187 if (!normalised) IntensityNormalization();
185 if (GPSData->GetIntensityVectorSize() > 1) <<
186 { <<
187 // Try to minimize locks <<
188 if (! normalised ) <<
189 { <<
190 // According to local variable, normal <<
191 // Check with underlying shared resour <<
192 // thread could have already normalize <<
193 GPSData->Lock(); <<
194 G4bool norm = GPSData->Normalised(); <<
195 if (!norm) <<
196 { <<
197 IntensityNormalization(); <<
198 } <<
199 // This takes care of the case in whic <<
200 // is False and the underlying resourc <<
201 normalised = GPSData->Normalised(); <<
202 GPSData->Unlock(); <<
203 } <<
204 G4double rndm = G4UniformRand(); 188 G4double rndm = G4UniformRand();
205 G4int i = 0 ; << 189 size_t i = 0 ;
206 if (! GPSData->GetFlatSampling() ) << 190 if (!flat_sampling) {
207 { << 191 while ( rndm > sourceProbability[i] ) i++;
208 while ( rndm > GPSData->GetSourceProba << 192 (currentSource = sourceVector[i]);
209 currentSource = GPSData->GetCurrentS << 193 } else {
210 } << 194 i = size_t (sourceIntensity.size()*rndm);
211 else << 195 currentSource = sourceVector[i];
212 { <<
213 i = G4int (GPSData->GetIntensityVector <<
214 currentSource = GPSData->GetCurrentSou <<
215 } 196 }
216 } 197 }
217 currentSource->GeneratePrimaryVertex(evt); << 198 currentSource-> GeneratePrimaryVertex(evt);
218 } 199 }
219 else << 200 else {
220 { << 201 for (size_t i = 0; i < sourceIntensity.size(); i++) {
221 for (G4int i = 0; i < GPSData->GetIntensi << 202 sourceVector[i]->GeneratePrimaryVertex(evt);
222 { <<
223 GPSData->GetCurrentSource(i)->GeneratePr <<
224 } 203 }
225 } 204 }
226 } 205 }
227 206