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 8.0)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                    <<   3 // * DISCLAIMER                                                       *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th <<   5 // * The following disclaimer summarizes all the specific disclaimers *
  6 // * the Geant4 Collaboration.  It is provided <<   6 // * of contributors to this software. The specific disclaimers,which *
  7 // * conditions of the Geant4 Software License <<   7 // * govern, are listed with their locations in:                      *
  8 // * LICENSE and available at  http://cern.ch/ <<   8 // *   http://cern.ch/geant4/license                                  *
  9 // * include a list of copyright holders.      << 
 10 // *                                                9 // *                                                                  *
 11 // * Neither the authors of this software syst     10 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     11 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     12 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     13 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  <<  14 // * use.                                                             *
 16 // * for the full disclaimer and the limitatio << 
 17 // *                                               15 // *                                                                  *
 18 // * This  code  implementation is the result  <<  16 // * This  code  implementation is the  intellectual property  of the *
 19 // * technical work of the GEANT4 collaboratio <<  17 // * GEANT4 collaboration.                                            *
 20 // * By using,  copying,  modifying or  distri <<  18 // * By copying,  distributing  or modifying the Program (or any work *
 21 // * any work based  on the software)  you  ag <<  19 // * based  on  the Program)  you indicate  your  acceptance of  this *
 22 // * use  in  resulting  scientific  publicati <<  20 // * statement, and all its terms.                                    *
 23 // * acceptance of all terms of the Geant4 Sof << 
 24 // *******************************************     21 // ********************************************************************
 25 //                                                 22 //
 26 //                                                 23 //
                                                   >>  24 // $Id: GammaRayTelPrimaryGeneratorAction.cc,v 1.8 2005/06/27 15:27:33 gunter Exp $
                                                   >>  25 // GEANT4 tag $Name: geant4-08-00 $
 27 // -------------------------------------------     26 // ------------------------------------------------------------
 28 //      GEANT 4 class implementation file          27 //      GEANT 4 class implementation file
 29 //      CERN Geneva Switzerland                    28 //      CERN Geneva Switzerland
 30 //                                                 29 //
 31 //                                                 30 //
 32 //      ------------ GammaRayTelPrimaryGenerat     31 //      ------------ GammaRayTelPrimaryGeneratorAction  ------
 33 //           by  G.Santin, F.Longo & R.Giannit     32 //           by  G.Santin, F.Longo & R.Giannitrapani (13 nov 2000)
 34 //                                                 33 //
 35 // *******************************************     34 // ************************************************************
 36                                                    35 
                                                   >>  36 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
                                                   >>  37 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
                                                   >>  38 
 37 #include "G4RunManager.hh"                         39 #include "G4RunManager.hh"
 38 #include "GammaRayTelPrimaryGeneratorAction.hh     40 #include "GammaRayTelPrimaryGeneratorAction.hh"
 39                                                    41 
 40 #include "GammaRayTelDetectorConstruction.hh"      42 #include "GammaRayTelDetectorConstruction.hh"
 41 #include "GammaRayTelPrimaryGeneratorMessenger     43 #include "GammaRayTelPrimaryGeneratorMessenger.hh"
 42                                                    44 
 43 #include "G4PhysicalConstants.hh"              << 
 44 #include "G4SystemOfUnits.hh"                  << 
 45 #include "G4Event.hh"                              45 #include "G4Event.hh"
 46 #include "G4ParticleGun.hh"                        46 #include "G4ParticleGun.hh"
 47 #include "G4GeneralParticleSource.hh"          << 
 48 #include "G4ParticleTable.hh"                      47 #include "G4ParticleTable.hh"
 49 #include "G4ParticleDefinition.hh"                 48 #include "G4ParticleDefinition.hh"
 50 #include "Randomize.hh"                            49 #include "Randomize.hh"
 51                                                    50 
 52 //....oooOO0OOooo........oooOO0OOooo........oo     51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 53                                                    52 
 54 GammaRayTelPrimaryGeneratorAction::GammaRayTel <<  53 GammaRayTelPrimaryGeneratorAction::GammaRayTelPrimaryGeneratorAction()
 55     detector = static_cast<const GammaRayTelDe <<  54   :rndmFlag("off"),nSourceType(0),nSpectrumType(0)
 56                                                <<  55 {
 57     // create a messenger for this class       <<  56   G4RunManager* runManager = G4RunManager::GetRunManager();
 58                                                <<  57   GammaRayTelDetector =
 59     gunMessenger = new GammaRayTelPrimaryGener <<  58     (GammaRayTelDetectorConstruction*)(runManager->GetUserDetectorConstruction());
                                                   >>  59 
                                                   >>  60   G4int n_particle = 1;
                                                   >>  61   particleGun  = new G4ParticleGun(n_particle);
                                                   >>  62   
                                                   >>  63   //create a messenger for this class
                                                   >>  64   gunMessenger = new GammaRayTelPrimaryGeneratorMessenger(this);
                                                   >>  65   
                                                   >>  66   // default particle kinematic
                                                   >>  67 
                                                   >>  68   G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
                                                   >>  69   G4String particleName;
                                                   >>  70   G4ParticleDefinition* particle
                                                   >>  71     = particleTable->FindParticle(particleName="e-");
                                                   >>  72   particleGun->SetParticleDefinition(particle);
                                                   >>  73   particleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,-1.));
                                                   >>  74   particleGun->SetParticleEnergy(30.*MeV);
                                                   >>  75   G4double position = 0.5*(GammaRayTelDetector->GetWorldSizeZ());
                                                   >>  76   particleGun->SetParticlePosition(G4ThreeVector(0.*cm,0.*cm,position));
                                                   >>  77   
                                                   >>  78 }
 60                                                    79 
 61     constexpr auto NUMBER_OF_PARTICLES{1};     <<  80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 62     particleGun = new G4ParticleGun(NUMBER_OF_ << 
 63                                                    81 
 64     // default particle kinematic              <<  82 GammaRayTelPrimaryGeneratorAction::~GammaRayTelPrimaryGeneratorAction()
                                                   >>  83 {
                                                   >>  84   delete particleGun;
                                                   >>  85   delete gunMessenger;
                                                   >>  86 }
 65                                                    87 
 66     auto *particleTable = G4ParticleTable::Get <<  88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 67     auto *particle = particleTable->FindPartic << 
 68     particleGun->SetParticleDefinition(particl << 
 69     particleGun->SetParticleMomentumDirection( << 
 70                                                    89 
 71     constexpr auto PARTICLE_ENERGY{30. * MeV}; <<  90 void GammaRayTelPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
 72     particleGun->SetParticleEnergy(PARTICLE_EN <<  91 {
                                                   >>  92   //this function is called at the begining of event
                                                   >>  93   // 
                                                   >>  94   G4double z0 = 0.5*(GammaRayTelDetector->GetWorldSizeZ());
                                                   >>  95   G4double x0 = 0.*cm, y0 = 0.*cm;
                                                   >>  96   
                                                   >>  97   G4ThreeVector pos0;
                                                   >>  98   G4ThreeVector dir0;
                                                   >>  99   G4ThreeVector vertex0 = G4ThreeVector(x0,y0,z0);
                                                   >> 100   
                                                   >> 101   dir0 = G4ThreeVector(0.,0.,-1.);
                                                   >> 102 
                                                   >> 103   G4double theta, phi, y, f;
                                                   >> 104   G4double theta0=0.;
                                                   >> 105   G4double phi0=0.;
                                                   >> 106   
                                                   >> 107   switch(nSourceType) {
                                                   >> 108   case 0:
                                                   >> 109     particleGun->SetParticlePosition(vertex0);
                                                   >> 110     particleGun->SetParticleMomentumDirection(dir0);
                                                   >> 111     break;
                                                   >> 112   case 1:
                                                   >> 113     // GS: Generate random position on the 4PIsphere to create a unif. distrib.
                                                   >> 114     // GS: on the sphere
                                                   >> 115     phi = G4UniformRand() * twopi;
                                                   >> 116     do {
                                                   >> 117       y = G4UniformRand()*1.0;
                                                   >> 118       theta = G4UniformRand() * pi;
                                                   >> 119       f = std::sin(theta);
                                                   >> 120     } while (y > f);
                                                   >> 121     vertex0 = G4ThreeVector(1.,0.,0.);
                                                   >> 122     vertex0.setMag(dVertexRadius);
                                                   >> 123     vertex0.setTheta(theta);
                                                   >> 124     vertex0.setPhi(phi);
                                                   >> 125     particleGun->SetParticlePosition(vertex0);
                                                   >> 126 
                                                   >> 127     dir0 = G4ThreeVector(1.,0.,0.);
                                                   >> 128     do {
                                                   >> 129       phi = G4UniformRand() * twopi;
                                                   >> 130       do {
                                                   >> 131   y = G4UniformRand()*1.0;
                                                   >> 132   theta = G4UniformRand() * pi;
                                                   >> 133   f = std::sin(theta);
                                                   >> 134       } while (y > f);
                                                   >> 135       dir0.setPhi(phi);
                                                   >> 136       dir0.setTheta(theta);
                                                   >> 137     } while (vertex0.dot(dir0) >= -0.7 * vertex0.mag());
                                                   >> 138     particleGun->SetParticleMomentumDirection((G4ParticleMomentum)dir0);
                                                   >> 139 
                                                   >> 140     break;
                                                   >> 141   case 2:
                                                   >> 142     // GS: Generate random position on the upper semi-sphere z>0 to create a unif. distrib.
                                                   >> 143     // GS: on a plane
                                                   >> 144     phi = G4UniformRand() * twopi;
                                                   >> 145     do {
                                                   >> 146       y = G4UniformRand()*1.0;
                                                   >> 147       theta = G4UniformRand() * halfpi;
                                                   >> 148       f = std::sin(theta) * std::cos(theta);
                                                   >> 149     } while (y > f);
                                                   >> 150     vertex0 = G4ThreeVector(1.,0.,0.);
                                                   >> 151     
                                                   >> 152     G4double xy = GammaRayTelDetector->GetWorldSizeXY();
                                                   >> 153     G4double z = GammaRayTelDetector->GetWorldSizeZ();
                                                   >> 154     
                                                   >> 155     if (dVertexRadius > xy*0.5)
                                                   >> 156       { 
                                                   >> 157   G4cout << "vertexRadius too big " << G4endl;
                                                   >> 158   G4cout << "vertexRadius setted to " << xy*0.45 << G4endl;
                                                   >> 159   dVertexRadius = xy*0.45;
                                                   >> 160       }
                                                   >> 161     
                                                   >> 162     if (dVertexRadius > z*0.5)
                                                   >> 163       { 
                                                   >> 164   G4cout << "vertexRadius too high " << G4endl;
                                                   >> 165   G4cout << "vertexRadius setted to " << z*0.45 << G4endl;
                                                   >> 166   dVertexRadius = z*0.45;
                                                   >> 167       }
                                                   >> 168     
                                                   >> 169 
                                                   >> 170     vertex0.setMag(dVertexRadius);
                                                   >> 171     vertex0.setTheta(theta);
                                                   >> 172     vertex0.setPhi(phi);
                                                   >> 173     
                                                   >> 174     // GS: Get the user defined direction for the primaries and
                                                   >> 175     // GS: Rotate the random position according to the user defined direction for the particle
                                                   >> 176 
                                                   >> 177     dir0 = particleGun->GetParticleMomentumDirection();
                                                   >> 178     if (dir0.mag() > 0.001) 
                                                   >> 179       {
                                                   >> 180   theta0 = dir0.theta();
                                                   >> 181   phi0   = dir0.phi();   
                                                   >> 182       }
                                                   >> 183     
                                                   >> 184     if (theta0!=0.) 
                                                   >> 185       {
                                                   >> 186   G4ThreeVector rotationAxis(1.,0.,0.);
                                                   >> 187   rotationAxis.setPhi(phi0+halfpi);
                                                   >> 188   vertex0.rotate(theta0+pi,rotationAxis);
                                                   >> 189       }
                                                   >> 190     particleGun->SetParticlePosition(vertex0);
                                                   >> 191     break;
                                                   >> 192   }
                                                   >> 193 
                                                   >> 194   
                                                   >> 195   G4double pEnergy;
                                                   >> 196   
                                                   >> 197   switch(nSpectrumType) {
                                                   >> 198   case 0:
                                                   >> 199     break;
                                                   >> 200   case 1:
                                                   >> 201     break;
                                                   >> 202   case 2:
                                                   >> 203     do {
                                                   >> 204       y = G4UniformRand()*100000.0;
                                                   >> 205       pEnergy = G4UniformRand() * 10. * GeV;
                                                   >> 206       f = std::pow(pEnergy * (1/GeV), -4.);
                                                   >> 207     } while (y > f);
                                                   >> 208     
                                                   >> 209     particleGun->SetParticleEnergy(pEnergy);
                                                   >> 210     
                                                   >> 211     break;
                                                   >> 212   case 3:
                                                   >> 213     break;
                                                   >> 214   }
 73                                                   215 
 74     auto position = 0.5 * (detector->GetWorldS << 216   particleGun->GeneratePrimaryVertex(anEvent);
 75     particleGun->SetParticlePosition(G4ThreeVe << 217   
 76     particleSource = new G4GeneralParticleSour << 
 77 }                                                 218 }
 78                                                   219 
 79 //....oooOO0OOooo........oooOO0OOooo........oo    220 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 80                                                   221 
 81 GammaRayTelPrimaryGeneratorAction::~GammaRayTe << 
 82     delete particleGun;                        << 
 83     delete particleSource;                     << 
 84     delete gunMessenger;                       << 
 85 }                                              << 
 86                                                   222 
 87 //....oooOO0OOooo........oooOO0OOooo........oo << 
 88                                                   223 
 89 void GammaRayTelPrimaryGeneratorAction::Genera << 224 
 90     if (sourceGun) {                           << 225 
 91         G4cout << "Using G4ParticleGun... " << << 226 
 92                                                << 227 
 93         // this function is called at the begi << 228 
 94         //                                     << 
 95         G4double x0 = 0. * cm;                 << 
 96         G4double y0 = 0. * cm;                 << 
 97         G4double z0 = 0.5 * (detector->GetWorl << 
 98                                                << 
 99         G4ThreeVector pos0;                    << 
100         auto vertex0 = G4ThreeVector(x0, y0, z << 
101         auto momentumDirection0 = G4ThreeVecto << 
102                                                << 
103         G4double theta;                        << 
104         G4double phi;                          << 
105         G4double y = 0.;                       << 
106         G4double f = 0.;                       << 
107         G4double theta0 = 0.;                  << 
108         G4double phi0 = 0.;                    << 
109                                                << 
110         switch (sourceType) {                  << 
111         case 0:                                << 
112             particleGun->SetParticlePosition(v << 
113             particleGun->SetParticleMomentumDi << 
114             break;                             << 
115         case 1:                                << 
116             // GS: Generate random position on << 
117             // GS: on the sphere               << 
118             phi = G4UniformRand() * twopi;     << 
119             do {                               << 
120                 y = G4UniformRand() * 1.0;     << 
121                 theta = G4UniformRand() * pi;  << 
122                 f = std::sin(theta);           << 
123             } while (y > f);                   << 
124             vertex0 = G4ThreeVector(1., 0., 0. << 
125             vertex0.setMag(vertexRadius);      << 
126             vertex0.setTheta(theta);           << 
127             vertex0.setPhi(phi);               << 
128             particleGun->SetParticlePosition(v << 
129                                                << 
130             momentumDirection0 = G4ThreeVector << 
131                                                << 
132             do {                               << 
133                 phi = G4UniformRand() * twopi; << 
134                 do {                           << 
135                     y = G4UniformRand() * 1.0; << 
136                     theta = G4UniformRand() *  << 
137                     f = std::sin(theta);       << 
138                 } while (y > f);               << 
139                 momentumDirection0.setPhi(phi) << 
140                 momentumDirection0.setTheta(th << 
141             } while (vertex0.dot(momentumDirec << 
142                                                << 
143             particleGun->SetParticleMomentumDi << 
144                                                << 
145             break;                             << 
146         case 2:                                << 
147             // GS: Generate random position on << 
148             // GS: on a plane                  << 
149             phi = G4UniformRand() * twopi;     << 
150                                                << 
151             do {                               << 
152                 y = G4UniformRand() * 1.0;     << 
153                 theta = G4UniformRand() * half << 
154                 f = std::sin(theta) * std::cos << 
155             } while (y > f);                   << 
156                                                << 
157             vertex0 = G4ThreeVector(1., 0., 0. << 
158                                                << 
159             auto xy = detector->GetWorldSizeXY << 
160             auto z = detector->GetWorldSizeZ() << 
161                                                << 
162             if (vertexRadius > xy * 0.5) {     << 
163                 G4cout << "vertexRadius too bi << 
164                 G4cout << "vertexRadius set to << 
165                 vertexRadius = xy * 0.45;      << 
166             }                                  << 
167                                                << 
168             if (vertexRadius > z * 0.5) {      << 
169                 G4cout << "vertexRadius too hi << 
170                 G4cout << "vertexRadius set to << 
171                 vertexRadius = z * 0.45;       << 
172             }                                  << 
173                                                << 
174             vertex0.setMag(vertexRadius);      << 
175             vertex0.setTheta(theta);           << 
176             vertex0.setPhi(phi);               << 
177                                                << 
178             // GS: Get the user defined direct << 
179             // GS: Rotate the random position  << 
180                                                << 
181             momentumDirection0 = particleGun-> << 
182             if (momentumDirection0.mag() > 0.0 << 
183                 theta0 = momentumDirection0.th << 
184                 phi0 = momentumDirection0.phi( << 
185             }                                  << 
186                                                << 
187             if (theta0 != 0.) {                << 
188                 G4ThreeVector rotationAxis(1., << 
189                 rotationAxis.setPhi(phi0 + hal << 
190                 vertex0.rotate(theta0 + pi, ro << 
191             }                                  << 
192             particleGun->SetParticlePosition(v << 
193             break;                             << 
194         }                                      << 
195                                                << 
196         constexpr auto INITIAL_PARTICLE_ENERGY << 
197         G4double particleEnergy = INITIAL_PART << 
198                                                << 
199         switch (spectrumType) {                << 
200         case 0: // Uniform energy (1 GeV - 10  << 
201             y = G4UniformRand();               << 
202             particleEnergy = y * 9.0 * GeV + 1 << 
203             G4cout << "Particle energy: " << p << 
204             break;                             << 
205         case 1: // Logarithmic energy          << 
206             y = G4UniformRand();               << 
207             particleEnergy = std::pow(10, y) * << 
208             G4cout << "Particle energy: " << p << 
209             break;                             << 
210         case 2: // Power law (-4)              << 
211             do {                               << 
212                 y = G4UniformRand() * 100000.0 << 
213                 particleEnergy = G4UniformRand << 
214                 f = std::pow(particleEnergy *  << 
215             } while (y > f);                   << 
216             // particleGun->SetParticleEnergy( << 
217             break;                             << 
218         case 3: // Monochromatic               << 
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     }                                          << 
230 }                                              << 
231                                                   229