LaserAlignmentSimulation/src/LaserBeamsBarrel.cc

0001 /** \file LaserBeamsBarrel.cc
0002  *
0003  *
0004  *  $Date: 2010/09/09 18:22:48 $
0005  *  $Revision: 1.7 $
0006  *  \author Maarten Thomas
0007  */
0008
0009 #include "Alignment/LaserAlignmentSimulation/interface/LaserBeamsBarrel.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 "G4ParticleDefinition.hh"
0017 #include "G4ParticleGun.hh"
0018 #include "G4SystemOfUnits.hh"
0019 #include "globals.hh"  // Global Constants and typedefs
0020
0021 LaserBeamsBarrel::LaserBeamsBarrel() : theParticleGun(nullptr), theDRand48Engine(nullptr) {
0022   G4int nPhotonsGun = 1;
0023   G4int nPhotonsBeam = 1;
0024   G4double Energy = 1.15 * eV;
0025   // call constructor with options
0026   LaserBeamsBarrel(nPhotonsGun, nPhotonsBeam, Energy);
0027 }
0028
0029 LaserBeamsBarrel::LaserBeamsBarrel(G4int nPhotonsInGun, G4int nPhotonsInBeam, G4double PhotonEnergy)
0030     : thenParticleInGun(0), thenParticle(0), thePhotonEnergy(0), theParticleGun(), theDRand48Engine() {
0031   /* *********************************************************************** */
0032   /*  initialize and configure the particle gun                              */
0033   /* *********************************************************************** */
0034
0035   // the Photon energy
0036   thePhotonEnergy = PhotonEnergy;
0037
0038   // number of particles in the Laser beam
0039   thenParticleInGun = nPhotonsInGun;
0040
0041   // number of particles in one beam. ATTENTION: each beam contains
0042   // nParticleInGun with the same startpoint and direction. nParticle gives the
0043   // number of particles in the beam with a different startpoint. They are used
0044   // to simulate the gaussian beamprofile of the Laser Beams.
0045   thenParticle = nPhotonsInBeam;
0046
0047   // create the particle gun
0048   theParticleGun = new G4ParticleGun(thenParticleInGun);
0049
0050   // default kinematics
0051   G4ParticleTable *theParticleTable = G4ParticleTable::GetParticleTable();
0052   G4ParticleDefinition *theOpticalPhoton = theParticleTable->FindParticle("opticalphoton");
0053
0054   theParticleGun->SetParticleDefinition(theOpticalPhoton);
0055   theParticleGun->SetParticleTime(0.0 * ns);
0056   theParticleGun->SetParticlePosition(G4ThreeVector(-500.0 * cm, 0.0 * cm, 0.0 * cm));
0057   theParticleGun->SetParticleMomentumDirection(G4ThreeVector(5.0, 3.0, 0.0));
0058   theParticleGun->SetParticleEnergy(10.0 * keV);
0059   setOptPhotonPolar(90.0);
0060
0061   // initialize the random number engine
0062   theDRand48Engine = new CLHEP::DRand48Engine();
0063 }
0064
0065 LaserBeamsBarrel::~LaserBeamsBarrel() {
0066   if (theParticleGun != nullptr) {
0067     delete theParticleGun;
0068   }
0069   if (theDRand48Engine != nullptr) {
0070     delete theDRand48Engine;
0071   }
0072 }
0073
0074 void LaserBeamsBarrel::GeneratePrimaries(G4Event *myEvent) {
0075   // this function is called at the beginning of an Event in
0076   // LaserAlignment::upDate(const BeginOfEvent * myEvent)
0077
0078   // use the random number generator service of the framework
0079   edm::Service<edm::RandomNumberGenerator> rng;
0080   unsigned int seed = rng->mySeed();
0081
0082   // set the seed
0083   theDRand48Engine->setSeed(seed);
0084
0085   // number of LaserBeams
0086   const G4int nLaserBeams = 8;
0087
0088   // z position of the Laserdiodes (value from design drawings)
0089   G4double LaserPositionZ = 1137.0 * mm;
0090
0091   // Radius of the Laser ring
0092   G4double LaserRingRadius = 564.0 * mm;
0093
0094   // phi positions of the Laserdiodes (from CMS Note 2001/053 or from
0095   // http://abbaneo.home.cern.ch/abbaneo/cms/layout)
0096   G4double LaserPhi[nLaserBeams] = {G4double(7.0 / 112.0) * G4double(2.0 * M_PI),
0097                                     G4double(23.0 / 112.0) * G4double(2.0 * M_PI),
0098                                     G4double(33.0 / 112.0) * G4double(2.0 * M_PI),
0099                                     G4double(49.0 / 112.0) * G4double(2.0 * M_PI),
0100                                     G4double(65.0 / 112.0) * G4double(2.0 * M_PI),
0101                                     G4double(77.0 / 112.0) * G4double(2.0 * M_PI),
0102                                     G4double(93.0 / 112.0) * G4double(2.0 * M_PI),
0103                                     G4double(103.0 / 112.0) * G4double(2.0 * M_PI)};
0104
0105   // width of the LaserBeams
0106   G4double LaserBeamSigmaX = 0.5 * mm;
0107   G4double LaserBeamSigmaY = 0.5 * mm;
0108
0109   // get the definition of the optical photon
0110   G4ParticleTable *theParticleTable = G4ParticleTable::GetParticleTable();
0111   G4ParticleDefinition *theOpticalPhoton = theParticleTable->FindParticle("opticalphoton");
0112
0113   // loop over the LaserBeams
0114   for (int theBeam = 0; theBeam < nLaserBeams; theBeam++) {
0115     // code for forward and backward beam
0116     // calculate x and y position of the current laser diode
0117     G4double LaserPositionX = cos(LaserPhi[theBeam]) * LaserRingRadius;
0118     G4double LaserPositionY = sin(LaserPhi[theBeam]) * LaserRingRadius;
0119
0120     // loop over all the particles in one beam
0121     for (int theParticle = 0; theParticle < thenParticle; theParticle++) {
0122       // get randomnumbers  and calculate the position
0123       CLHEP::RandGaussQ aGaussObjX(*theDRand48Engine, LaserPositionX, LaserBeamSigmaX);
0124       CLHEP::RandGaussQ aGaussObjY(*theDRand48Engine, LaserPositionY, LaserBeamSigmaY);
0125
0126       G4double theXPosition = aGaussObjX.fire();
0127       G4double theYPosition = aGaussObjY.fire();
0128       G4double theZPosition = LaserPositionZ;
0129
0130       // set the properties of the newly created particle
0131       theParticleGun->SetParticleDefinition(theOpticalPhoton);
0132       theParticleGun->SetParticleTime(0.0 * ns);
0133       theParticleGun->SetParticlePosition(G4ThreeVector(theXPosition, theYPosition, theZPosition));
0134       theParticleGun->SetParticleEnergy(thePhotonEnergy);
0135
0136       // loop over both directions of the beam
0137       for (int theDirection = 0; theDirection < 2; theDirection++) {
0138         // shoot in both beam directions ...
0139         if (theDirection == 0)  // shoot in forward direction (+z)
0140         {
0141           theParticleGun->SetParticleMomentumDirection(G4ThreeVector(0.0, 0.0, 1.0));
0142           // set the polarization
0143           setOptPhotonPolar(90.0);
0144           // generate the particle
0145           theParticleGun->GeneratePrimaryVertex(myEvent);
0146         } else if (theDirection == 1)  // shoot in backward direction (-z)
0147         {
0148           theParticleGun->SetParticleMomentumDirection(G4ThreeVector(0.0, 0.0, -1.0));
0149           // set the polarization
0150           setOptPhotonPolar(90.0);
0151           // generate the particle
0152           theParticleGun->GeneratePrimaryVertex(myEvent);
0153         }
0154       }  // end looop over both beam directions
0155     }    // end looop over particles in beam
0156   }      // end loop over beams
0157 }
0158
0159 void LaserBeamsBarrel::setOptPhotonPolar(G4double Angle) {
0160   /* *********************************************************************** */
0161   /*   to get optical processes working properly, you have to make sure      *
0162    *   that the photon polarisation is defined.                              */
0163   /* *********************************************************************** */
0164
0165   // first check if we have an optical photon
0166   if (theParticleGun->GetParticleDefinition()->GetParticleName() != "opticalphoton") {
0167     edm::LogWarning("SimLaserAlignment:LaserBeamsBarrel")
0168         << "<LaserBeamsBarrel::setOptPhotonPolar()>: WARNING! The ParticleGun "
0169            "is not an optical photon";
0170     return;
0171   }
0172
0173   //   G4cout << "  AC1CMS: The ParticleGun is an " <<
0174   //   theParticleGun->GetParticleDefinition()->GetParticleName();
0175   G4ThreeVector normal(1.0, 0.0, 0.0);
0176   G4ThreeVector kphoton = theParticleGun->GetParticleMomentumDirection();
0177   G4ThreeVector product = normal.cross(kphoton);
0178   G4double modul2 = product * product;
0179
0180   G4ThreeVector e_perpendicular(0.0, 0.0, 1.0);
0181
0182   if (modul2 > 0.0) {
0183     e_perpendicular = (1.0 / sqrt(modul2)) * product;
0184   }
0185
0186   G4ThreeVector e_parallel = e_perpendicular.cross(kphoton);
0187
0188   G4ThreeVector polar = cos(Angle) * e_parallel + sin(Angle) * e_perpendicular;
0189
0190   //   G4cout << ", the polarization = " << polar << G4endl;
0191   theParticleGun->SetParticlePolarization(polar);
0192 }