xray_fluorescence/src/XrayFluoPlanePrimaryGeneratorAction.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 // Author: Alfonso Mantero (Alfonso.Mantero@ge.infn.it)
0029 //
0030 // History:
0031 // -----------
0032 // 02 Sep 2003 Alfonso Mantero created
0033 //
0034 // -------------------------------------------------------------------
0035
0036 #include "XrayFluoPlanePrimaryGeneratorAction.hh"
0037 #include "XrayFluoPlaneDetectorConstruction.hh"
0038 #include "XrayFluoPlanePrimaryGeneratorMessenger.hh"
0039 #include "XrayFluoRunAction.hh"
0040 #include "XrayFluoAnalysisManager.hh"
0041 #include "XrayFluoDataSet.hh"
0042 #include "G4PhysicalConstants.hh"
0043 #include "G4SystemOfUnits.hh"
0044 #include "G4UnitsTable.hh"
0045 #include "G4DataVector.hh"
0046 #include "G4Event.hh"
0047 #include "G4ParticleGun.hh"
0048 #include "G4ParticleTable.hh"
0049 #include "G4ParticleDefinition.hh"
0050 #include "Randomize.hh"
0051
0052 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0053
0054 XrayFluoPlanePrimaryGeneratorAction::XrayFluoPlanePrimaryGeneratorAction(const XrayFluoPlaneDetectorConstruction* XrayFluoDC)
0055   :rndmFlag("on"),beam("off"),spectrum("off"),isoVert("off")
0056 {
0057
0058   XrayFluoDetector = XrayFluoDC;
0059
0060   G4int n_particle = 1;
0061   particleGun  = new G4ParticleGun(n_particle);
0062
0063   //create a messenger for this class
0064   gunMessenger = new XrayFluoPlanePrimaryGeneratorMessenger(this);
0065   runManager = new XrayFluoRunAction();
0066
0067   // default particle kinematic
0068
0069   G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
0070   G4String particleName;
0071   G4ParticleDefinition* particle
0072     = particleTable->FindParticle(particleName="gamma");
0073   particleGun->SetParticleDefinition(particle);
0074   particleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,1.));
0075
0076
0077   particleGun->SetParticleEnergy(10.*keV);
0078   G4double position = -0.5*(XrayFluoDetector->GetWorldSizeZ());
0079   particleGun->SetParticlePosition(G4ThreeVector(0.*cm,0.*cm,position));
0080
0081   G4cout << "XrayFluoPlanePrimaryGeneratorAction created  UUUUUUUUUUAAAAAAAAAAAAAAAAAAAAAAAaa" << G4endl;
0082
0083 }
0084
0085
0086 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0087
0088 XrayFluoPlanePrimaryGeneratorAction::~XrayFluoPlanePrimaryGeneratorAction()
0089 {
0090   delete particleGun;
0091   delete gunMessenger;
0092   delete runManager;
0093
0094   G4cout << "XrayFluoPlanePrimaryGeneratorAction deleted" << G4endl;
0095
0096 }
0097
0098 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0099
0100 void XrayFluoPlanePrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
0101 {
0102   //this function is called at the begining of event
0103   //
0104   G4double z0 = -0.5*(XrayFluoDetector->GetWorldSizeZ());
0105   G4double y0 = 0.*m, x0 = 0.*m;
0106   G4double dX = 0.5*(XrayFluoDetector->GetWorldSizeXY())-(XrayFluoDetector->GetPlaneSizeXY());
0107   if (rndmFlag == "on")
0108
0109     {y0 = (XrayFluoDetector->GetPlaneSizeXY())*(G4UniformRand()-0.5);
0110     x0 = (XrayFluoDetector->GetPlaneSizeXY())*(G4UniformRand()-0.5) + dX;
0111     }
0112
0113   z0 = -1 * dX;
0114
0115   particleGun->SetParticleMomentumDirection(G4ThreeVector(-1.,0.,1.));
0116
0117   particleGun->SetParticlePosition(G4ThreeVector(x0,y0,z0));
0118
0119   //randomize starting point
0120   if (beam == "on")
0121     {
0122       G4double radius = 0.5 * mm;
0123       G4double rho = radius*std::sqrt(G4UniformRand());
0124       G4double theta = 2*pi*G4UniformRand()*rad;
0125       G4double position = -0.5*(XrayFluoDetector->GetWorldSizeZ());
0126
0127       G4double y = rho * std::sin(theta);
0128       G4double x = rho * std::cos(theta);
0129
0130       particleGun->SetParticlePosition(G4ThreeVector(x,y,position));
0131     }
0132   //shoot particles according to a certain spectrum
0133   if (spectrum =="on")
0134     {
0135       G4String particle =  particleGun->GetParticleDefinition()
0136     ->GetParticleName();
0137       if(particle == "proton"|| particle == "alpha")
0138     {
0139       G4DataVector* energies =  runManager->GetEnergies();
0140       G4DataVector* data =  runManager->GetData();
0141
0142       G4double sum = runManager->GetDataSum();
0143       G4double partSum = 0;
0144       G4int j = 0;
0145       G4double random= sum*G4UniformRand();
0146       while (partSum<random)
0147         {
0148           partSum += (*data)[j];
0149           j++;
0150         }
0151
0152       particleGun->SetParticleEnergy((*energies)[j]);
0153
0154     }
0155       else if (particle == "gamma")
0156     {
0157       const XrayFluoDataSet* dataSet = runManager->GetGammaSet();
0158
0159       G4int i = 0;
0160       G4int id = 0;
0161       G4double minEnergy = 0. * keV;
0162       G4double particleEnergy= 0.;
0163       G4double maxEnergy = 10. * keV;
0164       G4double energyRange = maxEnergy - minEnergy;
0165
0166        while ( i == 0)
0167         {
0168           G4double random = G4UniformRand();
0169
0170           G4double randomNum = G4UniformRand(); //*5.0E6;
0171
0172           particleEnergy = (random*energyRange) + minEnergy;
0173
0174           if ((dataSet->FindValue(particleEnergy,id)) > randomNum)
0175         {
0176           i = 1;
0177
0178         }
0179         }
0180        particleGun->SetParticleEnergy(particleEnergy);
0181     }
0182     }
0183
0184   if (isoVert == "on")
0185     {
0186       G4double rho = 1. *m;
0187       //theta in [0;pi/2]
0188       G4double theta = (pi/2)*G4UniformRand();
0189       //phi in [-pi;pi]
0190       G4double phi = (G4UniformRand()*2*pi)- pi;
0191       G4double x = rho*std::sin(theta)*std::sin(phi);
0192       G4double y = rho*std::sin(theta)*std::cos(phi);
0193       G4double z = -(rho*std::cos(theta));
0194       particleGun->SetParticlePosition(G4ThreeVector(x,y,z));
0195
0196       G4double Xdim = XrayFluoDetector->GetPlaneSizeXY();
0197       G4double Ydim = XrayFluoDetector->GetPlaneSizeXY();
0198
0199       G4double Dx = Xdim*(G4UniformRand()-0.5);
0200
0201       G4double Dy = Ydim*(G4UniformRand()-0.5);
0202
0203       particleGun->SetParticleMomentumDirection(G4ThreeVector(-x+Dx,-y+Dy,-z));
0204
0205     }
0206 #ifdef G4ANALYSIS_USE
0207
0208   G4double partEnergy = particleGun->GetParticleEnergy();
0209   XrayFluoAnalysisManager* analysis =  XrayFluoAnalysisManager::getInstance();
0210   analysis->analysePrimaryGenerator(partEnergy/keV);
0211
0212 #endif
0213
0214   particleGun->GeneratePrimaryVertex(anEvent);
0215 }
0216
0217 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0218
0219
0220
0221
0222
0223
0224
0225
0226