xray_fluorescence/src/XrayFluoMercuryPrimaryGeneratorAction.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 "XrayFluoMercuryPrimaryGeneratorAction.hh"
0037 #include "XrayFluoMercuryDetectorConstruction.hh"
0038 #include "XrayFluoMercuryPrimaryGeneratorMessenger.hh"
0039 #include "XrayFluoRunAction.hh"
0040 #include "XrayFluoAnalysisManager.hh"
0041 #include "XrayFluoDataSet.hh"
0042 #include "G4PhysicalConstants.hh"
0043 #include "G4SystemOfUnits.hh"
0044 #include "G4DataVector.hh"
0045 #include "G4Event.hh"
0046 #include "G4ParticleGun.hh"
0047 #include "G4ParticleTable.hh"
0048 #include "G4ParticleDefinition.hh"
0049 #include "Randomize.hh"
0050
0051 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0052
0053 XrayFluoMercuryPrimaryGeneratorAction::XrayFluoMercuryPrimaryGeneratorAction(const XrayFluoMercuryDetectorConstruction* XrayFluoDC)
0054   :globalFlag(false),spectrum("off")
0055 {
0056
0057   XrayFluoDetector = XrayFluoDC;
0058
0059   G4int n_particle = 1;
0060   particleGun  = new G4ParticleGun(n_particle);
0061
0062   //create a messenger for this class
0063   gunMessenger = new XrayFluoMercuryPrimaryGeneratorMessenger(this);
0064   runManager = new XrayFluoRunAction();
0065
0066   // default particle kinematic
0067
0068   G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
0069   G4String particleName;
0070   G4ParticleDefinition* particle
0071     = particleTable->FindParticle(particleName="gamma");
0072   particleGun->SetParticleDefinition(particle);
0073   particleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,-1.));
0074
0075
0076   particleGun->SetParticleEnergy(10.*keV);
0077   G4double position = -0.5*(XrayFluoDetector->GetWorldSizeZ());
0078   particleGun->SetParticlePosition(G4ThreeVector(0.*cm,0.*cm,position));
0079
0080   G4cout << "XrayFluoMercuryPrimaryGeneratorAction created" << G4endl;
0081
0082 }
0083
0084
0085 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0086
0087 XrayFluoMercuryPrimaryGeneratorAction::~XrayFluoMercuryPrimaryGeneratorAction()
0088 {
0089   delete particleGun;
0090   delete gunMessenger;
0091   delete runManager;
0092
0093   G4cout << "XrayFluoMercuryPrimaryGeneratorAction deleted" << G4endl;
0094
0095 }
0096
0097 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0098
0099 void XrayFluoMercuryPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
0100 {
0101   //this function is called at the begining of event
0102   //
0103
0104   // Conidering the sunas a Poin-like source.
0105
0106   G4double z0 = -0.5*(XrayFluoDetector->GetWorldSizeZ());
0107   G4double y0 = 0.*m, x0 = 0.*m;
0108
0109
0110   // Let's try to illuminate only the prtion of Mercury surface that can be seen by the detector.
0111
0112   G4double spacecraftLatitude = XrayFluoDetector->GetOrbitInclination();
0113   G4double mercuryDia = XrayFluoDetector->GetMercuryDia();
0114   G4double sunDia = XrayFluoDetector->GetSunDia();
0115   G4double opticField = XrayFluoDetector->GetOpticAperture();
0116
0117
0118   G4double a = 2*std::tan(opticField/2);
0119
0120   //  if (!pointLikeFlag) {
0121
0122     // let's decide from wich point of the sun surface the particle is coming:
0123
0124   G4double theta = std::acos(2.*G4UniformRand() - 1.0);
0125     G4double phi = 2. * pi * G4UniformRand();
0126     G4double rho = sunDia/2;
0127
0128     G4double sunPosX = x0 + rho * std::sin(theta) * std::cos(phi);
0129     G4double sunPosY = y0 + rho * std::sin(theta) * std::sin(phi);
0130     G4double sunPosZ = z0 + rho * std::cos(theta);
0131
0132     particleGun->SetParticlePosition(G4ThreeVector(sunPosX,sunPosY,sunPosZ));
0133
0134     // the angle at the center of Mercury subtending the area seen by the optics:
0135     G4double alpha = 2 * a/mercuryDia;
0136
0137     if(!globalFlag){
0138       theta = alpha * G4UniformRand() + (180.*deg - spacecraftLatitude)-alpha/2.;
0139       phi = alpha * G4UniformRand() + 90. * deg - alpha/2.;
0140     }
0141
0142     else if(globalFlag){
0143       theta = pi/2. * rad * G4UniformRand() + 90.*deg ; //was 900., probably an error
0144       phi = 2*pi*rad * G4UniformRand() ;
0145     }
0146
0147     rho = mercuryDia/2.;
0148
0149     G4double mercuryPosX = rho * std::sin(theta) * std::cos(phi);
0150     G4double mercuryPosY = rho * std::sin(theta) * std::sin(phi);
0151     G4double mercuryPosZ = rho * std::cos(theta);
0152
0153     particleGun->SetParticleMomentumDirection(
0154                 G4ThreeVector(mercuryPosX-sunPosX ,mercuryPosY-sunPosY,mercuryPosZ-sunPosZ));
0155
0156     //  }
0157 //   if (pointLikeFlag) {
0158
0159 //   // theta is the angle that the mean direction of the incident light (on the desired
0160 //   // point of the surface of Mercury) makes with  the Z-axis
0161 //   G4double theta = std::asin( mercuryDia/2. * std::sin(spacecraftLatitude) /
0162 //           std::sqrt(std::pow(z0,2)+std::pow(mercuryDia/2.,2)-2*mercuryDia/2.*z0*std::cos(spacecraftLatitude)) );
0163
0164 //   // on the y axis, the light emitted from the Sun must be in [theta-phi;theta+phi]
0165 //   G4double phi = std::asin( mercuryDia/2.*std::sin(spacecraftLatitude) + a*std::cos(spacecraftLatitude) /
0166 //             std::sqrt( std::pow(mercuryDia/2.*std::sin(spacecraftLatitude) + a*std::cos(spacecraftLatitude) , 2) +
0167 //               std::pow(z0 - mercuryDia/2.*std::cos(spacecraftLatitude) - a*std::sin(spacecraftLatitude) , 2)) )
0168 //     - theta;
0169
0170 //   // on the x axis, the light emitted from the Sun must be in [-zeta;zeta]
0171 //   G4double zeta = std::atan( a/std::sqrt(std::pow(z0,2)+std::pow(mercuryDia,2)-2*mercuryDia*z0*std::cos(spacecraftLatitude)) );
0172
0173
0174
0175 //   //alpha in [-zeta;zeta]
0176 //   G4double alpha = (2*zeta)*G4UniformRand() - zeta;
0177 //   //beta in [theta-phi;theta+phi]
0178 //   G4double beta = (G4UniformRand()*2*phi) - phi + theta;
0179
0180 //   G4double dirY = std::sin(beta);
0181 //   G4double dirX = std::sin(alpha);
0182
0183 //   particleGun->SetParticleMomentumDirection(G4ThreeVector(dirX.,dirY,1.));
0184
0185 //   particleGun->SetParticlePosition(G4ThreeVector(x0,y0,z0));
0186
0187 //   }
0188
0189
0190
0191   //shoot particles according to a certain spectrum
0192   if (spectrum =="on")
0193     {
0194       G4String particle =  particleGun->GetParticleDefinition()
0195     ->GetParticleName();
0196       if(particle == "proton"|| particle == "alpha")
0197     {
0198       G4DataVector* energies =  runManager->GetEnergies();
0199       G4DataVector* data =  runManager->GetData();
0200
0201       G4double sum = runManager->GetDataSum();
0202       G4double partSum = 0;
0203       G4int j = 0;
0204       G4double random= sum*G4UniformRand();
0205       while (partSum<random)
0206         {
0207           partSum += (*data)[j];
0208           j++;
0209         }
0210
0211       particleGun->SetParticleEnergy((*energies)[j]);
0212
0213     }
0214       else if (particle == "gamma")
0215     {
0216       const XrayFluoDataSet* dataSet = runManager->GetGammaSet();
0217
0218       G4int i = 0;
0219       G4int id = 0;
0220       G4double minEnergy = 0. * keV;
0221       G4double particleEnergy= 0.;
0222       G4double maxEnergy = 10. * keV;
0223       G4double energyRange = maxEnergy - minEnergy;
0224
0225        while ( i == 0)
0226         {
0227           G4double random = G4UniformRand();
0228
0229           G4double randomNum = G4UniformRand(); //*5.0E6;
0230
0231           particleEnergy = (random*energyRange) + minEnergy;
0232
0233           if ((dataSet->FindValue(particleEnergy,id)) > randomNum)
0234         {
0235           i = 1;
0236
0237         }
0238         }
0239        particleGun->SetParticleEnergy(particleEnergy);
0240     }
0241     }
0242
0243
0244 #ifdef G4ANALYSIS_USE
0245
0246   G4double partEnergy = particleGun->GetParticleEnergy();
0247   XrayFluoAnalysisManager* analysis =  XrayFluoAnalysisManager::getInstance();
0248   analysis->analysePrimaryGenerator(partEnergy/keV);
0249
0250 #endif
0251
0252   particleGun->GeneratePrimaryVertex(anEvent);
0253 }
0254
0255 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0256
0257
0258
0259
0260
0261
0262
0263
0264