Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 // Please cite the following paper if you use this software 27 // Nucl.Instrum.Meth.B260:20-27, 2007 28 29 #include "PrimaryGeneratorAction.hh" 30 #include "PrimaryGeneratorMessenger.hh" 31 #include "G4SystemOfUnits.hh" 32 33 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 34 35 PrimaryGeneratorAction::PrimaryGeneratorAction(DetectorConstruction* DC) 36 :fDetector(DC) 37 { 38 fAngleMax = 0.09; 39 40 // Default 41 fEmission = 1; 42 43 G4int n_particle = 1; 44 fParticleGun = new G4ParticleGun(n_particle); 45 46 G4ParticleDefinition* particle = 47 G4ParticleTable::GetParticleTable()->FindParticle("proton"); 48 fParticleGun->SetParticleDefinition(particle); 49 50 fParticleGun->SetParticleEnergy(3*MeV); 51 52 fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,1.)); 53 54 fGunMessenger = new PrimaryGeneratorMessenger(this); 55 56 // Matrix 57 58 fBeamMatrix = CLHEP::HepMatrix(32,32); 59 60 } 61 62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 63 64 PrimaryGeneratorAction::~PrimaryGeneratorAction() 65 { 66 delete fParticleGun; 67 delete fGunMessenger; 68 } 69 70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 71 72 void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) 73 { 74 G4int numEvent; 75 76 numEvent=anEvent->GetEventID()+1; 77 78 G4double x0,y0,z0,theta,phi,xMom0,yMom0,zMom0,e0,de; 79 80 fShoot=false; 81 82 x0=0; y0=0; z0=0; theta=0; phi=0; e0=0; 83 84 // Coefficient computation 85 86 if (fEmission==1) 87 { 88 fDetector->SetCoef(1); 89 fShoot=true; 90 de = 0; 91 92 if (numEvent== 1) {theta = -0.00500; phi = -0.00500; de = -0.0050; } 93 if (numEvent== 2) {theta = -0.005/3; phi = -0.00500; de = -0.0050; } 94 if (numEvent== 3) {theta = 0.005/3; phi = -0.00500; de = -0.0050; } 95 if (numEvent== 4) {theta = 0.00500; phi = -0.00500; de = -0.0050; } 96 if (numEvent== 5) {theta = -0.00500; phi = -0.005/3; de = -0.0050; } 97 if (numEvent== 6) {theta = -0.005/3; phi = -0.005/3; de = -0.0050; } 98 if (numEvent== 7) {theta = 0.005/3; phi = -0.005/3; de = -0.0050; } 99 if (numEvent== 8) {theta = 0.00500; phi = -0.005/3; de = -0.0050; } 100 if (numEvent== 9) {theta = -0.00500; phi = 0.005/3; de = -0.0050; } 101 if (numEvent== 10) {theta = -0.005/3; phi = 0.005/3; de = -0.0050; } 102 if (numEvent== 11) {theta = 0.005/3; phi = 0.005/3; de = -0.0050; } 103 if (numEvent== 12) {theta = 0.00500; phi = 0.005/3; de = -0.0050; } 104 if (numEvent== 13) {theta = -0.00500; phi = 0.00500; de = -0.0050; } 105 if (numEvent== 14) {theta = -0.005/3; phi = 0.00500; de = -0.0050; } 106 if (numEvent== 15) {theta = 0.005/3; phi = 0.00500; de = -0.0050; } 107 if (numEvent== 16) {theta = 0.00500; phi = 0.00500; de = -0.0050; } 108 if (numEvent== 17) {theta = -0.00500; phi = -0.00500; de = 0.0050; } 109 if (numEvent== 18) {theta = -0.005/3; phi = -0.00500; de = 0.0050; } 110 if (numEvent== 19) {theta = 0.005/3; phi = -0.00500; de = 0.0050; } 111 if (numEvent== 20) {theta = 0.00500; phi = -0.00500; de = 0.0050; } 112 if (numEvent== 21) {theta = -0.00500; phi = -0.005/3; de = 0.0050; } 113 if (numEvent== 22) {theta = -0.005/3; phi = -0.005/3; de = 0.0050; } 114 if (numEvent== 23) {theta = 0.005/3; phi = -0.005/3; de = 0.0050; } 115 if (numEvent== 24) {theta = 0.00500; phi = -0.005/3; de = 0.0050; } 116 if (numEvent== 25) {theta = -0.00500; phi = 0.005/3; de = 0.0050; } 117 if (numEvent== 26) {theta = -0.005/3; phi = 0.005/3; de = 0.0050; } 118 if (numEvent== 27) {theta = 0.005/3; phi = 0.005/3; de = 0.0050; } 119 if (numEvent== 28) {theta = 0.00500; phi = 0.005/3; de = 0.0050; } 120 if (numEvent== 29) {theta = -0.00500; phi = 0.00500; de = 0.0050; } 121 if (numEvent== 30) {theta = -0.005/3; phi = 0.00500; de = 0.0050; } 122 if (numEvent== 31) {theta = 0.005/3; phi = 0.00500; de = 0.0050; } 123 if (numEvent== 32) {theta = 0.00500; phi = 0.00500; de = 0.0050; } 124 125 theta=theta*200*fAngleMax*1e-3; 126 phi=phi*200*fAngleMax*1e-3; 127 128 de=de/100; 129 e0 = 3*MeV + 2*de*3*MeV; 130 131 x0 = 0; 132 y0 = 0; 133 z0 = -8870*mm; 134 135 // triplet only 136 //z0 = -3230*mm; 137 138 G4float cte = 1 ; 139 G4float DE = 100*de; 140 141 phi = phi * 1000; 142 theta = theta * 1000; 143 144 // Matrix filling up 145 146 fBeamMatrix(numEvent,1) = cte; 147 148 fBeamMatrix(numEvent,2) = theta; 149 fBeamMatrix(numEvent,3) = phi; 150 fBeamMatrix(numEvent,4) = DE; 151 152 fBeamMatrix(numEvent,5) = theta * theta; 153 fBeamMatrix(numEvent,6) = theta * phi; 154 fBeamMatrix(numEvent,7) = phi * phi; 155 fBeamMatrix(numEvent,8) = theta * DE; 156 fBeamMatrix(numEvent,9) = phi * DE; 157 158 fBeamMatrix(numEvent,10) = theta * theta * theta; 159 fBeamMatrix(numEvent,11) = theta * theta * phi; 160 fBeamMatrix(numEvent,12) = theta * phi * phi; 161 fBeamMatrix(numEvent,13) = phi * phi * phi; 162 fBeamMatrix(numEvent,14) = theta * theta * DE; 163 fBeamMatrix(numEvent,15) = theta * phi * DE; 164 fBeamMatrix(numEvent,16) = phi * phi * DE; 165 166 fBeamMatrix(numEvent,17) = theta * theta * theta * phi; 167 fBeamMatrix(numEvent,18) = theta * theta * phi * phi; 168 fBeamMatrix(numEvent,19) = theta * phi * phi * phi; 169 fBeamMatrix(numEvent,20) = theta * theta * theta * DE; 170 fBeamMatrix(numEvent,21) = theta * theta * phi * DE; 171 fBeamMatrix(numEvent,22) = theta * phi * phi * DE; 172 fBeamMatrix(numEvent,23) = phi * phi * phi * DE; 173 174 fBeamMatrix(numEvent,24) = theta * theta * theta * phi * phi; 175 fBeamMatrix(numEvent,25) = theta * theta * phi * phi * phi; 176 fBeamMatrix(numEvent,26) = theta * theta * theta * phi * DE; 177 fBeamMatrix(numEvent,27) = theta * theta * phi * phi * DE; 178 fBeamMatrix(numEvent,28) = theta * phi * phi * phi * DE; 179 180 fBeamMatrix(numEvent,29) = theta * theta * theta * phi * phi * phi; 181 fBeamMatrix(numEvent,30) = theta * theta * theta * phi * phi * DE; 182 fBeamMatrix(numEvent,31) = theta * theta * phi * phi * phi * DE; 183 184 fBeamMatrix(numEvent,32) = theta * theta * theta * phi * phi * phi * DE; 185 186 // 187 188 phi = phi / 1000; 189 theta = theta / 1000; 190 191 /* 192 G4cout << fDetector->GetG1() << G4endl; 193 G4cout << fDetector->GetG2() << G4endl; 194 G4cout << fDetector->GetG3() << G4endl; 195 G4cout << fDetector->GetG4() << G4endl; 196 G4cout << fDetector->GetModel() << G4endl; 197 G4cout << fDetector->GetCoef() << G4endl; 198 G4cout << fDetector->GetGrid() << G4endl; 199 */ 200 201 } // end coefficient 202 203 // Full beam 204 205 if (fEmission==2) 206 { 207 fDetector->SetCoef(0); 208 fShoot=false; 209 210 G4double aR, angle, rR; 211 aR = -1; 212 213 e0= G4RandGauss::shoot(3*MeV,5.0955e-5*MeV); // AIFIRA ENERGY RESOLUTION 214 215 while (aR < 0) aR = G4RandGauss::shoot(0.10e-3 , 0.06e-3/2.35) * rad; // old =0.08e-3 displacement 216 217 angle = G4UniformRand() * 2 * CLHEP::pi *rad; 218 219 theta = aR * std::cos(angle); 220 221 phi = aR * std::sin(angle); 222 223 rR = XYofAngle(aR); 224 225 x0 = rR*std::cos(angle); 226 227 y0 = rR*std::sin(angle); 228 229 x0 = G4RandGauss::shoot(x0,220/2.35) *micrometer; 230 231 theta = G4RandGauss::shoot(theta,0.03e-3/2.35); 232 233 y0 = G4RandGauss::shoot(y0,220/2.35) *micrometer; 234 235 phi = G4RandGauss::shoot(phi,0.03e-3/2.35); 236 237 z0 = (-9120+250)*mm; //position C0 238 239 if (std::sqrt(x0*x0+y0*y0)/micrometer<2000) fShoot=true; 240 241 /* 242 xMom0 = std::sin(theta); 243 yMom0 = std::sin(phi); 244 zMom0 = std::cos(theta); 245 */ 246 247 } // end full beam 248 249 // shoot 250 251 if (fShoot) 252 { 253 254 G4cout << "-> Event= " << numEvent<< ": X0(mm)= " << x0/mm << " X0(mm) = " << y0/mm << " Z0(m) = " << z0/m 255 << " THETA(mrad)= " << theta/mrad << " PHI(mrad)= " << phi/mrad << " E0(MeV)= " << e0/MeV << G4endl; 256 257 xMom0 = std::sin(theta); 258 259 yMom0 = std::sin(phi); 260 261 zMom0 = std::sqrt(1.-xMom0*xMom0-yMom0*yMom0); 262 263 fParticleGun->SetParticleEnergy(e0); 264 265 fParticleGun->SetParticleMomentumDirection(G4ThreeVector(xMom0,yMom0,zMom0)); 266 267 fParticleGun->SetParticlePosition(G4ThreeVector(x0,y0,z0)); 268 269 G4ParticleDefinition* particle= 270 G4ParticleTable::GetParticleTable()->FindParticle("proton"); 271 272 fParticleGun->SetParticleDefinition(particle); 273 274 fParticleGun->GeneratePrimaryVertex(anEvent); 275 } 276 277 // end shoot 278 279 } 280 281 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 282 283 void PrimaryGeneratorAction::SetEmission(G4int value) 284 { 285 fEmission = value; 286 } 287 288 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 289 290 G4double PrimaryGeneratorAction::XYofAngle(G4double angle)// returns position in micrometers 291 { 292 return std::pow(20000*angle, 13); 293 } 294 295