nanobeam/src/PrimaryGeneratorAction.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 // Please cite the following paper if you use this software
0027 // Nucl.Instrum.Meth.B260:20-27, 2007
0028
0029 #include "PrimaryGeneratorAction.hh"
0030 #include "PrimaryGeneratorMessenger.hh"
0031 #include "G4SystemOfUnits.hh"
0032
0033 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0034
0035 PrimaryGeneratorAction::PrimaryGeneratorAction(DetectorConstruction* DC)
0036 :fDetector(DC)
0037 {
0038   fAngleMax = 0.09;
0039
0040   // Default
0041   fEmission = 1;
0042
0043   G4int n_particle = 1;
0044   fParticleGun  = new G4ParticleGun(n_particle);
0045
0046   G4ParticleDefinition* particle =
0047     G4ParticleTable::GetParticleTable()->FindParticle("proton");
0048       fParticleGun->SetParticleDefinition(particle);
0049
0050   fParticleGun->SetParticleEnergy(3*MeV);
0051
0052   fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,1.));
0053
0054   fGunMessenger = new PrimaryGeneratorMessenger(this);
0055
0056   // Matrix
0057
0058   fBeamMatrix = CLHEP::HepMatrix(32,32);
0059
0060 }
0061
0062 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0063
0064 PrimaryGeneratorAction::~PrimaryGeneratorAction()
0065 {
0066   delete fParticleGun;
0067   delete fGunMessenger;
0068 }
0069
0070 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0071
0072 void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
0073 {
0074   G4int numEvent;
0075
0076   numEvent=anEvent->GetEventID()+1;
0077
0078   G4double x0,y0,z0,theta,phi,xMom0,yMom0,zMom0,e0,de;
0079
0080   fShoot=false;
0081
0082   x0=0; y0=0; z0=0; theta=0; phi=0; e0=0;
0083
0084 // Coefficient computation
0085
0086 if (fEmission==1)
0087 {
0088     fDetector->SetCoef(1);
0089     fShoot=true;
0090     de = 0;
0091
0092     if (numEvent==  1)  {theta =    -0.00500; phi = -0.00500; de = -0.0050; }
0093     if (numEvent==  2)  {theta =    -0.005/3; phi = -0.00500; de = -0.0050; }
0094     if (numEvent==  3)  {theta =     0.005/3; phi = -0.00500; de = -0.0050; }
0095     if (numEvent==  4)  {theta =     0.00500; phi = -0.00500; de = -0.0050; }
0096     if (numEvent==  5)  {theta =    -0.00500; phi = -0.005/3; de = -0.0050; }
0097     if (numEvent==  6)  {theta =    -0.005/3; phi = -0.005/3; de = -0.0050; }
0098     if (numEvent==  7)  {theta =     0.005/3; phi = -0.005/3; de = -0.0050; }
0099     if (numEvent==  8)  {theta =     0.00500; phi = -0.005/3; de = -0.0050; }
0100     if (numEvent==  9)  {theta =    -0.00500; phi =  0.005/3; de = -0.0050; }
0101     if (numEvent==  10) {theta =    -0.005/3; phi =  0.005/3; de = -0.0050; }
0102     if (numEvent==  11) {theta =     0.005/3; phi =  0.005/3; de = -0.0050; }
0103     if (numEvent==  12) {theta =     0.00500; phi =  0.005/3; de = -0.0050; }
0104     if (numEvent==  13) {theta =    -0.00500; phi =  0.00500; de = -0.0050; }
0105     if (numEvent==  14) {theta =    -0.005/3; phi =  0.00500; de = -0.0050; }
0106     if (numEvent==  15) {theta =     0.005/3; phi =  0.00500; de = -0.0050; }
0107     if (numEvent==  16) {theta =     0.00500; phi =  0.00500; de = -0.0050; }
0108     if (numEvent==  17) {theta =    -0.00500; phi = -0.00500; de =  0.0050; }
0109     if (numEvent==  18) {theta =    -0.005/3; phi = -0.00500; de =  0.0050; }
0110     if (numEvent==  19) {theta =     0.005/3; phi = -0.00500; de =  0.0050; }
0111     if (numEvent==  20) {theta =     0.00500; phi = -0.00500; de =  0.0050; }
0112     if (numEvent==  21) {theta =    -0.00500; phi = -0.005/3; de =  0.0050; }
0113     if (numEvent==  22) {theta =    -0.005/3; phi = -0.005/3; de =  0.0050; }
0114     if (numEvent==  23) {theta =     0.005/3; phi = -0.005/3; de =  0.0050; }
0115     if (numEvent==  24) {theta =     0.00500; phi = -0.005/3; de =  0.0050; }
0116     if (numEvent==  25) {theta =    -0.00500; phi =  0.005/3; de =  0.0050; }
0117     if (numEvent==  26) {theta =    -0.005/3; phi =  0.005/3; de =  0.0050; }
0118     if (numEvent==  27) {theta =     0.005/3; phi =  0.005/3; de =  0.0050; }
0119     if (numEvent==  28) {theta =     0.00500; phi =  0.005/3; de =  0.0050; }
0120     if (numEvent==  29) {theta =    -0.00500; phi =  0.00500; de =  0.0050; }
0121     if (numEvent==  30) {theta =    -0.005/3; phi =  0.00500; de =  0.0050; }
0122     if (numEvent==  31) {theta =     0.005/3; phi =  0.00500; de =  0.0050; }
0123     if (numEvent==  32) {theta =     0.00500; phi =  0.00500; de =  0.0050; }
0124
0125     theta=theta*200*fAngleMax*1e-3;
0126     phi=phi*200*fAngleMax*1e-3;
0127
0128     de=de/100;
0129     e0 = 3*MeV + 2*de*3*MeV;
0130
0131     x0 = 0;
0132     y0 = 0;
0133     z0 = -8870*mm;
0134
0135     // triplet only
0136     //z0 = -3230*mm;
0137
0138         G4float cte = 1 ;
0139         G4float DE = 100*de;
0140
0141         phi = phi * 1000;
0142         theta = theta * 1000;
0143
0144         // Matrix filling up
0145
0146         fBeamMatrix(numEvent,1) = cte;
0147
0148         fBeamMatrix(numEvent,2) = theta;
0149         fBeamMatrix(numEvent,3) = phi;
0150         fBeamMatrix(numEvent,4) = DE;
0151
0152         fBeamMatrix(numEvent,5) = theta * theta;
0153         fBeamMatrix(numEvent,6) = theta * phi;
0154         fBeamMatrix(numEvent,7) = phi * phi;
0155         fBeamMatrix(numEvent,8) = theta * DE;
0156         fBeamMatrix(numEvent,9) = phi * DE;
0157
0158         fBeamMatrix(numEvent,10) = theta * theta * theta;
0159         fBeamMatrix(numEvent,11) = theta * theta * phi;
0160         fBeamMatrix(numEvent,12) = theta * phi * phi;
0161         fBeamMatrix(numEvent,13) = phi *  phi *  phi;
0162         fBeamMatrix(numEvent,14) = theta * theta * DE;
0163         fBeamMatrix(numEvent,15) = theta * phi * DE;
0164         fBeamMatrix(numEvent,16) = phi *  phi * DE;
0165
0166         fBeamMatrix(numEvent,17) = theta * theta * theta * phi;
0167         fBeamMatrix(numEvent,18) = theta * theta * phi * phi;
0168         fBeamMatrix(numEvent,19) = theta * phi * phi * phi;
0169         fBeamMatrix(numEvent,20) = theta * theta * theta * DE;
0170         fBeamMatrix(numEvent,21) = theta * theta * phi *  DE;
0171         fBeamMatrix(numEvent,22) = theta * phi * phi * DE;
0172         fBeamMatrix(numEvent,23) = phi * phi * phi * DE;
0173
0174         fBeamMatrix(numEvent,24) = theta * theta * theta * phi * phi;
0175         fBeamMatrix(numEvent,25) = theta * theta * phi * phi * phi;
0176         fBeamMatrix(numEvent,26) = theta * theta * theta * phi * DE;
0177         fBeamMatrix(numEvent,27) = theta * theta * phi * phi * DE;
0178         fBeamMatrix(numEvent,28) = theta * phi * phi * phi * DE;
0179
0180         fBeamMatrix(numEvent,29) = theta * theta * theta * phi * phi * phi;
0181         fBeamMatrix(numEvent,30) = theta * theta * theta * phi * phi * DE;
0182         fBeamMatrix(numEvent,31) = theta * theta * phi * phi * phi * DE;
0183
0184         fBeamMatrix(numEvent,32) = theta * theta * theta * phi * phi * phi * DE;
0185
0186        //
0187
0188         phi = phi / 1000;
0189         theta = theta / 1000;
0190
0191         /*
0192         G4cout << fDetector->GetG1() << G4endl;
0193         G4cout << fDetector->GetG2() << G4endl;
0194         G4cout << fDetector->GetG3() << G4endl;
0195         G4cout << fDetector->GetG4() << G4endl;
0196         G4cout << fDetector->GetModel() << G4endl;
0197         G4cout << fDetector->GetCoef() << G4endl;
0198         G4cout << fDetector->GetGrid() << G4endl;
0199         */
0200
0201 } // end coefficient
0202
0203 // Full beam
0204
0205 if (fEmission==2)
0206 {
0207     fDetector->SetCoef(0);
0208     fShoot=false;
0209
0210     G4double aR, angle, rR;
0211     aR = -1;
0212
0213     e0= G4RandGauss::shoot(3*MeV,5.0955e-5*MeV);  // AIFIRA ENERGY RESOLUTION
0214
0215     while (aR < 0) aR = G4RandGauss::shoot(0.10e-3 , 0.06e-3/2.35) * rad; // old =0.08e-3 displacement
0216
0217     angle = G4UniformRand() * 2 * CLHEP::pi *rad;
0218
0219     theta = aR * std::cos(angle);
0220
0221     phi = aR * std::sin(angle);
0222
0223     rR = XYofAngle(aR);
0224
0225     x0 = rR*std::cos(angle);
0226
0227     y0 = rR*std::sin(angle);
0228
0229     x0 = G4RandGauss::shoot(x0,220/2.35) *micrometer;
0230
0231     theta = G4RandGauss::shoot(theta,0.03e-3/2.35);
0232
0233     y0 = G4RandGauss::shoot(y0,220/2.35) *micrometer;
0234
0235     phi = G4RandGauss::shoot(phi,0.03e-3/2.35);
0236
0237     z0 = (-9120+250)*mm;  //position C0
0238
0239     if (std::sqrt(x0*x0+y0*y0)/micrometer<2000) fShoot=true;
0240
0241     /*
0242     xMom0 = std::sin(theta);
0243     yMom0 = std::sin(phi);
0244     zMom0 = std::cos(theta);
0245     */
0246
0247 } // end full beam
0248
0249 // shoot
0250
0251 if (fShoot)
0252 {
0253
0254   G4cout << "-> Event= " << numEvent<< ": X0(mm)= " << x0/mm << " X0(mm) = " << y0/mm << " Z0(m) = " << z0/m
0255          << " THETA(mrad)= " << theta/mrad << " PHI(mrad)= " << phi/mrad << " E0(MeV)= " << e0/MeV << G4endl;
0256
0257   xMom0 = std::sin(theta);
0258
0259   yMom0 = std::sin(phi);
0260
0261   zMom0 = std::sqrt(1.-xMom0*xMom0-yMom0*yMom0);
0262
0263   fParticleGun->SetParticleEnergy(e0);
0264
0265   fParticleGun->SetParticleMomentumDirection(G4ThreeVector(xMom0,yMom0,zMom0));
0266
0267   fParticleGun->SetParticlePosition(G4ThreeVector(x0,y0,z0));
0268
0269   G4ParticleDefinition* particle=
0270     G4ParticleTable::GetParticleTable()->FindParticle("proton");
0271
0272   fParticleGun->SetParticleDefinition(particle);
0273
0274   fParticleGun->GeneratePrimaryVertex(anEvent);
0275 }
0276
0277 // end shoot
0278
0279 }
0280
0281 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0282
0283 void PrimaryGeneratorAction::SetEmission(G4int value)
0284 {
0285   fEmission = value;
0286 }
0287
0288 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0289
0290 G4double PrimaryGeneratorAction::XYofAngle(G4double angle)//  returns position in micrometers
0291 {
0292   return std::pow(20000*angle, 13);
0293 }
0294