Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/gammaray_telescope/src/GammaRayTelPrimaryGeneratorAction.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/gammaray_telescope/src/GammaRayTelPrimaryGeneratorAction.cc (Version 11.3.0) and /examples/advanced/gammaray_telescope/src/GammaRayTelPrimaryGeneratorAction.cc (Version 9.4.p1)


  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: GammaRayTelPrimaryGeneratorAction.cc,v 1.11 2007-11-09 16:33:34 flongo Exp $
                                                   >>  28 // GEANT4 tag $Name: geant4-09-04-patch-01 $
 27 // -------------------------------------------     29 // ------------------------------------------------------------
 28 //      GEANT 4 class implementation file          30 //      GEANT 4 class implementation file
 29 //      CERN Geneva Switzerland                    31 //      CERN Geneva Switzerland
 30 //                                                 32 //
 31 //                                                 33 //
 32 //      ------------ GammaRayTelPrimaryGenerat     34 //      ------------ GammaRayTelPrimaryGeneratorAction  ------
 33 //           by  G.Santin, F.Longo & R.Giannit     35 //           by  G.Santin, F.Longo & R.Giannitrapani (13 nov 2000)
 34 //                                                 36 //
 35 // *******************************************     37 // ************************************************************
 36                                                    38 
                                                   >>  39 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
                                                   >>  40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
                                                   >>  41 
 37 #include "G4RunManager.hh"                         42 #include "G4RunManager.hh"
 38 #include "GammaRayTelPrimaryGeneratorAction.hh     43 #include "GammaRayTelPrimaryGeneratorAction.hh"
 39                                                    44 
 40 #include "GammaRayTelDetectorConstruction.hh"      45 #include "GammaRayTelDetectorConstruction.hh"
 41 #include "GammaRayTelPrimaryGeneratorMessenger     46 #include "GammaRayTelPrimaryGeneratorMessenger.hh"
 42                                                    47 
 43 #include "G4PhysicalConstants.hh"              << 
 44 #include "G4SystemOfUnits.hh"                  << 
 45 #include "G4Event.hh"                              48 #include "G4Event.hh"
 46 #include "G4ParticleGun.hh"                        49 #include "G4ParticleGun.hh"
 47 #include "G4GeneralParticleSource.hh"              50 #include "G4GeneralParticleSource.hh"
 48 #include "G4ParticleTable.hh"                      51 #include "G4ParticleTable.hh"
 49 #include "G4ParticleDefinition.hh"                 52 #include "G4ParticleDefinition.hh"
 50 #include "Randomize.hh"                            53 #include "Randomize.hh"
 51                                                    54 
 52 //....oooOO0OOooo........oooOO0OOooo........oo     55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 53                                                    56 
 54 GammaRayTelPrimaryGeneratorAction::GammaRayTel <<  57 GammaRayTelPrimaryGeneratorAction::GammaRayTelPrimaryGeneratorAction()
 55     detector = static_cast<const GammaRayTelDe <<  58   :rndmFlag("off"),nSourceType(0),nSpectrumType(0)
 56                                                <<  59 {
 57     // create a messenger for this class       <<  60   G4RunManager* runManager = G4RunManager::GetRunManager();
 58                                                <<  61   GammaRayTelDetector =
 59     gunMessenger = new GammaRayTelPrimaryGener <<  62     (GammaRayTelDetectorConstruction*)(runManager->GetUserDetectorConstruction());
 60                                                <<  63 
 61     constexpr auto NUMBER_OF_PARTICLES{1};     <<  64 
 62     particleGun = new G4ParticleGun(NUMBER_OF_ <<  65   //create a messenger for this class
 63                                                <<  66   
 64     // default particle kinematic              <<  67   gunMessenger = new GammaRayTelPrimaryGeneratorMessenger(this);
 65                                                <<  68 
 66     auto *particleTable = G4ParticleTable::Get <<  69   G4int n_particle = 1;
 67     auto *particle = particleTable->FindPartic <<  70 
 68     particleGun->SetParticleDefinition(particl <<  71   particleGun  = new G4ParticleGun(n_particle);     
 69     particleGun->SetParticleMomentumDirection( <<  72       // default particle kinematic
 70                                                <<  73   
 71     constexpr auto PARTICLE_ENERGY{30. * MeV}; <<  74   G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
 72     particleGun->SetParticleEnergy(PARTICLE_EN <<  75   G4String particleName;
 73                                                <<  76   G4ParticleDefinition* particle
 74     auto position = 0.5 * (detector->GetWorldS <<  77     = particleTable->FindParticle(particleName="e-");
 75     particleGun->SetParticlePosition(G4ThreeVe <<  78   particleGun->SetParticleDefinition(particle);
 76     particleSource = new G4GeneralParticleSour <<  79   particleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,-1.));
                                                   >>  80   particleGun->SetParticleEnergy(30.*MeV);
                                                   >>  81   G4double position = 0.5*(GammaRayTelDetector->GetWorldSizeZ());
                                                   >>  82   particleGun->SetParticlePosition(G4ThreeVector(0.*cm,0.*cm,position));
                                                   >>  83   particleSource = new G4GeneralParticleSource();
                                                   >>  84    
 77 }                                                  85 }
 78                                                    86 
 79 //....oooOO0OOooo........oooOO0OOooo........oo     87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 80                                                    88 
 81 GammaRayTelPrimaryGeneratorAction::~GammaRayTe <<  89 GammaRayTelPrimaryGeneratorAction::~GammaRayTelPrimaryGeneratorAction()
                                                   >>  90 {
                                                   >>  91   if (sourceGun)
 82     delete particleGun;                            92     delete particleGun;
                                                   >>  93   else
 83     delete particleSource;                         94     delete particleSource;
 84     delete gunMessenger;                       <<  95 
                                                   >>  96   delete gunMessenger;
 85 }                                                  97 }
 86                                                    98 
 87 //....oooOO0OOooo........oooOO0OOooo........oo     99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 88                                                   100 
 89 void GammaRayTelPrimaryGeneratorAction::Genera << 101 void GammaRayTelPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
 90     if (sourceGun) {                           << 102 {
 91         G4cout << "Using G4ParticleGun... " << << 103    if (sourceGun)
 92                                                << 104     {
 93         // this function is called at the begi << 105 
 94         //                                     << 106       //this function is called at the begining of event
 95         G4double x0 = 0. * cm;                 << 107       // 
 96         G4double y0 = 0. * cm;                 << 108       G4double z0 = 0.5*(GammaRayTelDetector->GetWorldSizeZ());
 97         G4double z0 = 0.5 * (detector->GetWorl << 109       G4double x0 = 0.*cm, y0 = 0.*cm;
 98                                                << 110   
 99         G4ThreeVector pos0;                    << 111       G4ThreeVector pos0;
100         auto vertex0 = G4ThreeVector(x0, y0, z << 112       G4ThreeVector dir0;
101         auto momentumDirection0 = G4ThreeVecto << 113       G4ThreeVector vertex0 = G4ThreeVector(x0,y0,z0);
102                                                << 114       
103         G4double theta;                        << 115       dir0 = G4ThreeVector(0.,0.,-1.);
104         G4double phi;                          << 116       
105         G4double y = 0.;                       << 117       G4double theta, phi, y, f;
106         G4double f = 0.;                       << 118       G4double theta0=0.;
107         G4double theta0 = 0.;                  << 119       G4double phi0=0.;
108         G4double phi0 = 0.;                    << 120       
109                                                << 121       switch(nSourceType) {
110         switch (sourceType) {                  << 122       case 0:
111         case 0:                                << 123   particleGun->SetParticlePosition(vertex0);
112             particleGun->SetParticlePosition(v << 124   particleGun->SetParticleMomentumDirection(dir0);
113             particleGun->SetParticleMomentumDi << 125   break;
114             break;                             << 126       case 1:
115         case 1:                                << 127   // GS: Generate random position on the 4PIsphere to create a unif. distrib.
116             // GS: Generate random position on << 128   // GS: on the sphere
117             // GS: on the sphere               << 129   phi = G4UniformRand() * twopi;
118             phi = G4UniformRand() * twopi;     << 130   do {
119             do {                               << 131     y = G4UniformRand()*1.0;
120                 y = G4UniformRand() * 1.0;     << 132     theta = G4UniformRand() * pi;
121                 theta = G4UniformRand() * pi;  << 133     f = std::sin(theta);
122                 f = std::sin(theta);           << 134   } while (y > f);
123             } while (y > f);                   << 135   vertex0 = G4ThreeVector(1.,0.,0.);
124             vertex0 = G4ThreeVector(1., 0., 0. << 136   vertex0.setMag(dVertexRadius);
125             vertex0.setMag(vertexRadius);      << 137   vertex0.setTheta(theta);
126             vertex0.setTheta(theta);           << 138   vertex0.setPhi(phi);
127             vertex0.setPhi(phi);               << 139   particleGun->SetParticlePosition(vertex0);
128             particleGun->SetParticlePosition(v << 140   
129                                                << 141   dir0 = G4ThreeVector(1.,0.,0.);
130             momentumDirection0 = G4ThreeVector << 142   do {
131                                                << 143     phi = G4UniformRand() * twopi;
132             do {                               << 144     do {
133                 phi = G4UniformRand() * twopi; << 145       y = G4UniformRand()*1.0;
134                 do {                           << 146       theta = G4UniformRand() * pi;
135                     y = G4UniformRand() * 1.0; << 147       f = std::sin(theta);
136                     theta = G4UniformRand() *  << 148     } while (y > f);
137                     f = std::sin(theta);       << 149     dir0.setPhi(phi);
138                 } while (y > f);               << 150     dir0.setTheta(theta);
139                 momentumDirection0.setPhi(phi) << 151   } while (vertex0.dot(dir0) >= -0.7 * vertex0.mag());
140                 momentumDirection0.setTheta(th << 152   particleGun->SetParticleMomentumDirection((G4ParticleMomentum)dir0);
141             } while (vertex0.dot(momentumDirec << 153   
142                                                << 154   break;
143             particleGun->SetParticleMomentumDi << 155       case 2:
144                                                << 156   // GS: Generate random position on the upper semi-sphere z>0 to create a unif. distrib.
145             break;                             << 157   // GS: on a plane
146         case 2:                                << 158   phi = G4UniformRand() * twopi;
147             // GS: Generate random position on << 159   do {
148             // GS: on a plane                  << 160     y = G4UniformRand()*1.0;
149             phi = G4UniformRand() * twopi;     << 161     theta = G4UniformRand() * halfpi;
150                                                << 162     f = std::sin(theta) * std::cos(theta);
151             do {                               << 163   } while (y > f);
152                 y = G4UniformRand() * 1.0;     << 164   vertex0 = G4ThreeVector(1.,0.,0.);
153                 theta = G4UniformRand() * half << 165   
154                 f = std::sin(theta) * std::cos << 166   G4double xy = GammaRayTelDetector->GetWorldSizeXY();
155             } while (y > f);                   << 167   G4double z = GammaRayTelDetector->GetWorldSizeZ();
156                                                << 168   
157             vertex0 = G4ThreeVector(1., 0., 0. << 169   if (dVertexRadius > xy*0.5)
158                                                << 170     { 
159             auto xy = detector->GetWorldSizeXY << 171       G4cout << "vertexRadius too big " << G4endl;
160             auto z = detector->GetWorldSizeZ() << 172       G4cout << "vertexRadius setted to " << xy*0.45 << G4endl;
161                                                << 173       dVertexRadius = xy*0.45;
162             if (vertexRadius > xy * 0.5) {     << 174     }
163                 G4cout << "vertexRadius too bi << 175   
164                 G4cout << "vertexRadius set to << 176   if (dVertexRadius > z*0.5)
165                 vertexRadius = xy * 0.45;      << 177     { 
166             }                                  << 178       G4cout << "vertexRadius too high " << G4endl;
167                                                << 179       G4cout << "vertexRadius setted to " << z*0.45 << G4endl;
168             if (vertexRadius > z * 0.5) {      << 180       dVertexRadius = z*0.45;
169                 G4cout << "vertexRadius too hi << 181     }
170                 G4cout << "vertexRadius set to << 182   
171                 vertexRadius = z * 0.45;       << 183   
172             }                                  << 184   vertex0.setMag(dVertexRadius);
173                                                << 185   vertex0.setTheta(theta);
174             vertex0.setMag(vertexRadius);      << 186   vertex0.setPhi(phi);
175             vertex0.setTheta(theta);           << 187   
176             vertex0.setPhi(phi);               << 188   // GS: Get the user defined direction for the primaries and
177                                                << 189   // GS: Rotate the random position according to the user defined direction for the particle
178             // GS: Get the user defined direct << 190   
179             // GS: Rotate the random position  << 191   dir0 = particleGun->GetParticleMomentumDirection();
180                                                << 192   if (dir0.mag() > 0.001) 
181             momentumDirection0 = particleGun-> << 193     {
182             if (momentumDirection0.mag() > 0.0 << 194       theta0 = dir0.theta();
183                 theta0 = momentumDirection0.th << 195       phi0   = dir0.phi();   
184                 phi0 = momentumDirection0.phi( << 196     }
185             }                                  << 197   
186                                                << 198   if (theta0!=0.) 
187             if (theta0 != 0.) {                << 199     {
188                 G4ThreeVector rotationAxis(1., << 200       G4ThreeVector rotationAxis(1.,0.,0.);
189                 rotationAxis.setPhi(phi0 + hal << 201       rotationAxis.setPhi(phi0+halfpi);
190                 vertex0.rotate(theta0 + pi, ro << 202       vertex0.rotate(theta0+pi,rotationAxis);
191             }                                  << 203     }
192             particleGun->SetParticlePosition(v << 204   particleGun->SetParticlePosition(vertex0);
193             break;                             << 205   break;
194         }                                      << 206       }
195                                                << 207       
196         constexpr auto INITIAL_PARTICLE_ENERGY << 208       
197         G4double particleEnergy = INITIAL_PART << 209       G4double pEnergy;
198                                                << 210       
199         switch (spectrumType) {                << 211       switch(nSpectrumType) {
200         case 0: // Uniform energy (1 GeV - 10  << 212       case 0:
201             y = G4UniformRand();               << 213   break;
202             particleEnergy = y * 9.0 * GeV + 1 << 214       case 1:
203             G4cout << "Particle energy: " << p << 215   break;
204             break;                             << 216       case 2:
205         case 1: // Logarithmic energy          << 217   do {
206             y = G4UniformRand();               << 218     y = G4UniformRand()*100000.0;
207             particleEnergy = std::pow(10, y) * << 219     pEnergy = G4UniformRand() * 10. * GeV;
208             G4cout << "Particle energy: " << p << 220     f = std::pow(pEnergy * (1/GeV), -4.);
209             break;                             << 221   } while (y > f);
210         case 2: // Power law (-4)              << 222   
211             do {                               << 223   particleGun->SetParticleEnergy(pEnergy);
212                 y = G4UniformRand() * 100000.0 << 224   
213                 particleEnergy = G4UniformRand << 225   break;
214                 f = std::pow(particleEnergy *  << 226       case 3:
215             } while (y > f);                   << 227   break;
216             // particleGun->SetParticleEnergy( << 228       }
217             break;                             << 229       
218         case 3: // Monochromatic               << 230       particleGun->GeneratePrimaryVertex(anEvent);
219             particleEnergy = particleGun->GetP << 
220             // 100 MeV;                        << 
221             G4cout << "Particle energy: " << p << 
222             break;                             << 
223         }                                      << 
224         particleGun->SetParticleEnergy(particl << 
225         G4cout << "Particle: " << particleGun- << 
226         particleGun->GeneratePrimaryVertex(eve << 
227     } else {                                   << 
228         particleSource->GeneratePrimaryVertex( << 
229     }                                             231     }
                                                   >> 232    else
                                                   >> 233      {
                                                   >> 234        particleSource->GeneratePrimaryVertex(anEvent);
                                                   >> 235      }
                                                   >> 236    
230 }                                                 237 }
                                                   >> 238 
                                                   >> 239 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
                                                   >> 240 
                                                   >> 241 
                                                   >> 242 
                                                   >> 243 
                                                   >> 244 
                                                   >> 245 
                                                   >> 246 
                                                   >> 247 
231                                                   248