| 1 | // |
|---|
| 2 | // ******************************************************************** |
|---|
| 3 | // * License and Disclaimer * |
|---|
| 4 | // * * |
|---|
| 5 | // * The Geant4 software is copyright of the Copyright Holders of * |
|---|
| 6 | // * the Geant4 Collaboration. It is provided under the terms and * |
|---|
| 7 | // * conditions of the Geant4 Software License, included in the file * |
|---|
| 8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
|---|
| 9 | // * include a list of copyright holders. * |
|---|
| 10 | // * * |
|---|
| 11 | // * Neither the authors of this software system, nor their employing * |
|---|
| 12 | // * institutes,nor the agencies providing financial support for this * |
|---|
| 13 | // * work make any representation or warranty, express or implied, * |
|---|
| 14 | // * regarding this software system or assume any liability for its * |
|---|
| 15 | // * use. Please see the license in the file LICENSE and URL above * |
|---|
| 16 | // * for the full disclaimer and the limitation of liability. * |
|---|
| 17 | // * * |
|---|
| 18 | // * This code implementation is the result of the scientific and * |
|---|
| 19 | // * technical work of the GEANT4 collaboration. * |
|---|
| 20 | // * By using, copying, modifying or distributing the software (or * |
|---|
| 21 | // * any work based on the software) you agree to acknowledge its * |
|---|
| 22 | // * use in resulting scientific publications, and indicate your * |
|---|
| 23 | // * acceptance of all terms of the Geant4 Software license. * |
|---|
| 24 | // ******************************************************************** |
|---|
| 25 | // |
|---|
| 26 | |
|---|
| 27 | #include "ExamplePrimaryGeneratorAction.hh" |
|---|
| 28 | #include "G4Event.hh" |
|---|
| 29 | #include "G4ParticleTable.hh" |
|---|
| 30 | #include "G4ParticleDefinition.hh" |
|---|
| 31 | #include "globals.hh" |
|---|
| 32 | #include "G4ThreeVector.hh" |
|---|
| 33 | |
|---|
| 34 | #include "G4GeneralParticleSource.hh" |
|---|
| 35 | |
|---|
| 36 | ExamplePrimaryGeneratorAction::ExamplePrimaryGeneratorAction() |
|---|
| 37 | { |
|---|
| 38 | //DEBUG |
|---|
| 39 | DebugXmin = -999999.; |
|---|
| 40 | |
|---|
| 41 | particleGun = new G4GeneralParticleSource (); |
|---|
| 42 | } |
|---|
| 43 | |
|---|
| 44 | ExamplePrimaryGeneratorAction::~ExamplePrimaryGeneratorAction() |
|---|
| 45 | { |
|---|
| 46 | // if(verbosityLevel == 2) |
|---|
| 47 | //{ |
|---|
| 48 | // Always give the user debug info. |
|---|
| 49 | |
|---|
| 50 | G4cout << "Output of DEBUG stuff" << G4endl; |
|---|
| 51 | G4cout << "positional stuff" << G4endl; |
|---|
| 52 | G4cout << "Scale, X, Scale, Y, Scale, Z" << G4endl; |
|---|
| 53 | G4double scalex, scaley, scalez; |
|---|
| 54 | G4int i; |
|---|
| 55 | for( i=0; i<100; i++) |
|---|
| 56 | { |
|---|
| 57 | scalex = DebugXmin + (i+1)*DebugXStep; |
|---|
| 58 | scaley = DebugYmin + (i+1)*DebugYStep; |
|---|
| 59 | scalez = DebugZmin + (i+1)*DebugZStep; |
|---|
| 60 | G4cout << scalex << " " << debugx[i] << " " << scaley << " " << debugy[i] << " " << scalez << " " << debugz[i] << G4endl; |
|---|
| 61 | } |
|---|
| 62 | |
|---|
| 63 | G4cout << "Scale, Number Px, Py, Pz, Scale, Number Theta, Scale, Number Phi" << G4endl; |
|---|
| 64 | G4double scalep, scalet, scalephi; |
|---|
| 65 | for(i=0; i<100; i++) |
|---|
| 66 | { |
|---|
| 67 | scalep = -1 + (i+1)*0.02; |
|---|
| 68 | scalet = (i+1) * (pi/100.); |
|---|
| 69 | scalephi = (i+1) * (twopi/100.); |
|---|
| 70 | G4cout << scalep << " " << debugpx[i] << " " << debugpy[i] << " " << debugpz[i] << " " << scalet << " " << debugtheta[i] << " " << scalephi << " " << debugphi[i] << G4endl; |
|---|
| 71 | } |
|---|
| 72 | |
|---|
| 73 | G4cout << "Initial Energy, No. Of Events" << G4endl; |
|---|
| 74 | G4double ene_out = 0.; |
|---|
| 75 | for(i=0; i<100; i++) |
|---|
| 76 | { |
|---|
| 77 | if(EneDisType == "Mono") |
|---|
| 78 | ene_out = emin/2. + (i+1)*debug_energy_step; |
|---|
| 79 | // else if(EneDisType == "Arb") |
|---|
| 80 | //{ |
|---|
| 81 | // if(IntType == "Spline") |
|---|
| 82 | // ene_out = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(0)) + (i+1)*debug_energy_step; |
|---|
| 83 | // else |
|---|
| 84 | // ene_out = ArbEnergyH.GetLowEdgeEnergy(size_t(0)) + (i+1)*debug_energy_step; |
|---|
| 85 | //} |
|---|
| 86 | //else if(EnergyDisType == "Epn") |
|---|
| 87 | //{ |
|---|
| 88 | // ene_out = IPDFEnergyH.GetLowEdgeEnergy(size_t(0)) + (i+1)*debug_energy_step; |
|---|
| 89 | //} |
|---|
| 90 | else |
|---|
| 91 | ene_out = emin + (i+1)*debug_energy_step; |
|---|
| 92 | G4cout << ene_out << " " << debugenergy[i] << G4endl; |
|---|
| 93 | } |
|---|
| 94 | //} |
|---|
| 95 | // G4cout << "About to delete particleGun " << G4endl; |
|---|
| 96 | delete particleGun; |
|---|
| 97 | //G4cout << "Deleted particleGun " << G4endl; |
|---|
| 98 | } |
|---|
| 99 | |
|---|
| 100 | void ExamplePrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) |
|---|
| 101 | { |
|---|
| 102 | //G4cout << "About to generate primary vertex" << G4endl; |
|---|
| 103 | particleGun->GeneratePrimaryVertex(anEvent); |
|---|
| 104 | //G4cout << "Back " << DebugXmin << G4endl; |
|---|
| 105 | |
|---|
| 106 | if(DebugXmin == -999999.) |
|---|
| 107 | { |
|---|
| 108 | SourceType = particleGun->GetPosDisType(); |
|---|
| 109 | //G4cout << "SourceType " << SourceType << G4endl; |
|---|
| 110 | SourceShape = particleGun->GetPosDisShape(); |
|---|
| 111 | //G4cout << "SourceShape " << SourceShape << G4endl; |
|---|
| 112 | radius = particleGun->GetRadius(); |
|---|
| 113 | //G4cout << "radius " << radius << G4endl; |
|---|
| 114 | halfx = particleGun->GetHalfX(); |
|---|
| 115 | //G4cout << "halfx " << halfx << G4endl; |
|---|
| 116 | halfy = particleGun->GetHalfY(); |
|---|
| 117 | //G4cout << "halfy " << halfy << G4endl; |
|---|
| 118 | halfz = particleGun->GetHalfZ(); |
|---|
| 119 | //G4cout << "halfz " << halfz << G4endl; |
|---|
| 120 | centre = particleGun->GetCentreCoords(); |
|---|
| 121 | //G4cout << "centre " << centre << G4endl; |
|---|
| 122 | |
|---|
| 123 | // for use with energy |
|---|
| 124 | EneDisType = particleGun->GetEnergyDisType(); |
|---|
| 125 | //G4cout << "EneDisType " << EneDisType << G4endl; |
|---|
| 126 | InterpolationType = particleGun->GetIntType(); |
|---|
| 127 | //G4cout << "InterpolationType " << InterpolationType << G4endl; |
|---|
| 128 | if(EneDisType == "Arb") |
|---|
| 129 | { |
|---|
| 130 | emin = particleGun->GetArbEmin(); |
|---|
| 131 | //G4cout << "emin " << emin << G4endl; |
|---|
| 132 | emax = particleGun->GetArbEmax(); |
|---|
| 133 | //G4cout << "emax " << emax << G4endl; |
|---|
| 134 | } |
|---|
| 135 | else |
|---|
| 136 | { |
|---|
| 137 | emin = particleGun->GetEmin(); |
|---|
| 138 | //G4cout << "emin " << emin << G4endl; |
|---|
| 139 | emax = particleGun->GetEmax(); |
|---|
| 140 | //G4cout << "emax " << emax << G4endl; |
|---|
| 141 | } |
|---|
| 142 | } |
|---|
| 143 | |
|---|
| 144 | G4ThreeVector ParticlePos = particleGun->GetParticlePosition(); |
|---|
| 145 | //G4cout << "ParticlePos " << ParticlePos << G4endl; |
|---|
| 146 | G4ThreeVector ParticleMomDir = particleGun->GetParticleMomentumDirection(); |
|---|
| 147 | //G4cout << "ParticleMomDir " << ParticleMomDir << G4endl; |
|---|
| 148 | |
|---|
| 149 | //G4cout << "Starting DEBUG stuff " << DebugXmin << G4endl; |
|---|
| 150 | // DEBUG SECTION |
|---|
| 151 | // if(verbosityLevel == 2) |
|---|
| 152 | // G4cout << "Collecting DEBUG info" <<G4endl; |
|---|
| 153 | G4int idebug = 0; |
|---|
| 154 | // position |
|---|
| 155 | if(DebugXmin == -999999.) |
|---|
| 156 | { |
|---|
| 157 | //G4cout << "Here 11 " << SourceType << " " << centre <<G4endl; |
|---|
| 158 | if(SourceType == "Point") |
|---|
| 159 | { |
|---|
| 160 | // DEBUG - make Xmin etc 0.5 * point and Xmax etc 1.5 * point |
|---|
| 161 | if(centre.x() == 0.0) |
|---|
| 162 | { |
|---|
| 163 | DebugXmin = -2.; |
|---|
| 164 | DebugXmax = 2.; |
|---|
| 165 | } |
|---|
| 166 | else |
|---|
| 167 | { |
|---|
| 168 | DebugXmin = centre.x() * 0.5; |
|---|
| 169 | DebugXmax = centre.x() * 1.5; |
|---|
| 170 | } |
|---|
| 171 | |
|---|
| 172 | if(centre.y() == 0.0) |
|---|
| 173 | { |
|---|
| 174 | DebugYmin = -2.; |
|---|
| 175 | DebugYmax = 2.; |
|---|
| 176 | } |
|---|
| 177 | else |
|---|
| 178 | { |
|---|
| 179 | DebugYmin = centre.y() * 0.5; |
|---|
| 180 | DebugYmax = centre.y() * 1.5; |
|---|
| 181 | } |
|---|
| 182 | |
|---|
| 183 | if(centre.z() == 0.0) |
|---|
| 184 | { |
|---|
| 185 | DebugZmin = -2.; |
|---|
| 186 | DebugZmax = 2.; |
|---|
| 187 | } |
|---|
| 188 | else |
|---|
| 189 | { |
|---|
| 190 | DebugZmin = centre.z() * 0.5; |
|---|
| 191 | DebugZmax = centre.z() * 1.5; |
|---|
| 192 | } |
|---|
| 193 | } |
|---|
| 194 | else |
|---|
| 195 | { |
|---|
| 196 | //G4cout << "Here 11a " << SourceShape << G4endl; |
|---|
| 197 | if((SourceShape == "Circle") || (SourceShape == "Annulus") || (SourceShape == "Sphere")) |
|---|
| 198 | { |
|---|
| 199 | DebugZmax = radius; |
|---|
| 200 | } |
|---|
| 201 | else if((SourceShape == "Ellipse") || (SourceShape == "Ellipsoid")) |
|---|
| 202 | { |
|---|
| 203 | DebugZmax = halfx; |
|---|
| 204 | if(halfy > DebugZmax) |
|---|
| 205 | DebugZmax = halfy; |
|---|
| 206 | if(halfz > DebugZmax) |
|---|
| 207 | DebugZmax = halfz; |
|---|
| 208 | } |
|---|
| 209 | else if(SourceShape == "Square") |
|---|
| 210 | { |
|---|
| 211 | DebugZmax = halfx; |
|---|
| 212 | } |
|---|
| 213 | else if(SourceShape == "Rectangle") |
|---|
| 214 | { |
|---|
| 215 | DebugZmax = halfx; |
|---|
| 216 | if(DebugZmax < halfy) |
|---|
| 217 | DebugZmax = halfy; |
|---|
| 218 | } |
|---|
| 219 | else if(SourceShape == "Cylinder") |
|---|
| 220 | { |
|---|
| 221 | if(radius >= halfz) |
|---|
| 222 | DebugZmax = radius; |
|---|
| 223 | else |
|---|
| 224 | DebugZmax = halfz; |
|---|
| 225 | } |
|---|
| 226 | else if(SourceShape == "Para") |
|---|
| 227 | { |
|---|
| 228 | DebugZmax = halfx; |
|---|
| 229 | if(DebugZmax < halfy) |
|---|
| 230 | DebugZmax = halfy; |
|---|
| 231 | if(DebugZmax < halfz) |
|---|
| 232 | DebugZmax = halfz; |
|---|
| 233 | } |
|---|
| 234 | DebugZmax = 3 * DebugZmax; |
|---|
| 235 | DebugXmin = centre.x() - DebugZmax; |
|---|
| 236 | DebugYmin = centre.y() - DebugZmax; |
|---|
| 237 | DebugZmin = centre.z() - DebugZmax; |
|---|
| 238 | DebugXmax = centre.x() + DebugZmax; |
|---|
| 239 | DebugYmax = centre.y() + DebugZmax; |
|---|
| 240 | DebugZmax = centre.z() + DebugZmax; |
|---|
| 241 | } |
|---|
| 242 | DebugXStep = (DebugXmax - DebugXmin)/100.; |
|---|
| 243 | DebugYStep = (DebugYmax - DebugYmin)/100.; |
|---|
| 244 | DebugZStep = (DebugZmax - DebugZmin)/100.; |
|---|
| 245 | } |
|---|
| 246 | |
|---|
| 247 | //G4cout << "out of first bit" << G4endl; |
|---|
| 248 | idebug = 0; |
|---|
| 249 | G4double X_edge = DebugXmin; |
|---|
| 250 | //G4cout << "entering while loop 1 " << ParticlePos.x() << G4endl; |
|---|
| 251 | while (ParticlePos.x() > X_edge) |
|---|
| 252 | { |
|---|
| 253 | //G4cout << idebug << " " << ParticlePos.x() << " " << X_edge << G4endl; |
|---|
| 254 | X_edge = DebugXmin + (idebug+1)*DebugXStep; |
|---|
| 255 | idebug++; |
|---|
| 256 | } |
|---|
| 257 | debugx[idebug-1] = debugx[idebug-1] + 1; |
|---|
| 258 | idebug = 0; |
|---|
| 259 | G4double Y_edge = DebugYmin; |
|---|
| 260 | //G4cout << "entering while loop 2" << G4endl; |
|---|
| 261 | while (ParticlePos.y() > Y_edge) |
|---|
| 262 | { |
|---|
| 263 | Y_edge = DebugYmin + (idebug+1)*DebugYStep; |
|---|
| 264 | idebug++; |
|---|
| 265 | } |
|---|
| 266 | debugy[idebug-1] = debugy[idebug-1] + 1; |
|---|
| 267 | idebug = 0; |
|---|
| 268 | G4double Z_edge = DebugZmin; |
|---|
| 269 | //G4cout << "entering while loop 3" << G4endl; |
|---|
| 270 | while (ParticlePos.z() > Z_edge) |
|---|
| 271 | { |
|---|
| 272 | Z_edge = DebugZmin + (idebug+1)*DebugZStep; |
|---|
| 273 | idebug++; |
|---|
| 274 | } |
|---|
| 275 | debugz[idebug-1] = debugz[idebug-1] + 1; |
|---|
| 276 | |
|---|
| 277 | //G4cout << "Here 12" << G4endl; |
|---|
| 278 | // trajectory |
|---|
| 279 | // px, py, pz are unit vectors so arrays run -1 to 1 |
|---|
| 280 | idebug = 0; |
|---|
| 281 | G4double Px_edge = -1.; |
|---|
| 282 | //G4cout << "Loop 1" << G4endl; |
|---|
| 283 | while (ParticleMomDir.x() > Px_edge) |
|---|
| 284 | { |
|---|
| 285 | Px_edge = -1 + (idebug+1)*0.02; |
|---|
| 286 | idebug++; |
|---|
| 287 | } |
|---|
| 288 | debugpx[idebug-1] = debugpx[idebug-1] + 1; |
|---|
| 289 | idebug = 0; |
|---|
| 290 | G4double Py_edge = -1.; |
|---|
| 291 | //G4cout << "Loop 2" << G4endl; |
|---|
| 292 | while (ParticleMomDir.y() > Py_edge) |
|---|
| 293 | { |
|---|
| 294 | Py_edge = -1 + (idebug+1)*0.02; |
|---|
| 295 | idebug++; |
|---|
| 296 | } |
|---|
| 297 | debugpy[idebug-1] = debugpy[idebug-1] + 1; |
|---|
| 298 | idebug = 0; |
|---|
| 299 | G4double Pz_edge = -1.; |
|---|
| 300 | //G4cout << "Loop 3" << G4endl; |
|---|
| 301 | while (ParticleMomDir.z() > Pz_edge) |
|---|
| 302 | { |
|---|
| 303 | Pz_edge = -1 + (idebug+1)*0.02; |
|---|
| 304 | idebug++; |
|---|
| 305 | } |
|---|
| 306 | debugpz[idebug-1] = debugpz[idebug-1] + 1; |
|---|
| 307 | |
|---|
| 308 | G4double theta = particleGun->GetTheta(); |
|---|
| 309 | G4double phi = particleGun->GetPhi(); |
|---|
| 310 | |
|---|
| 311 | //G4cout << "Here 13" << G4endl; |
|---|
| 312 | // Theta ranges 0-Pi, and Phi goes 0-two pi. |
|---|
| 313 | idebug = 0; |
|---|
| 314 | G4double Theta_edge = 0; |
|---|
| 315 | while (theta > Theta_edge) |
|---|
| 316 | { |
|---|
| 317 | Theta_edge = (idebug+1) * (pi/100.); |
|---|
| 318 | idebug++; |
|---|
| 319 | } |
|---|
| 320 | debugtheta[idebug-1] = debugtheta[idebug-1] + 1; |
|---|
| 321 | idebug = 0; |
|---|
| 322 | G4double Phi_edge = 0.; |
|---|
| 323 | while (phi > Phi_edge) |
|---|
| 324 | { |
|---|
| 325 | Phi_edge = (idebug+1) * (twopi/100.); |
|---|
| 326 | idebug++; |
|---|
| 327 | } |
|---|
| 328 | debugphi[idebug-1] = debugphi[idebug-1] + 1; |
|---|
| 329 | |
|---|
| 330 | //G4cout << "Here 14" << G4endl; |
|---|
| 331 | // Energy |
|---|
| 332 | if(EneDisType == "Mono") |
|---|
| 333 | debug_energy_step = emin/100.; |
|---|
| 334 | // else if(EneDisType == "Arb") |
|---|
| 335 | // { |
|---|
| 336 | // if(InterpolationType == "Spline") |
|---|
| 337 | //{ |
|---|
| 338 | // G4int len = G4int(IPDFArbEnergyH.GetVectorLength()); |
|---|
| 339 | // debug_energy_step = (IPDFArbEnergyH.GetLowEdgeEnergy(size_t(len-1)) - IPDFArbEnergyH.GetLowEdgeEnergy(size_t(0)))/100.; |
|---|
| 340 | //} |
|---|
| 341 | // else |
|---|
| 342 | //{ |
|---|
| 343 | // G4int len = G4int(ArbEnergyH.GetVectorLength()); |
|---|
| 344 | // debug_energy_step = (ArbEnergyH.GetLowEdgeEnergy(size_t(len-1)) - ArbEnergyH.GetLowEdgeEnergy(size_t(0)))/100.; |
|---|
| 345 | //} |
|---|
| 346 | //} |
|---|
| 347 | //else if(EneDisType == "Epn") |
|---|
| 348 | // { |
|---|
| 349 | // G4int len = G4int(IPDFEnergyH.GetVectorLength()); |
|---|
| 350 | // debug_energy_step = (IPDFEnergyH.GetLowEdgeEnergy(size_t(len-1)) - IPDFEnergyH.GetLowEdgeEnergy(size_t(0)))/100.; |
|---|
| 351 | //} |
|---|
| 352 | else |
|---|
| 353 | debug_energy_step = (emax - emin)/100.; |
|---|
| 354 | |
|---|
| 355 | //G4cout << "Here 15" << G4endl; |
|---|
| 356 | G4double PartEnergy = particleGun->GetParticleEnergy(); |
|---|
| 357 | //G4cout << "Energy is " << PartEnergy << " " << emin << " " << emax << G4endl; |
|---|
| 358 | |
|---|
| 359 | G4double Ebin_edge = 0.; |
|---|
| 360 | idebug = 0; |
|---|
| 361 | while (PartEnergy > Ebin_edge) |
|---|
| 362 | { |
|---|
| 363 | if(EneDisType == "Mono") |
|---|
| 364 | Ebin_edge = emin/2. + (idebug+1)*debug_energy_step; |
|---|
| 365 | // else if(EneDisType == "Arb") |
|---|
| 366 | //{ |
|---|
| 367 | // if(InterpolationType == "Spline") |
|---|
| 368 | // Ebin_edge = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(0)) + (idebug+1)*debug_energy_step; |
|---|
| 369 | // else |
|---|
| 370 | // Ebin_edge = ArbEnergyH.GetLowEdgeEnergy(size_t(0)) + (idebug+1)*debug_energy_step; |
|---|
| 371 | //} |
|---|
| 372 | // else if(EneDisType == "Epn") |
|---|
| 373 | //{ |
|---|
| 374 | // Ebin_edge = IPDFEnergyH.GetLowEdgeEnergy(size_t(0)) + (idebug+1)*debug_energy_step; |
|---|
| 375 | //} |
|---|
| 376 | else |
|---|
| 377 | Ebin_edge = emin + (idebug+1)*debug_energy_step; |
|---|
| 378 | idebug++; |
|---|
| 379 | } |
|---|
| 380 | debugenergy[idebug-1] = debugenergy[idebug-1] + 1; |
|---|
| 381 | |
|---|
| 382 | // G4cout << "debug-energy_step " << debug_energy_step << " " << Ebin_edge << G4endl; |
|---|
| 383 | //G4cout << "Ending thingy" << G4endl; |
|---|
| 384 | |
|---|
| 385 | } |
|---|
| 386 | |
|---|
| 387 | |
|---|
| 388 | |
|---|
