Line Code
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192
/** \file LaserBeamsBarrel.cc
 *
 *
 *  $Date: 2010/09/09 18:22:48 $
 *  $Revision: 1.7 $
 *  \author Maarten Thomas
 */

#include "Alignment/LaserAlignmentSimulation/interface/LaserBeamsBarrel.h"

#include "FWCore/AbstractServices/interface/RandomNumberGenerator.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ServiceRegistry/interface/Service.h"

#include <CLHEP/Random/RandGaussQ.h>
#include "G4ParticleDefinition.hh"
#include "G4ParticleGun.hh"
#include <CLHEP/Units/SystemOfUnits.h>
#include "globals.hh"  // Global Constants and typedefs

LaserBeamsBarrel::LaserBeamsBarrel() : theParticleGun(nullptr), theDRand48Engine(nullptr) {
  G4int nPhotonsGun = 1;
  G4int nPhotonsBeam = 1;
  G4double Energy = 1.15 * CLHEP::eV;
  // call constructor with options
  LaserBeamsBarrel(nPhotonsGun, nPhotonsBeam, Energy);
}

LaserBeamsBarrel::LaserBeamsBarrel(G4int nPhotonsInGun, G4int nPhotonsInBeam, G4double PhotonEnergy)
    : thenParticleInGun(0), thenParticle(0), thePhotonEnergy(0), theParticleGun(), theDRand48Engine() {
  /* *********************************************************************** */
  /*  initialize and configure the particle gun                              */
  /* *********************************************************************** */

  // the Photon energy
  thePhotonEnergy = PhotonEnergy;

  // number of particles in the Laser beam
  thenParticleInGun = nPhotonsInGun;

  // number of particles in one beam. ATTENTION: each beam contains
  // nParticleInGun with the same startpoint and direction. nParticle gives the
  // number of particles in the beam with a different startpoint. They are used
  // to simulate the gaussian beamprofile of the Laser Beams.
  thenParticle = nPhotonsInBeam;

  // create the particle gun
  theParticleGun = new G4ParticleGun(thenParticleInGun);

  // default kinematics
  G4ParticleTable *theParticleTable = G4ParticleTable::GetParticleTable();
  G4ParticleDefinition *theOpticalPhoton = theParticleTable->FindParticle("opticalphoton");

  theParticleGun->SetParticleDefinition(theOpticalPhoton);
  theParticleGun->SetParticleTime(0.0 * CLHEP::ns);
  theParticleGun->SetParticlePosition(G4ThreeVector(-500.0 * CLHEP::cm, 0.0 * CLHEP::cm, 0.0 * CLHEP::cm));
  theParticleGun->SetParticleMomentumDirection(G4ThreeVector(5.0, 3.0, 0.0));
  theParticleGun->SetParticleEnergy(10.0 * CLHEP::keV);
  setOptPhotonPolar(90.0);

  // initialize the random number engine
  theDRand48Engine = new CLHEP::DRand48Engine();
}

LaserBeamsBarrel::~LaserBeamsBarrel() {
  if (theParticleGun != nullptr) {
    delete theParticleGun;
  }
  if (theDRand48Engine != nullptr) {
    delete theDRand48Engine;
  }
}

