xray_fluorescence/src/XrayFluoPrimaryGeneratorAction.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: Elena Guardincerri (Elena.Guardincerri@ge.infn.it)
0029 //
0030 // History:
0031 // -----------
0032 // 28 Nov 2001 Elena Guardincerri     Created
0033 //
0034 // -------------------------------------------------------------------
0035
0036 #include "XrayFluoPrimaryGeneratorAction.hh"
0037 #include "XrayFluoDetectorConstruction.hh"
0038 #include "XrayFluoPrimaryGeneratorMessenger.hh"
0039 #include "XrayFluoRunAction.hh"
0040 #include "G4Event.hh"
0041 #include "G4Gamma.hh"
0042 #include "G4ParticleGun.hh"
0043 #include "G4ParticleTable.hh"
0044 #include "G4ParticleDefinition.hh"
0045 #include "G4MTRunManager.hh"
0046 #include "Randomize.hh"
0047 #include "XrayFluoAnalysisManager.hh"
0048 #include "XrayFluoDataSet.hh"
0049 #include "G4PhysicalConstants.hh"
0050 #include "G4SystemOfUnits.hh"
0051
0052 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0053
0054 XrayFluoPrimaryGeneratorAction::XrayFluoPrimaryGeneratorAction(const
0055                                    XrayFluoDetectorConstruction* XrayFluoDC)
0056   :rndmFlag("off"),beam("off"),spectrum("off"),isoVert("off"),phaseSpaceGunFlag(false),
0057    rayleighFlag(true), detectorPosition(0)
0058 {
0059   runAction = 0;
0060   XrayFluoDetector = XrayFluoDC;
0061
0062   G4int n_particle = 1;
0063   particleGun  = new G4ParticleGun(n_particle);
0064
0065   //create a messenger for this class
0066   gunMessenger = new XrayFluoPrimaryGeneratorMessenger(this);
0067
0068   // default particle kinematic
0069   G4ParticleDefinition* particle
0070     = G4Gamma::Definition();
0071   particleGun->SetParticleDefinition(particle);
0072   particleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,1.));
0073   particleGun->SetParticleEnergy(10. * keV);
0074
0075   G4double position = -0.5*(XrayFluoDetector->GetWorldSizeZ());
0076   particleGun->SetParticlePosition(G4ThreeVector(0.*cm,0.*cm,position));
0077
0078   G4cout << "XrayFluoPrimaryGeneratorAction created" << G4endl;
0079
0080 }
0081
0082 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0083
0084 void XrayFluoPrimaryGeneratorAction::ActivatePhaseSpace(G4String fileName) {
0085
0086   // load phase-space
0087   phaseSpaceGunFlag = true;
0088
0089   // reads the data stored on disk form previous runs
0090   // and get these data to data members
0091
0092   XrayFluoAnalysisManager* analysis =  XrayFluoAnalysisManager::getInstance();
0093   analysis->LoadGunData(fileName, rayleighFlag);
0094   detectorPosition = XrayFluoDetector->GetDetectorPosition();
0095   detectorPosition.setR(detectorPosition.r()-(5.*cm)); // 5 cm before the detector, so in front of it.
0096
0097 }
0098
0099 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0100
0101 void XrayFluoPrimaryGeneratorAction::SetRayleighFlag (G4bool value)
0102 {
0103   rayleighFlag = value;
0104 }
0105
0106
0107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0108
0109 XrayFluoPrimaryGeneratorAction::~XrayFluoPrimaryGeneratorAction()
0110 {
0111   delete particleGun;
0112   delete gunMessenger;
0113 }
0114
0115 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0116
0117 void XrayFluoPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
0118 {
0119   //retrieve runAction, if not done
0120   if (!runAction)
0121     {
0122       //Sequential runaction
0123       if (G4RunManager::GetRunManager()->GetRunManagerType() ==
0124       G4RunManager::sequentialRM)
0125     runAction = static_cast<const XrayFluoRunAction*>
0126       (G4RunManager::GetRunManager()->GetUserRunAction());
0127       else //MT master runaction
0128     runAction = static_cast<const XrayFluoRunAction*>
0129       (G4MTRunManager::GetMasterRunManager()->GetUserRunAction());
0130       if (!runAction)
0131     G4cout << "Something wrong here!" << G4endl;
0132     }
0133
0134   //this function is called at the begining of event
0135   //
0136   G4double z0 = -0.5*(XrayFluoDetector->GetWorldSizeZ());
0137   G4double y0 = 0.*cm, x0 = 0.*cm;
0138   if (rndmFlag == "on")
0139     {
0140       y0 = (XrayFluoDetector->GetDia3SizeXY())/std::sqrt(2.)*(G4UniformRand()-0.5); // it was GetSampleSizeXY(),
0141       x0 = (XrayFluoDetector->GetDia3SizeXY())/std::sqrt(2.)*(G4UniformRand()-0.5); // not divided by std::sqrt(2.)
0142     }
0143   particleGun->SetParticlePosition(G4ThreeVector(x0,y0,z0));
0144
0145   //randomize starting point
0146   if (beam == "on")
0147     {
0148       G4double radius = 0.5 * mm;
0149       G4double rho = radius*std::sqrt(G4UniformRand());
0150       G4double theta = 2*pi*G4UniformRand()*rad;
0151       G4double position = -0.5*(XrayFluoDetector->GetWorldSizeZ());
0152
0153       G4double y = rho * std::sin(theta);
0154       G4double x = rho * std::cos(theta);
0155
0156       particleGun->SetParticlePosition(G4ThreeVector(x,y,position));
0157     }
0158   //shoot particles according to a certain spectrum
0159   if (spectrum =="on")
0160     {
0161       G4String particle =  particleGun->GetParticleDefinition()
0162     ->GetParticleName();
0163       if(particle == "proton"|| particle == "alpha")
0164     {
0165       G4DataVector* energies =  runAction->GetEnergies();
0166       G4DataVector* data =  runAction->GetData();
0167
0168       G4double sum = runAction->GetDataSum();
0169       G4double partSum = 0;
0170       G4int j = 0;
0171       G4double random= sum*G4UniformRand();
0172       while (partSum<random)
0173         {
0174           partSum += (*data)[j];
0175           j++;
0176         }
0177
0178       particleGun->SetParticleEnergy((*energies)[j]);
0179
0180     }
0181       else if (particle == "gamma")
0182     {
0183       const XrayFluoDataSet* dataSet = runAction->GetGammaSet();
0184
0185       G4int i = 0;
0186       G4int id = 0;
0187       G4double minEnergy = 0. * keV;
0188       G4double particleEnergy= 0.;
0189       G4double maxEnergy = 10. * keV;
0190       G4double energyRange = maxEnergy - minEnergy;
0191
0192        while ( i == 0)
0193         {
0194           G4double random = G4UniformRand();
0195
0196           G4double randomNum = G4UniformRand(); //*5.0E6;
0197
0198           particleEnergy = (random*energyRange) + minEnergy;
0199
0200           if ((dataSet->FindValue(particleEnergy,id)) > randomNum)
0201         {
0202           i = 1;
0203
0204         }
0205         }
0206        particleGun->SetParticleEnergy(particleEnergy);
0207     }
0208     }
0209
0210   // Randomize starting point and direction
0211
0212   if (isoVert == "on")
0213     {
0214       G4double rho = 1. *m;
0215       //theta in [0;pi/2]
0216       G4double theta = (pi/2)*G4UniformRand();
0217       //phi in [-pi;pi]
0218       G4double phi = (G4UniformRand()*2*pi)- pi;
0219       G4double x = rho*std::sin(theta)*std::sin(phi);
0220       G4double y = rho*std::sin(theta)*std::cos(phi);
0221       G4double z = -(rho*std::cos(theta));
0222       particleGun->SetParticlePosition(G4ThreeVector(x,y,z));
0223
0224       G4double Xdim = XrayFluoDetector->GetSampleSizeXY();
0225       G4double Ydim = XrayFluoDetector->GetSampleSizeXY();
0226
0227       G4double Dx = Xdim*(G4UniformRand()-0.5);
0228
0229       G4double Dy = Ydim*(G4UniformRand()-0.5);
0230
0231       particleGun->SetParticleMomentumDirection(G4ThreeVector(-x+Dx,-y+Dy,-z));
0232
0233     }
0234
0235   // using prevoiously genereated emissions from sample.....
0236
0237   if (phaseSpaceGunFlag){
0238
0239     particleGun->SetParticlePosition(detectorPosition);
0240     particleGun->SetParticleMomentumDirection(detectorPosition);
0241
0242     G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
0243
0244     const std::pair<G4double,G4String> kine =
0245       XrayFluoAnalysisManager::getInstance()->GetEmittedParticleEnergyAndType();
0246
0247     G4double energy = kine.first;
0248     G4ParticleDefinition* particle = particleTable->FindParticle(kine.second);
0249
0250     particleGun->SetParticleEnergy(energy);
0251     particleGun->SetParticleDefinition(particle);
0252
0253
0254   }
0255
0256   G4double partEnergy = particleGun->GetParticleEnergy();
0257   XrayFluoAnalysisManager* analysis =  XrayFluoAnalysisManager::getInstance();
0258   analysis->analysePrimaryGenerator(partEnergy/keV);
0259
0260
0261   particleGun->GeneratePrimaryVertex(anEvent);
0262 }
0263
0264 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0265
0266
0267
0268
0269
0270
0271
0272
0273