| 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 // 26 //
27 // ------------------------------------------- 27 // --------------------------------------------------------------
28 // GEANT 4 - Underground Dark Matter Detecto 28 // GEANT 4 - Underground Dark Matter Detector Advanced Example
29 // 29 //
30 // For information related to this code c 30 // For information related to this code contact: Alex Howard
31 // e-mail: alexander.howard@cern.ch 31 // e-mail: alexander.howard@cern.ch
32 // ------------------------------------------- 32 // --------------------------------------------------------------
33 // Comments 33 // Comments
34 // 34 //
35 // Underground Advanced 35 // Underground Advanced
36 // by A. Howard and H. Araujo 36 // by A. Howard and H. Araujo
37 // (27th November 2001) 37 // (27th November 2001)
38 // 38 //
39 // 39 //
40 // ParticleSource program 40 // ParticleSource program
41 // ------------------------------------------- 41 // --------------------------------------------------------------
42 ////////////////////////////////////////////// 42 //////////////////////////////////////////////////////////////////////////////
43 // This particle source is a shortened version 43 // This particle source is a shortened version of G4GeneralParticleSource by
44 // C Ferguson, F Lei & P Truscott (University 44 // C Ferguson, F Lei & P Truscott (University of Southampton / DERA), with
45 // some minor modifications. 45 // some minor modifications.
46 ////////////////////////////////////////////// 46 //////////////////////////////////////////////////////////////////////////////
47 47
48 #include <cmath> <<
49 <<
50 #include "DMXParticleSource.hh" 48 #include "DMXParticleSource.hh"
51 49
52 #include "G4PhysicalConstants.hh" <<
53 #include "G4SystemOfUnits.hh" <<
54 #include "G4PrimaryParticle.hh" 50 #include "G4PrimaryParticle.hh"
55 #include "G4Event.hh" 51 #include "G4Event.hh"
56 #include "Randomize.hh" 52 #include "Randomize.hh"
>> 53 #include <cmath>
57 #include "G4TransportationManager.hh" 54 #include "G4TransportationManager.hh"
58 #include "G4VPhysicalVolume.hh" 55 #include "G4VPhysicalVolume.hh"
59 #include "G4PhysicalVolumeStore.hh" 56 #include "G4PhysicalVolumeStore.hh"
60 #include "G4ParticleTable.hh" 57 #include "G4ParticleTable.hh"
61 #include "G4ParticleDefinition.hh" 58 #include "G4ParticleDefinition.hh"
62 #include "G4IonTable.hh" 59 #include "G4IonTable.hh"
63 #include "G4Ions.hh" 60 #include "G4Ions.hh"
64 #include "G4TrackingManager.hh" 61 #include "G4TrackingManager.hh"
65 #include "G4Track.hh" 62 #include "G4Track.hh"
66 63
67 64
68 DMXParticleSource::DMXParticleSource() { 65 DMXParticleSource::DMXParticleSource() {
69 66
70 NumberOfParticlesToBeGenerated = 1; 67 NumberOfParticlesToBeGenerated = 1;
71 particle_definition = nullptr; << 68 particle_definition = NULL;
72 G4ThreeVector zero(0., 0., 0.); 69 G4ThreeVector zero(0., 0., 0.);
73 particle_momentum_direction = G4ParticleMome 70 particle_momentum_direction = G4ParticleMomentum(1., 0., 0.);
74 particle_energy = 1.0*MeV; 71 particle_energy = 1.0*MeV;
75 particle_position = zero; 72 particle_position = zero;
76 particle_time = 0.0; 73 particle_time = 0.0;
77 particle_polarization = zero; 74 particle_polarization = zero;
78 particle_charge = 0.0; 75 particle_charge = 0.0;
79 76
80 SourcePosType = "Volume"; 77 SourcePosType = "Volume";
81 Shape = "NULL"; 78 Shape = "NULL";
82 halfz = 0.; 79 halfz = 0.;
83 Radius = 0.; 80 Radius = 0.;
84 CentreCoords = zero; 81 CentreCoords = zero;
85 Confine = false; 82 Confine = false;
86 VolName = "NULL"; 83 VolName = "NULL";
87 84
88 AngDistType = "iso"; 85 AngDistType = "iso";
89 MinTheta = 0.; 86 MinTheta = 0.;
90 MaxTheta = pi; 87 MaxTheta = pi;
91 MinPhi = 0.; 88 MinPhi = 0.;
92 MaxPhi = twopi; 89 MaxPhi = twopi;
93 90
94 EnergyDisType = "Mono"; 91 EnergyDisType = "Mono";
95 MonoEnergy = 1*MeV; 92 MonoEnergy = 1*MeV;
96 93
97 verbosityLevel = 0; 94 verbosityLevel = 0;
98 95
99 theMessenger = new DMXParticleSourceMessenge 96 theMessenger = new DMXParticleSourceMessenger(this);
100 gNavigator = G4TransportationManager::GetTra 97 gNavigator = G4TransportationManager::GetTransportationManager()
101 ->GetNavigatorForTracking(); 98 ->GetNavigatorForTracking();
102 } 99 }
103 100
104 DMXParticleSource::~DMXParticleSource() 101 DMXParticleSource::~DMXParticleSource()
105 { 102 {
106 delete theMessenger; 103 delete theMessenger;
107 } 104 }
108 105
109 void DMXParticleSource::SetPosDisType(G4String 106 void DMXParticleSource::SetPosDisType(G4String PosType)
110 { 107 {
111 SourcePosType = PosType; 108 SourcePosType = PosType;
112 } 109 }
113 110
114 void DMXParticleSource::SetPosDisShape(G4Strin 111 void DMXParticleSource::SetPosDisShape(G4String shapeType)
115 { 112 {
116 Shape = shapeType; 113 Shape = shapeType;
117 } 114 }
118 115
119 void DMXParticleSource::SetCentreCoords(G4Thre 116 void DMXParticleSource::SetCentreCoords(G4ThreeVector coordsOfCentre)
120 { 117 {
121 CentreCoords = coordsOfCentre; 118 CentreCoords = coordsOfCentre;
122 } 119 }
123 120
124 void DMXParticleSource::SetHalfZ(G4double zhal 121 void DMXParticleSource::SetHalfZ(G4double zhalf)
125 { 122 {
126 halfz = zhalf; 123 halfz = zhalf;
127 } 124 }
128 125
129 void DMXParticleSource::SetRadius(G4double rad << 126 void DMXParticleSource::SetRadius(G4double rad)
130 { 127 {
131 Radius = radius; << 128 Radius = rad;
132 } 129 }
133 130
134 void DMXParticleSource::ConfineSourceToVolume( 131 void DMXParticleSource::ConfineSourceToVolume(G4String Vname)
135 { 132 {
136 VolName = Vname; 133 VolName = Vname;
137 if(verbosityLevel == 2) G4cout << VolName << 134 if(verbosityLevel == 2) G4cout << VolName << G4endl;
138 135
139 // checks if selected volume exists 136 // checks if selected volume exists
140 G4VPhysicalVolume *tempPV = nullptr; << 137 G4VPhysicalVolume *tempPV = NULL;
141 G4PhysicalVolumeStore *PVStore = nullptr; << 138 G4PhysicalVolumeStore *PVStore = 0;
142 G4String theRequiredVolumeName = VolName; 139 G4String theRequiredVolumeName = VolName;
143 PVStore = G4PhysicalVolumeStore::GetInstance 140 PVStore = G4PhysicalVolumeStore::GetInstance();
>> 141 G4int i = 0;
144 G4bool found = false; 142 G4bool found = false;
145 if(verbosityLevel == 2) G4cout << PVStore->s 143 if(verbosityLevel == 2) G4cout << PVStore->size() << G4endl;
146 144
147 tempPV = PVStore->GetVolume(theRequiredVolum << 145 // recasting required since PVStore->size() is actually a signed int...
148 if (tempPV != nullptr) { found = true; } << 146 while (!found && i<(G4int)PVStore->size())
>> 147 {
>> 148 tempPV = (*PVStore)[i];
>> 149 found = tempPV->GetName() == theRequiredVolumeName;
>> 150 if(verbosityLevel == 2)
>> 151 G4cout << i << " " << " " << tempPV->GetName()
>> 152 << " " << theRequiredVolumeName << " " << found << G4endl;
>> 153 if (!found)
>> 154 {i++;}
>> 155 }
149 156
150 // found = true then the volume exists else 157 // found = true then the volume exists else it doesnt.
151 if(found == true) { 158 if(found == true) {
152 if(verbosityLevel >= 1) 159 if(verbosityLevel >= 1)
153 G4cout << "Volume " << VolName << " exis 160 G4cout << "Volume " << VolName << " exists" << G4endl;
154 Confine = true; 161 Confine = true;
155 } 162 }
156 else if(VolName=="NULL") 163 else if(VolName=="NULL")
157 Confine = false; 164 Confine = false;
158 else { 165 else {
159 G4cout << " **** Error: Volume does not ex 166 G4cout << " **** Error: Volume does not exist **** " << G4endl;
160 G4cout << " Ignoring confine condition" << 167 G4cout << " Ignoring confine condition" << G4endl;
161 VolName = "NULL"; 168 VolName = "NULL";
162 Confine = false; 169 Confine = false;
163 } 170 }
>> 171
164 } 172 }
165 173
>> 174
166 void DMXParticleSource::SetAngDistType(G4Strin 175 void DMXParticleSource::SetAngDistType(G4String atype)
167 { 176 {
168 AngDistType = atype; 177 AngDistType = atype;
169 } 178 }
170 179
>> 180
171 void DMXParticleSource::GeneratePointSource() 181 void DMXParticleSource::GeneratePointSource()
172 { 182 {
173 // Generates Points given the point source. 183 // Generates Points given the point source.
174 if(SourcePosType == "Point") 184 if(SourcePosType == "Point")
175 particle_position = CentreCoords; 185 particle_position = CentreCoords;
176 else 186 else
177 if(verbosityLevel >= 1) 187 if(verbosityLevel >= 1)
178 G4cout << "Error SourcePosType is not se 188 G4cout << "Error SourcePosType is not set to Point" << G4endl;
179 } 189 }
180 190
>> 191
181 void DMXParticleSource::GeneratePointsInVolume 192 void DMXParticleSource::GeneratePointsInVolume()
182 { 193 {
183 G4ThreeVector RandPos; 194 G4ThreeVector RandPos;
184 G4double x=0., y=0., z=0.; 195 G4double x=0., y=0., z=0.;
185 196
186 if(SourcePosType != "Volume" && verbosityLev 197 if(SourcePosType != "Volume" && verbosityLevel >= 1)
187 G4cout << "Error SourcePosType not Volume" 198 G4cout << "Error SourcePosType not Volume" << G4endl;
188 199
189 if(Shape == "Sphere") { 200 if(Shape == "Sphere") {
190 x = Radius*2.; 201 x = Radius*2.;
191 y = Radius*2.; 202 y = Radius*2.;
192 z = Radius*2.; 203 z = Radius*2.;
193 while(((x*x)+(y*y)+(z*z)) > (Radius*Radius 204 while(((x*x)+(y*y)+(z*z)) > (Radius*Radius)) {
194 x = G4UniformRand(); 205 x = G4UniformRand();
195 y = G4UniformRand(); 206 y = G4UniformRand();
196 z = G4UniformRand(); 207 z = G4UniformRand();
197 208
198 x = (x*2.*Radius) - Radius; 209 x = (x*2.*Radius) - Radius;
199 y = (y*2.*Radius) - Radius; 210 y = (y*2.*Radius) - Radius;
200 z = (z*2.*Radius) - Radius; 211 z = (z*2.*Radius) - Radius;
201 } 212 }
202 } 213 }
203 214
204 else if(Shape == "Cylinder") { 215 else if(Shape == "Cylinder") {
205 x = Radius*2.; 216 x = Radius*2.;
206 y = Radius*2.; 217 y = Radius*2.;
207 while(((x*x)+(y*y)) > (Radius*Radius)) { 218 while(((x*x)+(y*y)) > (Radius*Radius)) {
208 x = G4UniformRand(); 219 x = G4UniformRand();
209 y = G4UniformRand(); 220 y = G4UniformRand();
210 z = G4UniformRand(); 221 z = G4UniformRand();
211 x = (x*2.*Radius) - Radius; 222 x = (x*2.*Radius) - Radius;
212 y = (y*2.*Radius) - Radius; 223 y = (y*2.*Radius) - Radius;
213 z = (z*2.*halfz) - halfz; 224 z = (z*2.*halfz) - halfz;
214 } 225 }
215 } 226 }
216 227
217 else 228 else
218 G4cout << "Error: Volume Shape Does Not Ex 229 G4cout << "Error: Volume Shape Does Not Exist" << G4endl;
219 230
220 RandPos.setX(x); 231 RandPos.setX(x);
221 RandPos.setY(y); 232 RandPos.setY(y);
222 RandPos.setZ(z); 233 RandPos.setZ(z);
223 particle_position = CentreCoords + RandPos; 234 particle_position = CentreCoords + RandPos;
224 235
225 } 236 }
226 237
>> 238
227 G4bool DMXParticleSource::IsSourceConfined() 239 G4bool DMXParticleSource::IsSourceConfined()
228 { 240 {
229 241
230 // Method to check point is within the volum 242 // Method to check point is within the volume specified
231 if(Confine == false) 243 if(Confine == false)
232 G4cout << "Error: Confine is false" << G4e 244 G4cout << "Error: Confine is false" << G4endl;
233 G4ThreeVector null_vec(0.,0.,0.); << 245 G4ThreeVector null(0.,0.,0.);
234 G4ThreeVector *ptr = &null_vec; << 246 G4ThreeVector *ptr;
>> 247 ptr = &null;
235 248
236 // Check particle_position is within VolName 249 // Check particle_position is within VolName
237 G4VPhysicalVolume *theVolume; 250 G4VPhysicalVolume *theVolume;
238 theVolume=gNavigator->LocateGlobalPointAndSe 251 theVolume=gNavigator->LocateGlobalPointAndSetup(particle_position,ptr,true);
239 G4String theVolName = theVolume->GetName(); 252 G4String theVolName = theVolume->GetName();
240 if(theVolName == VolName) { 253 if(theVolName == VolName) {
241 if(verbosityLevel >= 1) 254 if(verbosityLevel >= 1)
242 G4cout << "Particle is in volume " << Vo 255 G4cout << "Particle is in volume " << VolName << G4endl;
243 return(true); 256 return(true);
244 } 257 }
245 else 258 else
246 return(false); 259 return(false);
247 } 260 }
248 261
>> 262
249 void DMXParticleSource::SetParticleMomentumDir 263 void DMXParticleSource::SetParticleMomentumDirection
250 (G4ParticleMomentum aDirection) { 264 (G4ParticleMomentum aDirection) {
251 265
252 particle_momentum_direction = aDirection.un 266 particle_momentum_direction = aDirection.unit();
253 } 267 }
254 268
>> 269
255 void DMXParticleSource::GenerateIsotropicFlux( 270 void DMXParticleSource::GenerateIsotropicFlux()
256 { 271 {
257 272
258 G4double rndm, rndm2; 273 G4double rndm, rndm2;
259 G4double px, py, pz; 274 G4double px, py, pz;
260 275
261 G4double sintheta, sinphi, costheta, cosphi; 276 G4double sintheta, sinphi, costheta, cosphi;
262 rndm = G4UniformRand(); 277 rndm = G4UniformRand();
263 costheta = std::cos(MinTheta) - rndm * (std: << 278 costheta = std::cos(MinTheta) - rndm * (std::cos(MinTheta) - std::cos(MaxTheta));
264 - std::cos(Max <<
265 sintheta = std::sqrt(1. - costheta*costheta) 279 sintheta = std::sqrt(1. - costheta*costheta);
266 280
267 rndm2 = G4UniformRand(); 281 rndm2 = G4UniformRand();
268 Phi = MinPhi + (MaxPhi - MinPhi) * rndm2; 282 Phi = MinPhi + (MaxPhi - MinPhi) * rndm2;
269 sinphi = std::sin(Phi); 283 sinphi = std::sin(Phi);
270 cosphi = std::cos(Phi); 284 cosphi = std::cos(Phi);
271 285
272 px = -sintheta * cosphi; 286 px = -sintheta * cosphi;
273 py = -sintheta * sinphi; 287 py = -sintheta * sinphi;
274 pz = -costheta; 288 pz = -costheta;
275 289
276 G4double ResMag = std::sqrt((px*px) + (py*py 290 G4double ResMag = std::sqrt((px*px) + (py*py) + (pz*pz));
277 px = px/ResMag; 291 px = px/ResMag;
278 py = py/ResMag; 292 py = py/ResMag;
279 pz = pz/ResMag; 293 pz = pz/ResMag;
280 294
281 particle_momentum_direction.setX(px); 295 particle_momentum_direction.setX(px);
282 particle_momentum_direction.setY(py); 296 particle_momentum_direction.setY(py);
283 particle_momentum_direction.setZ(pz); 297 particle_momentum_direction.setZ(pz);
284 298
285 // particle_momentum_direction now holds uni 299 // particle_momentum_direction now holds unit momentum vector.
286 if(verbosityLevel >= 2) 300 if(verbosityLevel >= 2)
287 G4cout << "Generating isotropic vector: " << 301 G4cout << "Generating isotropic vector: " << particle_momentum_direction << G4endl;
288 << particle_momentum_direction << G <<
289 } 302 }
290 303
>> 304
291 void DMXParticleSource::SetEnergyDisType(G4Str 305 void DMXParticleSource::SetEnergyDisType(G4String DisType)
292 { 306 {
293 EnergyDisType = DisType; 307 EnergyDisType = DisType;
294 } 308 }
295 309
296 void DMXParticleSource::SetMonoEnergy(G4double 310 void DMXParticleSource::SetMonoEnergy(G4double menergy)
297 { 311 {
298 MonoEnergy = menergy; 312 MonoEnergy = menergy;
299 } 313 }
300 314
301 void DMXParticleSource::GenerateMonoEnergetic( 315 void DMXParticleSource::GenerateMonoEnergetic()
302 { 316 {
303 particle_energy = MonoEnergy; 317 particle_energy = MonoEnergy;
304 } 318 }
305 319
>> 320
306 void DMXParticleSource::SetVerbosity(int vL) 321 void DMXParticleSource::SetVerbosity(int vL)
307 { 322 {
308 verbosityLevel = vL; 323 verbosityLevel = vL;
309 G4cout << "Verbosity Set to: " << verbosityL 324 G4cout << "Verbosity Set to: " << verbosityLevel << G4endl;
310 } 325 }
311 326
312 void DMXParticleSource::SetParticleDefinition 327 void DMXParticleSource::SetParticleDefinition
313 (G4ParticleDefinition* aParticleDefinition) 328 (G4ParticleDefinition* aParticleDefinition)
314 { 329 {
315 particle_definition = aParticleDefinition; 330 particle_definition = aParticleDefinition;
316 particle_charge = particle_definition->GetPD 331 particle_charge = particle_definition->GetPDGCharge();
317 } 332 }
318 333
>> 334
319 void DMXParticleSource::GeneratePrimaryVertex( 335 void DMXParticleSource::GeneratePrimaryVertex(G4Event *evt)
320 { 336 {
321 337
322 if(particle_definition==nullptr) { << 338 if(particle_definition==NULL) {
323 G4cout << "No particle has been defined!" 339 G4cout << "No particle has been defined!" << G4endl;
324 return; 340 return;
325 } 341 }
326 342
327 // Position 343 // Position
328 G4bool srcconf = false; 344 G4bool srcconf = false;
329 G4int LoopCount = 0; 345 G4int LoopCount = 0;
330 346
331 while(srcconf == false) { 347 while(srcconf == false) {
332 if(SourcePosType == "Point") 348 if(SourcePosType == "Point")
333 GeneratePointSource(); 349 GeneratePointSource();
334 else if(SourcePosType == "Volume") 350 else if(SourcePosType == "Volume")
335 GeneratePointsInVolume(); 351 GeneratePointsInVolume();
336 else { 352 else {
337 G4cout << "Error: SourcePosType undefine 353 G4cout << "Error: SourcePosType undefined" << G4endl;
338 G4cout << "Generating point source" << G 354 G4cout << "Generating point source" << G4endl;
339 GeneratePointSource(); 355 GeneratePointSource();
340 } 356 }
341 if(Confine == true) { 357 if(Confine == true) {
342 srcconf = IsSourceConfined(); 358 srcconf = IsSourceConfined();
343 // if source in confined srcconf = true 359 // if source in confined srcconf = true terminating the loop
344 // if source isnt confined srcconf = fal 360 // if source isnt confined srcconf = false and loop continues
345 } 361 }
346 else if(Confine == false) 362 else if(Confine == false)
347 srcconf = true; // terminate loop 363 srcconf = true; // terminate loop
348 364
349 ++LoopCount; << 365 LoopCount++;
350 if(LoopCount == 100000) { 366 if(LoopCount == 100000) {
351 G4cout << "***************************** 367 G4cout << "*************************************" << G4endl;
352 G4cout << "LoopCount = 100000" << G4en 368 G4cout << "LoopCount = 100000" << G4endl;
353 G4cout << "Either the source distribut 369 G4cout << "Either the source distribution >> confinement" << G4endl;
354 G4cout << "or any confining volume may 370 G4cout << "or any confining volume may not overlap with" << G4endl;
355 G4cout << "the source distribution or 371 G4cout << "the source distribution or any confining volumes" << G4endl;
356 G4cout << "may not exist"<< G4endl; 372 G4cout << "may not exist"<< G4endl;
357 G4cout << "If you have set confine the 373 G4cout << "If you have set confine then this will be ignored" <<G4endl;
358 G4cout << "for this event." << G4endl; 374 G4cout << "for this event." << G4endl;
359 G4cout << "*************************** 375 G4cout << "*************************************" << G4endl;
360 srcconf = true; //Avoids an infinite l 376 srcconf = true; //Avoids an infinite loop
361 } 377 }
362 } 378 }
363 379
364 // Angular stuff 380 // Angular stuff
365 if(AngDistType == "iso") 381 if(AngDistType == "iso")
366 GenerateIsotropicFlux(); 382 GenerateIsotropicFlux();
367 else if(AngDistType == "direction") 383 else if(AngDistType == "direction")
368 SetParticleMomentumDirection(particle_mome 384 SetParticleMomentumDirection(particle_momentum_direction);
369 else 385 else
370 G4cout << "Error: AngDistType has unusual 386 G4cout << "Error: AngDistType has unusual value" << G4endl;
371 // Energy stuff 387 // Energy stuff
372 if(EnergyDisType == "Mono") 388 if(EnergyDisType == "Mono")
373 GenerateMonoEnergetic(); 389 GenerateMonoEnergetic();
374 else 390 else
375 G4cout << "Error: EnergyDisType has unusua 391 G4cout << "Error: EnergyDisType has unusual value" << G4endl;
376 392
377 // create a new vertex 393 // create a new vertex
378 G4PrimaryVertex* vertex = 394 G4PrimaryVertex* vertex =
379 new G4PrimaryVertex(particle_position,part 395 new G4PrimaryVertex(particle_position,particle_time);
380 396
381 if(verbosityLevel >= 2) 397 if(verbosityLevel >= 2)
382 G4cout << "Creating primaries and assignin 398 G4cout << "Creating primaries and assigning to vertex" << G4endl;
383 // create new primaries and set them to the 399 // create new primaries and set them to the vertex
384 G4double mass = particle_definition->GetPDG 400 G4double mass = particle_definition->GetPDGMass();
385 G4double energy = particle_energy + mass; 401 G4double energy = particle_energy + mass;
386 G4double pmom = std::sqrt(energy*energy-mass 402 G4double pmom = std::sqrt(energy*energy-mass*mass);
387 G4double px = pmom*particle_momentum_directi 403 G4double px = pmom*particle_momentum_direction.x();
388 G4double py = pmom*particle_momentum_directi 404 G4double py = pmom*particle_momentum_direction.y();
389 G4double pz = pmom*particle_momentum_directi 405 G4double pz = pmom*particle_momentum_direction.z();
390 406
391 if(verbosityLevel >= 1){ 407 if(verbosityLevel >= 1){
392 G4cout << "Particle name: " 408 G4cout << "Particle name: "
393 << particle_definition->GetParticle << 409 << particle_definition->GetParticleName() << G4endl;
394 G4cout << " Energy: "<<particle_ener 410 G4cout << " Energy: "<<particle_energy << G4endl;
395 G4cout << " Position: "<<particle_posi 411 G4cout << " Position: "<<particle_position<< G4endl;
396 G4cout << " Direction: "<<particle_mome 412 G4cout << " Direction: "<<particle_momentum_direction << G4endl;
397 G4cout << " NumberOfParticlesToBeGenerated 413 G4cout << " NumberOfParticlesToBeGenerated: "
398 << NumberOfParticlesToBeGenerated < << 414 << NumberOfParticlesToBeGenerated << G4endl;
399 } 415 }
400 416
401 for( G4int i=0; i<NumberOfParticlesToBeGener << 417
>> 418 for( G4int i=0; i<NumberOfParticlesToBeGenerated; i++ ) {
402 G4PrimaryParticle* particle = 419 G4PrimaryParticle* particle =
403 new G4PrimaryParticle(particle_definitio 420 new G4PrimaryParticle(particle_definition,px,py,pz);
404 particle->SetMass( mass ); 421 particle->SetMass( mass );
405 particle->SetCharge( particle_charge ); 422 particle->SetCharge( particle_charge );
406 particle->SetPolarization(particle_polariz 423 particle->SetPolarization(particle_polarization.x(),
407 particle_polariz << 424 particle_polarization.y(),
408 particle_polariz << 425 particle_polarization.z());
409 vertex->SetPrimary( particle ); 426 vertex->SetPrimary( particle );
410 } 427 }
411 evt->AddPrimaryVertex( vertex ); 428 evt->AddPrimaryVertex( vertex );
412 if(verbosityLevel > 1) 429 if(verbosityLevel > 1)
413 G4cout << " Primary Vetex generated "<< G4 430 G4cout << " Primary Vetex generated "<< G4endl;
414 } 431 }
>> 432
>> 433
>> 434
>> 435
415 436