void LaserBeamsBarrel::GeneratePrimaries(G4Event *myEvent) {
  // this function is called at the beginning of an Event in
  // LaserAlignment::upDate(const BeginOfEvent * myEvent)

  // use the random number generator service of the framework
  edm::Service<edm::RandomNumberGenerator> rng;
  unsigned int seed = rng->mySeed();

  // set the seed
  theDRand48Engine->setSeed(seed);

  // number of LaserBeams
  const G4int nLaserBeams = 8;

  // z position of the Laserdiodes (value from design drawings)
  G4double LaserPositionZ = 1137.0 * CLHEP::mm;

  // Radius of the Laser ring
  G4double LaserRingRadius = 564.0 * CLHEP::mm;

  // phi positions of the Laserdiodes (from CMS Note 2001/053 or from
  // http://abbaneo.home.cern.ch/abbaneo/cms/layout)
  G4double LaserPhi[nLaserBeams] = {G4double(7.0 / 112.0) * G4double(2.0 * M_PI),
                                    G4double(23.0 / 112.0) * G4double(2.0 * M_PI),
                                    G4double(33.0 / 112.0) * G4double(2.0 * M_PI),
                                    G4double(49.0 / 112.0) * G4double(2.0 * M_PI),
                                    G4double(65.0 / 112.0) * G4double(2.0 * M_PI),
                                    G4double(77.0 / 112.0) * G4double(2.0 * M_PI),
                                    G4double(93.0 / 112.0) * G4double(2.0 * M_PI),
                                    G4double(103.0 / 112.0) * G4double(2.0 * M_PI)};

  // width of the LaserBeams
  G4double LaserBeamSigmaX = 0.5 * CLHEP::mm;
  G4double LaserBeamSigmaY = 0.5 * CLHEP::mm;

  // get the definition of the optical photon
  G4ParticleTable *theParticleTable = G4ParticleTable::GetParticleTable();
  G4ParticleDefinition *theOpticalPhoton = theParticleTable->FindParticle("opticalphoton");

  // loop over the LaserBeams
  for (int theBeam = 0; theBeam < nLaserBeams; theBeam++) {
    // code for forward and backward beam
    // calculate x and y position of the current laser diode
    G4double LaserPositionX = cos(LaserPhi[theBeam]) * LaserRingRadius;
    G4double LaserPositionY = sin(LaserPhi[theBeam]) * LaserRingRadius;

    // loop over all the particles in one beam
    for (int theParticle = 0; theParticle < thenParticle; theParticle++) {
      // get randomnumbers  and calculate the position
      CLHEP::RandGaussQ aGaussObjX(*theDRand48Engine, LaserPositionX, LaserBeamSigmaX);
      CLHEP::RandGaussQ aGaussObjY(*theDRand48Engine, LaserPositionY, LaserBeamSigmaY);

      G4double theXPosition = aGaussObjX.fire();
      G4double theYPosition = aGaussObjY.fire();
      G4double theZPosition = LaserPositionZ;

      // set the properties of the newly created particle
      theParticleGun->SetParticleDefinition(theOpticalPhoton);
      theParticleGun->SetParticleTime(0.0 * CLHEP::ns);
      theParticleGun->SetParticlePosition(G4ThreeVector(theXPosition, theYPosition, theZPosition));
      theParticleGun->SetParticleEnergy(thePhotonEnergy);

      // loop over both directions of the beam
      for (int theDirection = 0; theDirection < 2; theDirection++) {
        // shoot in both beam directions ...
        if (theDirection == 0)  // shoot in forward direction (+z)
        {
          theParticleGun->SetParticleMomentumDirection(G4ThreeVector(0.0, 0.0, 1.0));
          // set the polarization
          setOptPhotonPolar(90.0);
          // generate the particle
          theParticleGun->GeneratePrimaryVertex(myEvent);
        } else if (theDirection == 1)  // shoot in backward direction (-z)
        {
          theParticleGun->SetParticleMomentumDirection(G4ThreeVector(0.0, 0.0, -1.0));
          // set the polarization
          setOptPhotonPolar(90.0);
          // generate the particle
          theParticleGun->GeneratePrimaryVertex(myEvent);
        }
      }  // end looop over both beam directions
    }  // end looop over particles in beam
  }  // end loop over beams
}

void LaserBeamsBarrel::setOptPhotonPolar(G4double Angle) {
  /* *********************************************************************** */
  /*   to get optical processes working properly, you have to make sure      *
   *   that the photon polarisation is defined.                              */
  /* *********************************************************************** */

  // first check if we have an optical photon
  if (theParticleGun->GetParticleDefinition()->GetParticleName() != "opticalphoton") {
    edm::LogWarning("SimLaserAlignment:LaserBeamsBarrel")
        << "<LaserBeamsBarrel::setOptPhotonPolar()>: WARNING! The ParticleGun "
           "is not an optical photon";
    return;
  }

  //   G4cout << "  AC1CMS: The ParticleGun is an " <<
  //   theParticleGun->GetParticleDefinition()->GetParticleName();
  G4ThreeVector normal(1.0, 0.0, 0.0);
  G4ThreeVector kphoton = theParticleGun->GetParticleMomentumDirection();
  G4ThreeVector product = normal.cross(kphoton);
  G4double modul2 = product * product;

  G4ThreeVector e_perpendicular(0.0, 0.0, 1.0);

  if (modul2 > 0.0) {
    e_perpendicular = (1.0 / sqrt(modul2)) * product;
  }

  G4ThreeVector e_parallel = e_perpendicular.cross(kphoton);

  G4ThreeVector polar = cos(Angle) * e_parallel + sin(Angle) * e_perpendicular;

  //   G4cout << ", the polarization = " << polar << G4endl;
  theParticleGun->SetParticlePolarization(polar);
}