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