| 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 |
|
|---|
