File indexing completed on 2024-09-07 04:34:31
0001
0002
0003
0004
0005
0006
0007
0008
0009 #include "Alignment/LaserAlignmentSimulation/interface/LaserBeamsTEC2.h"
0010
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 #include "FWCore/ServiceRegistry/interface/Service.h"
0013 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0014
0015 #include <CLHEP/Random/RandGaussQ.h>
0016 #include <CLHEP/Units/SystemOfUnits.h>
0017 #include "G4ParticleDefinition.hh"
0018 #include "G4ParticleGun.hh"
0019 #include "globals.hh" // Global Constants and typedefs
0020
0021 LaserBeamsTEC2::LaserBeamsTEC2() : theParticleGun(nullptr), theDRand48Engine(nullptr) {
0022 G4int nPhotonsGun = 1;
0023 G4int nPhotonsBeam = 1;
0024 G4double Energy = 1.15 * CLHEP::eV;
0025
0026 LaserBeamsTEC2(nPhotonsGun, nPhotonsBeam, Energy);
0027 }
0028
0029 LaserBeamsTEC2::LaserBeamsTEC2(G4int nPhotonsInGun, G4int nPhotonsInBeam, G4double PhotonEnergy)
0030 : thenParticleInGun(0), thenParticle(0), thePhotonEnergy(0), theParticleGun(), theDRand48Engine() {
0031
0032
0033
0034
0035
0036 thePhotonEnergy = PhotonEnergy;
0037
0038
0039 thenParticleInGun = nPhotonsInGun;
0040
0041
0042
0043
0044
0045 thenParticle = nPhotonsInBeam;
0046
0047
0048 theParticleGun = new G4ParticleGun(thenParticleInGun);
0049
0050
0051 G4ParticleTable *theParticleTable = G4ParticleTable::GetParticleTable();
0052 G4ParticleDefinition *theOpticalPhoton = theParticleTable->FindParticle("opticalphoton");
0053
0054 theParticleGun->SetParticleDefinition(theOpticalPhoton);
0055 theParticleGun->SetParticleTime(0.0 * CLHEP::ns);
0056 theParticleGun->SetParticlePosition(G4ThreeVector(-500.0 * CLHEP::cm, 0.0 * CLHEP::cm, 0.0 * CLHEP::cm));
0057 theParticleGun->SetParticleMomentumDirection(G4ThreeVector(5.0, 3.0, 0.0));
0058 theParticleGun->SetParticleEnergy(10.0 * CLHEP::keV);
0059 setOptPhotonPolar(90.0);
0060
0061
0062 theDRand48Engine = new CLHEP::DRand48Engine();
0063 }
0064
0065 LaserBeamsTEC2::~LaserBeamsTEC2() {
0066 if (theParticleGun != nullptr) {
0067 delete theParticleGun;
0068 }
0069 if (theDRand48Engine != nullptr) {
0070 delete theDRand48Engine;
0071 }
0072 }
0073
0074 void LaserBeamsTEC2::GeneratePrimaries(G4Event *myEvent) {
0075
0076
0077
0078
0079 edm::Service<edm::RandomNumberGenerator> rng;
0080 unsigned int seed = rng->mySeed();
0081
0082
0083 theDRand48Engine->setSeed(seed);
0084
0085
0086 const G4int nLaserRings = 2;
0087 const G4int nLaserBeams = 8;
0088
0089
0090
0091 using CLHEP::mm;
0092 G4double LaserPositionZ = -2057.5 * mm;
0093
0094
0095 G4double LaserRingRadius[nLaserRings] = {564.0 * mm, 840.0 * mm};
0096
0097
0098 G4double LaserPhi0 = 0.392699082;
0099
0100
0101 G4double LaserRingSigmaX[nLaserRings] = {0.5 * mm, 0.5 * mm};
0102 G4double LaserRingSigmaY[nLaserRings] = {0.5 * mm, 0.5 * mm};
0103
0104
0105 G4ParticleTable *theParticleTable = G4ParticleTable::GetParticleTable();
0106 G4ParticleDefinition *theOpticalPhoton = theParticleTable->FindParticle("opticalphoton");
0107
0108
0109 for (int theRing = 0; theRing < nLaserRings; theRing++) {
0110
0111 for (int theBeam = 0; theBeam < nLaserBeams; theBeam++) {
0112
0113
0114 G4double LaserPositionPhi = LaserPhi0 + G4double(theBeam * G4double(G4double(2 * M_PI) / nLaserBeams));
0115
0116
0117 G4double LaserPositionX = cos(LaserPositionPhi) * LaserRingRadius[theRing];
0118 G4double LaserPositionY = sin(LaserPositionPhi) * LaserRingRadius[theRing];
0119
0120
0121 for (int theParticle = 0; theParticle < thenParticle; theParticle++) {
0122
0123 CLHEP::RandGaussQ aGaussObjX(*theDRand48Engine, LaserPositionX, LaserRingSigmaX[theRing]);
0124 CLHEP::RandGaussQ aGaussObjY(*theDRand48Engine, LaserPositionY, LaserRingSigmaY[theRing]);
0125
0126 G4double theXPosition = aGaussObjX.fire();
0127 G4double theYPosition = aGaussObjY.fire();
0128 G4double theZPosition = LaserPositionZ;
0129
0130
0131 theParticleGun->SetParticleDefinition(theOpticalPhoton);
0132 theParticleGun->SetParticleTime(0.0 * CLHEP::ns);
0133 theParticleGun->SetParticlePosition(G4ThreeVector(theXPosition, theYPosition, theZPosition));
0134 theParticleGun->SetParticleEnergy(thePhotonEnergy);
0135
0136
0137 for (int theDirection = 0; theDirection < 2; theDirection++) {
0138
0139 if (theDirection == 0)
0140 {
0141 theParticleGun->SetParticleMomentumDirection(G4ThreeVector(0.0, 0.0, 1.0));
0142
0143 setOptPhotonPolar(90.0);
0144
0145 theParticleGun->GeneratePrimaryVertex(myEvent);
0146 }
0147
0148 if (theDirection == 1)
0149 {
0150 theParticleGun->SetParticleMomentumDirection(G4ThreeVector(0.0, 0.0, -1.0));
0151
0152 setOptPhotonPolar(90.0);
0153
0154 theParticleGun->GeneratePrimaryVertex(myEvent);
0155 }
0156 }
0157 }
0158 }
0159 }
0160 }
0161
0162 void LaserBeamsTEC2::setOptPhotonPolar(G4double Angle) {
0163
0164
0165
0166
0167
0168
0169 if (theParticleGun->GetParticleDefinition()->GetParticleName() != "opticalphoton") {
0170 edm::LogWarning("SimLaserAlignment:LaserBeamsTEC2")
0171 << "<LaserBeamsTEC2::setOptPhotonPolar()>: WARNING! The ParticleGun is "
0172 "not an optical photon";
0173 return;
0174 }
0175
0176
0177
0178 G4ThreeVector normal(1.0, 0.0, 0.0);
0179 G4ThreeVector kphoton = theParticleGun->GetParticleMomentumDirection();
0180 G4ThreeVector product = normal.cross(kphoton);
0181 G4double modul2 = product * product;
0182
0183 G4ThreeVector e_perpendicular(0.0, 0.0, 1.0);
0184
0185 if (modul2 > 0.0) {
0186 e_perpendicular = (1.0 / sqrt(modul2)) * product;
0187 }
0188
0189 G4ThreeVector e_parallel = e_perpendicular.cross(kphoton);
0190
0191 G4ThreeVector polar = cos(Angle) * e_parallel + sin(Angle) * e_perpendicular;
0192
0193
0194 theParticleGun->SetParticlePolarization(polar);
0195 }