underground_physics/src/DMXParticleSource.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 // --------------------------------------------------------------
0028 //   GEANT 4 - Underground Dark Matter Detector Advanced Example
0029 //
0030 //      For information related to this code contact: Alex Howard
0031 //      e-mail: alexander.howard@cern.ch
0032 // --------------------------------------------------------------
0033 // Comments
0034 //
0035 //                  Underground Advanced
0036 //               by A. Howard and H. Araujo
0037 //                    (27th November 2001)
0038 //
0039 //
0040 // ParticleSource program
0041 // --------------------------------------------------------------
0042 //////////////////////////////////////////////////////////////////////////////
0043 // This particle source is a shortened version of G4GeneralParticleSource by
0044 // C Ferguson, F Lei & P Truscott (University of Southampton / DERA), with
0045 // some minor modifications.
0046 //////////////////////////////////////////////////////////////////////////////
0047
0048 #include <cmath>
0049
0050 #include "DMXParticleSource.hh"
0051
0052 #include "G4PhysicalConstants.hh"
0053 #include "G4SystemOfUnits.hh"
0054 #include "G4PrimaryParticle.hh"
0055 #include "G4Event.hh"
0056 #include "Randomize.hh"
0057 #include "G4TransportationManager.hh"
0058 #include "G4VPhysicalVolume.hh"
0059 #include "G4PhysicalVolumeStore.hh"
0060 #include "G4ParticleTable.hh"
0061 #include "G4ParticleDefinition.hh"
0062 #include "G4IonTable.hh"
0063 #include "G4Ions.hh"
0064 #include "G4TrackingManager.hh"
0065 #include "G4Track.hh"
0066
0067
0068 DMXParticleSource::DMXParticleSource() {
0069
0070   NumberOfParticlesToBeGenerated = 1;
0071   particle_definition = nullptr;
0072   G4ThreeVector zero(0., 0., 0.);
0073   particle_momentum_direction = G4ParticleMomentum(1., 0., 0.);
0074   particle_energy = 1.0*MeV;
0075   particle_position = zero;
0076   particle_time = 0.0;
0077   particle_polarization = zero;
0078   particle_charge = 0.0;
0079
0080   SourcePosType = "Volume";
0081   Shape = "NULL";
0082   halfz = 0.;
0083   Radius = 0.;
0084   CentreCoords = zero;
0085   Confine = false;
0086   VolName = "NULL";
0087
0088   AngDistType = "iso";
0089   MinTheta = 0.;
0090   MaxTheta = pi;
0091   MinPhi = 0.;
0092   MaxPhi = twopi;
0093
0094   EnergyDisType = "Mono";
0095   MonoEnergy = 1*MeV;
0096
0097   verbosityLevel = 0;
0098
0099   theMessenger = new DMXParticleSourceMessenger(this);
0100   gNavigator = G4TransportationManager::GetTransportationManager()
0101     ->GetNavigatorForTracking();
0102 }
0103
0104 DMXParticleSource::~DMXParticleSource()
0105 {
0106   delete theMessenger;
0107 }
0108
0109 void DMXParticleSource::SetPosDisType(G4String PosType)
0110 {
0111   SourcePosType = PosType;
0112 }
0113
0114 void DMXParticleSource::SetPosDisShape(G4String shapeType)
0115 {
0116   Shape = shapeType;
0117 }
0118
0119 void DMXParticleSource::SetCentreCoords(G4ThreeVector coordsOfCentre)
0120 {
0121   CentreCoords = coordsOfCentre;
0122 }
0123
0124 void DMXParticleSource::SetHalfZ(G4double zhalf)
0125 {
0126   halfz = zhalf;
0127 }
0128
0129 void DMXParticleSource::SetRadius(G4double radius)
0130 {
0131   Radius = radius;
0132 }
0133
0134 void DMXParticleSource::ConfineSourceToVolume(G4String Vname)
0135 {
0136   VolName = Vname;
0137   if(verbosityLevel == 2) G4cout << VolName << G4endl;
0138
0139   // checks if selected volume exists
0140   G4VPhysicalVolume *tempPV      = nullptr;
0141   G4PhysicalVolumeStore *PVStore = nullptr;
0142   G4String theRequiredVolumeName = VolName;
0143   PVStore = G4PhysicalVolumeStore::GetInstance();
0144   G4bool found = false;
0145   if(verbosityLevel == 2) G4cout << PVStore->size() << G4endl;
0146
0147   tempPV = PVStore->GetVolume(theRequiredVolumeName, false);
0148   if (tempPV != nullptr) { found = true; }
0149
0150   // found = true then the volume exists else it doesnt.
0151   if(found == true) {
0152     if(verbosityLevel >= 1)
0153       G4cout << "Volume " << VolName << " exists" << G4endl;
0154     Confine = true;
0155   }
0156   else if(VolName=="NULL")
0157     Confine = false;
0158   else {
0159     G4cout << " **** Error: Volume does not exist **** " << G4endl;
0160     G4cout << " Ignoring confine condition" << G4endl;
0161     VolName = "NULL";
0162     Confine = false;
0163   }
0164 }
0165
0166 void DMXParticleSource::SetAngDistType(G4String atype)
0167 {
0168   AngDistType = atype;
0169 }
0170
0171 void DMXParticleSource::GeneratePointSource()
0172 {
0173   // Generates Points given the point source.
0174   if(SourcePosType == "Point")
0175     particle_position = CentreCoords;
0176   else
0177     if(verbosityLevel >= 1)
0178       G4cout << "Error SourcePosType is not set to Point" << G4endl;
0179 }
0180
0181 void DMXParticleSource::GeneratePointsInVolume()
0182 {
0183   G4ThreeVector RandPos;
0184   G4double x=0., y=0., z=0.;
0185
0186   if(SourcePosType != "Volume" && verbosityLevel >= 1)
0187     G4cout << "Error SourcePosType not Volume" << G4endl;
0188
0189   if(Shape == "Sphere") {
0190     x = Radius*2.;
0191     y = Radius*2.;
0192     z = Radius*2.;
0193     while(((x*x)+(y*y)+(z*z)) > (Radius*Radius)) {
0194       x = G4UniformRand();
0195       y = G4UniformRand();
0196       z = G4UniformRand();
0197
0198       x = (x*2.*Radius) - Radius;
0199       y = (y*2.*Radius) - Radius;
0200       z = (z*2.*Radius) - Radius;
0201     }
0202   }
0203
0204   else if(Shape == "Cylinder") {
0205     x = Radius*2.;
0206     y = Radius*2.;
0207     while(((x*x)+(y*y)) > (Radius*Radius)) {
0208       x = G4UniformRand();
0209       y = G4UniformRand();
0210       z = G4UniformRand();
0211       x = (x*2.*Radius) - Radius;
0212       y = (y*2.*Radius) - Radius;
0213       z = (z*2.*halfz) - halfz;
0214     }
0215   }
0216
0217   else
0218     G4cout << "Error: Volume Shape Does Not Exist" << G4endl;
0219
0220   RandPos.setX(x);
0221   RandPos.setY(y);
0222   RandPos.setZ(z);
0223   particle_position = CentreCoords + RandPos;
0224
0225 }
0226
0227 G4bool DMXParticleSource::IsSourceConfined()
0228 {
0229
0230   // Method to check point is within the volume specified
0231   if(Confine == false)
0232     G4cout << "Error: Confine is false" << G4endl;
0233   G4ThreeVector null_vec(0.,0.,0.);
0234   G4ThreeVector *ptr = &null_vec;
0235
0236   // Check particle_position is within VolName
0237   G4VPhysicalVolume *theVolume;
0238   theVolume=gNavigator->LocateGlobalPointAndSetup(particle_position,ptr,true);
0239   G4String theVolName = theVolume->GetName();
0240   if(theVolName == VolName) {
0241     if(verbosityLevel >= 1)
0242       G4cout << "Particle is in volume " << VolName << G4endl;
0243     return(true);
0244   }
0245   else
0246     return(false);
0247 }
0248
0249 void DMXParticleSource::SetParticleMomentumDirection
0250    (G4ParticleMomentum aDirection) {
0251
0252   particle_momentum_direction =  aDirection.unit();
0253 }
0254
0255 void DMXParticleSource::GenerateIsotropicFlux()
0256 {
0257
0258   G4double rndm, rndm2;
0259   G4double px, py, pz;
0260
0261   G4double sintheta, sinphi, costheta, cosphi;
0262   rndm = G4UniformRand();
0263   costheta = std::cos(MinTheta) - rndm * (std::cos(MinTheta)
0264                                 - std::cos(MaxTheta));
0265   sintheta = std::sqrt(1. - costheta*costheta);
0266
0267   rndm2 = G4UniformRand();
0268   Phi = MinPhi + (MaxPhi - MinPhi) * rndm2;
0269   sinphi = std::sin(Phi);
0270   cosphi = std::cos(Phi);
0271
0272   px = -sintheta * cosphi;
0273   py = -sintheta * sinphi;
0274   pz = -costheta;
0275
0276   G4double ResMag = std::sqrt((px*px) + (py*py) + (pz*pz));
0277   px = px/ResMag;
0278   py = py/ResMag;
0279   pz = pz/ResMag;
0280
0281   particle_momentum_direction.setX(px);
0282   particle_momentum_direction.setY(py);
0283   particle_momentum_direction.setZ(pz);
0284
0285   // particle_momentum_direction now holds unit momentum vector.
0286   if(verbosityLevel >= 2)
0287     G4cout << "Generating isotropic vector: "
0288            << particle_momentum_direction << G4endl;
0289 }
0290
0291 void DMXParticleSource::SetEnergyDisType(G4String DisType)
0292 {
0293   EnergyDisType = DisType;
0294 }
0295
0296 void DMXParticleSource::SetMonoEnergy(G4double menergy)
0297 {
0298   MonoEnergy = menergy;
0299 }
0300
0301 void DMXParticleSource::GenerateMonoEnergetic()
0302 {
0303   particle_energy = MonoEnergy;
0304 }
0305
0306 void DMXParticleSource::SetVerbosity(int vL)
0307 {
0308   verbosityLevel = vL;
0309   G4cout << "Verbosity Set to: " << verbosityLevel << G4endl;
0310 }
0311
0312 void DMXParticleSource::SetParticleDefinition
0313   (G4ParticleDefinition* aParticleDefinition)
0314 {
0315   particle_definition = aParticleDefinition;
0316   particle_charge = particle_definition->GetPDGCharge();
0317 }
0318
0319 void DMXParticleSource::GeneratePrimaryVertex(G4Event *evt)
0320 {
0321
0322   if(particle_definition==nullptr) {
0323     G4cout << "No particle has been defined!" << G4endl;
0324     return;
0325   }
0326
0327   // Position
0328   G4bool srcconf = false;
0329   G4int LoopCount = 0;
0330
0331   while(srcconf == false)  {
0332     if(SourcePosType == "Point")
0333       GeneratePointSource();
0334     else if(SourcePosType == "Volume")
0335       GeneratePointsInVolume();
0336     else {
0337       G4cout << "Error: SourcePosType undefined" << G4endl;
0338       G4cout << "Generating point source" << G4endl;
0339       GeneratePointSource();
0340     }
0341     if(Confine == true) {
0342       srcconf = IsSourceConfined();
0343       // if source in confined srcconf = true terminating the loop
0344       // if source isnt confined srcconf = false and loop continues
0345     }
0346     else if(Confine == false)
0347       srcconf = true; // terminate loop
0348
0349     ++LoopCount;
0350     if(LoopCount == 100000) {
0351       G4cout << "*************************************" << G4endl;
0352         G4cout << "LoopCount = 100000" << G4endl;
0353         G4cout << "Either the source distribution >> confinement" << G4endl;
0354         G4cout << "or any confining volume may not overlap with" << G4endl;
0355         G4cout << "the source distribution or any confining volumes" << G4endl;
0356         G4cout << "may not exist"<< G4endl;
0357         G4cout << "If you have set confine then this will be ignored" <<G4endl;
0358         G4cout << "for this event." << G4endl;
0359         G4cout << "*************************************" << G4endl;
0360         srcconf = true; //Avoids an infinite loop
0361       }
0362   }
0363
0364   // Angular stuff
0365   if(AngDistType == "iso")
0366     GenerateIsotropicFlux();
0367   else if(AngDistType == "direction")
0368     SetParticleMomentumDirection(particle_momentum_direction);
0369   else
0370     G4cout << "Error: AngDistType has unusual value" << G4endl;
0371   // Energy stuff
0372   if(EnergyDisType == "Mono")
0373     GenerateMonoEnergetic();
0374   else
0375     G4cout << "Error: EnergyDisType has unusual value" << G4endl;
0376
0377   // create a new vertex
0378   G4PrimaryVertex* vertex =
0379     new G4PrimaryVertex(particle_position,particle_time);
0380
0381   if(verbosityLevel >= 2)
0382     G4cout << "Creating primaries and assigning to vertex" << G4endl;
0383   // create new primaries and set them to the vertex
0384   G4double mass =  particle_definition->GetPDGMass();
0385   G4double energy = particle_energy + mass;
0386   G4double pmom = std::sqrt(energy*energy-mass*mass);
0387   G4double px = pmom*particle_momentum_direction.x();
0388   G4double py = pmom*particle_momentum_direction.y();
0389   G4double pz = pmom*particle_momentum_direction.z();
0390
0391   if(verbosityLevel >= 1){
0392     G4cout << "Particle name: "
0393            << particle_definition->GetParticleName() << G4endl;
0394     G4cout << "       Energy: "<<particle_energy << G4endl;
0395     G4cout << "     Position: "<<particle_position<< G4endl;
0396     G4cout << "    Direction: "<<particle_momentum_direction << G4endl;
0397     G4cout << " NumberOfParticlesToBeGenerated: "
0398            << NumberOfParticlesToBeGenerated << G4endl;
0399   }
0400
0401   for( G4int i=0; i<NumberOfParticlesToBeGenerated; ++i ) {
0402     G4PrimaryParticle* particle =
0403       new G4PrimaryParticle(particle_definition,px,py,pz);
0404     particle->SetMass( mass );
0405     particle->SetCharge( particle_charge );
0406     particle->SetPolarization(particle_polarization.x(),
0407                               particle_polarization.y(),
0408                               particle_polarization.z());
0409     vertex->SetPrimary( particle );
0410   }
0411   evt->AddPrimaryVertex( vertex );
0412   if(verbosityLevel > 1)
0413     G4cout << " Primary Vetex generated "<< G4endl;
0414 }