Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/xray_fluorescence/src/XrayFluoMercuryPrimaryGeneratorAction.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /examples/advanced/xray_fluorescence/src/XrayFluoMercuryPrimaryGeneratorAction.cc (Version 11.3.0) and /examples/advanced/xray_fluorescence/src/XrayFluoMercuryPrimaryGeneratorAction.cc (Version 10.0.p2)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 //                                                 26 //
                                                   >>  27 // $Id: XrayFluoPlanePrimaryGeneratorAction.cc
                                                   >>  28 // GEANT4 tag $Name: 
 27 //                                                 29 //
 28 // Author: Alfonso Mantero (Alfonso.Mantero@ge     30 // Author: Alfonso Mantero (Alfonso.Mantero@ge.infn.it)
 29 //                                                 31 //
 30 // History:                                        32 // History:
 31 // -----------                                     33 // -----------
 32 // 02 Sep 2003 Alfonso Mantero created             34 // 02 Sep 2003 Alfonso Mantero created
 33 //                                                 35 //
 34 // -------------------------------------------     36 // -------------------------------------------------------------------
 35                                                    37 
 36 #include "XrayFluoMercuryPrimaryGeneratorActio     38 #include "XrayFluoMercuryPrimaryGeneratorAction.hh"
 37 #include "XrayFluoMercuryDetectorConstruction.     39 #include "XrayFluoMercuryDetectorConstruction.hh"
 38 #include "XrayFluoMercuryPrimaryGeneratorMesse     40 #include "XrayFluoMercuryPrimaryGeneratorMessenger.hh"
 39 #include "XrayFluoRunAction.hh"                    41 #include "XrayFluoRunAction.hh"
 40 #include "XrayFluoAnalysisManager.hh"              42 #include "XrayFluoAnalysisManager.hh"
 41 #include "XrayFluoDataSet.hh"                      43 #include "XrayFluoDataSet.hh"
 42 #include "G4PhysicalConstants.hh"                  44 #include "G4PhysicalConstants.hh"
 43 #include "G4SystemOfUnits.hh"                      45 #include "G4SystemOfUnits.hh"
 44 #include "G4DataVector.hh"                         46 #include "G4DataVector.hh"
 45 #include "G4Event.hh"                              47 #include "G4Event.hh"
 46 #include "G4ParticleGun.hh"                        48 #include "G4ParticleGun.hh"
 47 #include "G4ParticleTable.hh"                      49 #include "G4ParticleTable.hh"
 48 #include "G4ParticleDefinition.hh"                 50 #include "G4ParticleDefinition.hh"
 49 #include "Randomize.hh"                            51 #include "Randomize.hh"
 50                                                    52 
 51 //....oooOO0OOooo........oooOO0OOooo........oo     53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 52                                                    54 
 53 XrayFluoMercuryPrimaryGeneratorAction::XrayFlu <<  55 XrayFluoMercuryPrimaryGeneratorAction::XrayFluoMercuryPrimaryGeneratorAction(XrayFluoMercuryDetectorConstruction* XrayFluoDC)
 54   :globalFlag(false),spectrum("off")               56   :globalFlag(false),spectrum("off")
 55 {                                                  57 {
 56                                                    58 
 57   XrayFluoDetector = XrayFluoDC;                   59   XrayFluoDetector = XrayFluoDC;
 58                                                    60 
 59   G4int n_particle = 1;                            61   G4int n_particle = 1;
 60   particleGun  = new G4ParticleGun(n_particle)     62   particleGun  = new G4ParticleGun(n_particle);
 61                                                    63   
 62   //create a messenger for this class              64   //create a messenger for this class
 63   gunMessenger = new XrayFluoMercuryPrimaryGen     65   gunMessenger = new XrayFluoMercuryPrimaryGeneratorMessenger(this);
 64   runManager = new XrayFluoRunAction();            66   runManager = new XrayFluoRunAction();
 65                                                    67   
 66   // default particle kinematic                    68   // default particle kinematic
 67                                                    69   
 68   G4ParticleTable* particleTable = G4ParticleT     70   G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
 69   G4String particleName;                           71   G4String particleName;
 70   G4ParticleDefinition* particle                   72   G4ParticleDefinition* particle
 71     = particleTable->FindParticle(particleName     73     = particleTable->FindParticle(particleName="gamma");
 72   particleGun->SetParticleDefinition(particle)     74   particleGun->SetParticleDefinition(particle);
 73   particleGun->SetParticleMomentumDirection(G4     75   particleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,-1.));
 74                                                    76   
 75                                                    77 
 76   particleGun->SetParticleEnergy(10.*keV);         78   particleGun->SetParticleEnergy(10.*keV);
 77   G4double position = -0.5*(XrayFluoDetector->     79   G4double position = -0.5*(XrayFluoDetector->GetWorldSizeZ());
 78   particleGun->SetParticlePosition(G4ThreeVect     80   particleGun->SetParticlePosition(G4ThreeVector(0.*cm,0.*cm,position));
 79                                                    81 
 80   G4cout << "XrayFluoMercuryPrimaryGeneratorAc     82   G4cout << "XrayFluoMercuryPrimaryGeneratorAction created" << G4endl;
 81                                                    83   
 82 }                                                  84 }
 83                                                    85 
 84                                                    86 
 85 //....oooOO0OOooo........oooOO0OOooo........oo     87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 86                                                    88 
 87 XrayFluoMercuryPrimaryGeneratorAction::~XrayFl     89 XrayFluoMercuryPrimaryGeneratorAction::~XrayFluoMercuryPrimaryGeneratorAction()
 88 {                                                  90 {
 89   delete particleGun;                              91   delete particleGun;
 90   delete gunMessenger;                             92   delete gunMessenger;
 91   delete runManager;                               93   delete runManager;
 92                                                    94 
 93   G4cout << "XrayFluoMercuryPrimaryGeneratorAc     95   G4cout << "XrayFluoMercuryPrimaryGeneratorAction deleted" << G4endl;
 94                                                    96 
 95 }                                                  97 }
 96                                                    98 
 97 //....oooOO0OOooo........oooOO0OOooo........oo     99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 98                                                   100 
 99 void XrayFluoMercuryPrimaryGeneratorAction::Ge    101 void XrayFluoMercuryPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
100 {                                                 102 {
101   //this function is called at the begining of    103   //this function is called at the begining of event
102   //                                              104   // 
103                                                   105 
104   // Conidering the sunas a Poin-like source.     106   // Conidering the sunas a Poin-like source.
105                                                   107 
106   G4double z0 = -0.5*(XrayFluoDetector->GetWor    108   G4double z0 = -0.5*(XrayFluoDetector->GetWorldSizeZ());
107   G4double y0 = 0.*m, x0 = 0.*m;                  109   G4double y0 = 0.*m, x0 = 0.*m;
108                                                   110 
109                                                   111 
110   // Let's try to illuminate only the prtion o    112   // Let's try to illuminate only the prtion of Mercury surface that can be seen by the detector.
111                                                   113 
112   G4double spacecraftLatitude = XrayFluoDetect    114   G4double spacecraftLatitude = XrayFluoDetector->GetOrbitInclination();
113   G4double mercuryDia = XrayFluoDetector->GetM    115   G4double mercuryDia = XrayFluoDetector->GetMercuryDia();
114   G4double sunDia = XrayFluoDetector->GetSunDi    116   G4double sunDia = XrayFluoDetector->GetSunDia();
115   G4double opticField = XrayFluoDetector->GetO    117   G4double opticField = XrayFluoDetector->GetOpticAperture();
116                                                   118   
117                                                   119  
118   G4double a = 2*std::tan(opticField/2);          120   G4double a = 2*std::tan(opticField/2);
119                                                   121   
120   //  if (!pointLikeFlag) {                       122   //  if (!pointLikeFlag) {
121                                                   123 
122     // let's decide from wich point of the sun    124     // let's decide from wich point of the sun surface the particle is coming:
123                                                   125 
124   G4double theta = std::acos(2.*G4UniformRand(    126   G4double theta = std::acos(2.*G4UniformRand() - 1.0);
125     G4double phi = 2. * pi * G4UniformRand();     127     G4double phi = 2. * pi * G4UniformRand();
126     G4double rho = sunDia/2;                      128     G4double rho = sunDia/2;                         
127                                                   129 
128     G4double sunPosX = x0 + rho * std::sin(the    130     G4double sunPosX = x0 + rho * std::sin(theta) * std::cos(phi);
129     G4double sunPosY = y0 + rho * std::sin(the    131     G4double sunPosY = y0 + rho * std::sin(theta) * std::sin(phi);
130     G4double sunPosZ = z0 + rho * std::cos(the    132     G4double sunPosZ = z0 + rho * std::cos(theta);           
131                                                   133 
132     particleGun->SetParticlePosition(G4ThreeVe    134     particleGun->SetParticlePosition(G4ThreeVector(sunPosX,sunPosY,sunPosZ));
133                                                   135 
134     // the angle at the center of Mercury subt    136     // the angle at the center of Mercury subtending the area seen by the optics:
135     G4double alpha = 2 * a/mercuryDia;            137     G4double alpha = 2 * a/mercuryDia; 
136                                                   138     
137     if(!globalFlag){                              139     if(!globalFlag){
138       theta = alpha * G4UniformRand() + (180.*    140       theta = alpha * G4UniformRand() + (180.*deg - spacecraftLatitude)-alpha/2.;
139       phi = alpha * G4UniformRand() + 90. * de    141       phi = alpha * G4UniformRand() + 90. * deg - alpha/2.;     
140     }                                             142     }    
141                                                   143     
142     else if(globalFlag){                          144     else if(globalFlag){
143       theta = pi/2. * rad * G4UniformRand() +     145       theta = pi/2. * rad * G4UniformRand() + 90.*deg ; //was 900., probably an error
144       phi = 2*pi*rad * G4UniformRand() ;          146       phi = 2*pi*rad * G4UniformRand() ;
145     }                                             147     }
146                                                   148     
147     rho = mercuryDia/2.;                          149     rho = mercuryDia/2.;                        
148                                                   150     
149     G4double mercuryPosX = rho * std::sin(thet    151     G4double mercuryPosX = rho * std::sin(theta) * std::cos(phi);
150     G4double mercuryPosY = rho * std::sin(thet    152     G4double mercuryPosY = rho * std::sin(theta) * std::sin(phi);
151     G4double mercuryPosZ = rho * std::cos(thet    153     G4double mercuryPosZ = rho * std::cos(theta);           
152                                                   154     
153     particleGun->SetParticleMomentumDirection(    155     particleGun->SetParticleMomentumDirection(
154           G4ThreeVector(mercuryPosX-sunPosX ,m    156           G4ThreeVector(mercuryPosX-sunPosX ,mercuryPosY-sunPosY,mercuryPosZ-sunPosZ));
155                                                   157 
156     //  }                                         158     //  }
157 //   if (pointLikeFlag) {                         159 //   if (pointLikeFlag) {
158                                                   160 
159 //   // theta is the angle that the mean direc    161 //   // theta is the angle that the mean direction of the incident light (on the desired 
160 //   // point of the surface of Mercury) makes    162 //   // point of the surface of Mercury) makes with  the Z-axis
161 //   G4double theta = std::asin( mercuryDia/2.    163 //   G4double theta = std::asin( mercuryDia/2. * std::sin(spacecraftLatitude) / 
162 //       std::sqrt(std::pow(z0,2)+std::pow(mer    164 //       std::sqrt(std::pow(z0,2)+std::pow(mercuryDia/2.,2)-2*mercuryDia/2.*z0*std::cos(spacecraftLatitude)) );  
163                                                   165 
164 //   // on the y axis, the light emitted from     166 //   // on the y axis, the light emitted from the Sun must be in [theta-phi;theta+phi]
165 //   G4double phi = std::asin( mercuryDia/2.*s    167 //   G4double phi = std::asin( mercuryDia/2.*std::sin(spacecraftLatitude) + a*std::cos(spacecraftLatitude) /
166 //           std::sqrt( std::pow(mercuryDia/2.    168 //           std::sqrt( std::pow(mercuryDia/2.*std::sin(spacecraftLatitude) + a*std::cos(spacecraftLatitude) , 2) +
167 //           std::pow(z0 - mercuryDia/2.*std::    169 //           std::pow(z0 - mercuryDia/2.*std::cos(spacecraftLatitude) - a*std::sin(spacecraftLatitude) , 2)) ) 
168 //     - theta;                                   170 //     - theta;  
169                                                   171   
170 //   // on the x axis, the light emitted from     172 //   // on the x axis, the light emitted from the Sun must be in [-zeta;zeta]
171 //   G4double zeta = std::atan( a/std::sqrt(st    173 //   G4double zeta = std::atan( a/std::sqrt(std::pow(z0,2)+std::pow(mercuryDia,2)-2*mercuryDia*z0*std::cos(spacecraftLatitude)) );
172                                                   174   
173                                                   175   
174                                                   176   
175 //   //alpha in [-zeta;zeta]                      177 //   //alpha in [-zeta;zeta]
176 //   G4double alpha = (2*zeta)*G4UniformRand()    178 //   G4double alpha = (2*zeta)*G4UniformRand() - zeta;
177 //   //beta in [theta-phi;theta+phi]              179 //   //beta in [theta-phi;theta+phi]
178 //   G4double beta = (G4UniformRand()*2*phi) -    180 //   G4double beta = (G4UniformRand()*2*phi) - phi + theta;
179                                                   181   
180 //   G4double dirY = std::sin(beta);              182 //   G4double dirY = std::sin(beta);
181 //   G4double dirX = std::sin(alpha);             183 //   G4double dirX = std::sin(alpha);
182                                                   184   
183 //   particleGun->SetParticleMomentumDirection    185 //   particleGun->SetParticleMomentumDirection(G4ThreeVector(dirX.,dirY,1.));
184                                                   186   
185 //   particleGun->SetParticlePosition(G4ThreeV    187 //   particleGun->SetParticlePosition(G4ThreeVector(x0,y0,z0));
186                                                   188 
187 //   }                                            189 //   }
188                                                   190 
189                                                   191 
190                                                   192   
191   //shoot particles according to a certain spe    193   //shoot particles according to a certain spectrum
192   if (spectrum =="on")                            194   if (spectrum =="on")
193     {                                             195     {
194       G4String particle =  particleGun->GetPar    196       G4String particle =  particleGun->GetParticleDefinition()
195   ->GetParticleName();                            197   ->GetParticleName();
196       if(particle == "proton"|| particle == "a    198       if(particle == "proton"|| particle == "alpha")
197   {                                               199   {
198     G4DataVector* energies =  runManager->GetE    200     G4DataVector* energies =  runManager->GetEnergies();
199     G4DataVector* data =  runManager->GetData(    201     G4DataVector* data =  runManager->GetData();
200                                                   202    
201     G4double sum = runManager->GetDataSum();      203     G4double sum = runManager->GetDataSum();
202     G4double partSum = 0;                         204     G4double partSum = 0;
203     G4int j = 0;                                  205     G4int j = 0;
204     G4double random= sum*G4UniformRand();         206     G4double random= sum*G4UniformRand();
205     while (partSum<random)                        207     while (partSum<random)
206       {                                           208       {
207         partSum += (*data)[j];                    209         partSum += (*data)[j];
208         j++;                                      210         j++;
209       }                                           211       }
210                                                   212    
211     particleGun->SetParticleEnergy((*energies)    213     particleGun->SetParticleEnergy((*energies)[j]);
212                                                   214   
213   }                                               215   }
214       else if (particle == "gamma")               216       else if (particle == "gamma")
215   {                                               217   {
216     const XrayFluoDataSet* dataSet = runManage    218     const XrayFluoDataSet* dataSet = runManager->GetGammaSet();
217                                                   219     
218     G4int i = 0;                                  220     G4int i = 0;
219     G4int id = 0;                                 221     G4int id = 0;
220     G4double minEnergy = 0. * keV;                222     G4double minEnergy = 0. * keV;
221     G4double particleEnergy= 0.;                  223     G4double particleEnergy= 0.;
222     G4double maxEnergy = 10. * keV;               224     G4double maxEnergy = 10. * keV;
223     G4double energyRange = maxEnergy - minEner    225     G4double energyRange = maxEnergy - minEnergy;
224                                                   226 
225      while ( i == 0)                              227      while ( i == 0)
226       {                                           228       {
227         G4double random = G4UniformRand();        229         G4double random = G4UniformRand();
228                                                   230         
229         G4double randomNum = G4UniformRand();     231         G4double randomNum = G4UniformRand(); //*5.0E6;
230                                                   232         
231         particleEnergy = (random*energyRange)     233         particleEnergy = (random*energyRange) + minEnergy;
232                                                   234         
233         if ((dataSet->FindValue(particleEnergy    235         if ((dataSet->FindValue(particleEnergy,id)) > randomNum)
234     {                                             236     {
235       i = 1;                                      237       i = 1;
236                                                   238       
237     }                                             239     }
238       }                                           240       }
239      particleGun->SetParticleEnergy(particleEn    241      particleGun->SetParticleEnergy(particleEnergy);
240   }                                               242   }
241     }                                             243     }
242                                                   244   
243                                                   245 
244 #ifdef G4ANALYSIS_USE                             246 #ifdef G4ANALYSIS_USE 
245                                                   247 
246   G4double partEnergy = particleGun->GetPartic    248   G4double partEnergy = particleGun->GetParticleEnergy();
247   XrayFluoAnalysisManager* analysis =  XrayFlu    249   XrayFluoAnalysisManager* analysis =  XrayFluoAnalysisManager::getInstance();
248   analysis->analysePrimaryGenerator(partEnergy    250   analysis->analysePrimaryGenerator(partEnergy/keV);
249                                                   251 
250 #endif                                            252 #endif
251                                                   253 
252   particleGun->GeneratePrimaryVertex(anEvent);    254   particleGun->GeneratePrimaryVertex(anEvent);
253 }                                                 255 }
254                                                   256 
255 //....oooOO0OOooo........oooOO0OOooo........oo    257 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
256                                                   258 
257                                                   259 
258                                                   260 
259                                                   261 
260                                                   262 
261                                                   263 
262                                                   264 
263                                                   265 
264                                                   266 
265                                                   267