CSCDigitizer/src/CSCDriftSim.cc

0001 #include "DataFormats/GeometryVector/interface/LocalVector.h"
0002 #include "DataFormats/Math/interface/Point3D.h"
0003 #include "DataFormats/Math/interface/Vector3D.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "Geometry/CSCGeometry/interface/CSCChamberSpecs.h"
0006 #include "Geometry/CSCGeometry/interface/CSCLayer.h"
0007 #include "Geometry/CSCGeometry/interface/CSCLayerGeometry.h"
0008 #include "MagneticField/Engine/interface/MagneticField.h"
0009 #include "Math/GenVector/RotationZ.h"
0010 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0011 #include "SimMuon/CSCDigitizer/src/CSCDetectorHit.h"
0012 #include "SimMuon/CSCDigitizer/src/CSCDriftSim.h"
0013
0014 #include <CLHEP/Random/RandFlat.h>
0015 #include <CLHEP/Random/RandGaussQ.h>
0016 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0017 #include <CLHEP/Units/SystemOfUnits.h>
0018
0019 #include <cmath>
0020 #include <iostream>
0021
0022 static const int N_INTEGRAL_STEPS = 700;
0023
0024 CSCDriftSim::CSCDriftSim()
0025     : bz(0.),  // should make these local variables
0026       ycell(0.),
0027       zcell(0.),
0028       dNdEIntegral(N_INTEGRAL_STEPS, 0.),
0029       STEP_SIZE(0.01),
0030       ELECTRON_DIFFUSION_COEFF(0.0161),
0031       theMagneticField(nullptr) {
0032   // just initialize avalanche sim.  There has to be a better
0033   // way to take the integral of a function!
0034   double sum = 0.;
0035   int i;
0036   for (i = 0; i < N_INTEGRAL_STEPS; ++i) {
0037     if (i > 1) {
0038       double xx = STEP_SIZE * (double(i) - 0.5);
0039       double dNdE = pow(xx, 0.38) * exp(-1.38 * xx);
0040
0041       sum += dNdE;
0042     }
0043     // store this value in the map
0044     dNdEIntegral[i] = sum;
0045   }
0046
0047   // now normalize the whole map
0048   for (i = 0; i < N_INTEGRAL_STEPS; ++i) {
0049     dNdEIntegral[i] /= sum;
0050   }
0051 }
0052
0053 CSCDriftSim::~CSCDriftSim() {}
0054
0055 CSCDetectorHit CSCDriftSim::getWireHit(const Local3DPoint &pos,
0056                                        const CSCLayer *layer,
0057                                        int nearestWire,
0058                                        const PSimHit &simHit,
0059                                        CLHEP::HepRandomEngine *engine) {
0060   const CSCChamberSpecs *specs = layer->chamber()->specs();
0061   const CSCLayerGeometry *geom = layer->geometry();
0062   math::LocalPoint clusterPos(pos.x(), pos.y(), pos.z());
0063   LogTrace("CSCDriftSim") << "CSCDriftSim: ionization cluster at: " << pos;
0064   // set the coordinate system with the x-axis along the nearest wire,
0065   // with the origin in the center of the chamber, on that wire.
0066   math::LocalVector yShift(0, -1. * geom->yOfWire(nearestWire), 0.);
0067   ROOT::Math::RotationZ rotation(-1. * geom->wireAngle());
0068
0069   clusterPos = yShift + clusterPos;
0070   clusterPos = rotation * clusterPos;
0071   GlobalPoint globalPosition = layer->surface().toGlobal(pos);
0072   assert(theMagneticField != nullptr);
0073
0074   //  bz = theMagneticField->inTesla(globalPosition).z() * 10.;
0075
0076   // We need magnetic field in _local_ coordinates
0077   // Interface now allows access in kGauss directly.
0078   bz = layer->toLocal(theMagneticField->inKGauss(globalPosition)).z();
0079
0080   // these subroutines label the coordinates in GEANT coords...
0081   ycell = clusterPos.z() / specs->anodeCathodeSpacing();
0082   zcell = 2. * clusterPos.y() / specs->wireSpacing();
0083
0084   LogTrace("CSCDriftSim") << "CSCDriftSim: bz " << bz << " avgDrift " << avgDrift() << " wireAngle "
0085                           << geom->wireAngle() << " ycell " << ycell << " zcell " << zcell;
0086
0087   double avgPathLength, pathSigma, avgDriftTime, driftTimeSigma;
0088   static const float B_FIELD_CUT = 15.f;
0089   if (fabs(bz) < B_FIELD_CUT) {
0090     avgPathLength = avgPathLengthLowB();
0091     pathSigma = pathSigmaLowB();
0092     avgDriftTime = avgDriftTimeLowB();
0093     driftTimeSigma = driftTimeSigmaLowB();
0094   } else {
0095     avgPathLength = avgPathLengthHighB();
0096     pathSigma = pathSigmaHighB();
0097     avgDriftTime = avgDriftTimeHighB();
0098     driftTimeSigma = driftTimeSigmaHighB();
0099   }
0100
0101   // electron drift path length
0102   double pathLength = std::max(CLHEP::RandGaussQ::shoot(engine, avgPathLength, pathSigma), 0.);
0103
0104   // electron drift distance along the anode wire, including diffusion
0105   double diffusionSigma = ELECTRON_DIFFUSION_COEFF * sqrt(pathLength);
0106   double x = clusterPos.x() + CLHEP::RandGaussQ::shoot(engine, avgDrift(), driftSigma()) +
0107              CLHEP::RandGaussQ::shoot(engine, 0., diffusionSigma);
0108
0109   // electron drift time
0110   double driftTime = std::max(CLHEP::RandGaussQ::shoot(engine, avgDriftTime, driftTimeSigma), 0.);
0111
0112   //@@ Parameters which should be defined outside the code
0113   // f_att is the fraction of drift electrons lost due to attachment
0114   // static const double f_att = 0.5;
0115   static const double f_collected = 0.82;
0116
0117   // Avalanche charge, with fluctuation ('avalancheCharge()' is the fluctuation
0118   // generator!)
0119   // double charge = avalancheCharge() * f_att * f_collected *
0120   // gasGain(layer->id()) * e_SI * 1.e15;
0121   // doing fattachment by random chance of killing
0122   double charge = avalancheCharge(engine) * f_collected * gasGain(layer->id()) * CLHEP::e_SI * 1.e15;
0123
0124   float t = simHit.tof() + driftTime;
0125   LogTrace("CSCDriftSim") << "CSCDriftSim: tof = " << simHit.tof() << " driftTime = " << driftTime
0126                           << " MEDH = " << CSCDetectorHit(nearestWire, charge, x, t, &simHit);
0127   return CSCDetectorHit(nearestWire, charge, x, t, &simHit);
0128 }
0129
0130 // Generate avalanche fluctuation
0131 #include <algorithm>
0132 double CSCDriftSim::avalancheCharge(CLHEP::HepRandomEngine *engine) {
0133   double returnVal = 0.;
0134   // pick a random value along the dNdE integral
0135   double x = CLHEP::RandFlat::shoot(engine);
0136   size_t i;
0137   size_t isiz = dNdEIntegral.size();
0138   /*
0139   for(i = 0; i < isiz-1; ++i) {
0140     if(dNdEIntegral[i] > x) break;
0141   }
0142   */
0143   // return position of first element with a value >= x
0144   std::vector<double>::const_iterator p = lower_bound(dNdEIntegral.begin(), dNdEIntegral.end(), x);
0145   if (p == dNdEIntegral.end())
0146     i = isiz - 1;
0147   else
0148     i = p - dNdEIntegral.begin();
0149
0150   // now extrapolate between values
0151   if (i == isiz - 1) {
0152     // edm::LogInfo("CSCDriftSim") << "Funky integral in CSCDriftSim " << x;
0153     returnVal = STEP_SIZE * double(i) * dNdEIntegral[i];
0154   } else {
0155     double x1 = dNdEIntegral[i];
0156     double x2 = dNdEIntegral[i + 1];
0157     returnVal = STEP_SIZE * (double(i) + (x - x1) / (x2 - x1));
0158   }
0159   LogTrace("CSCDriftSim") << "CSCDriftSim: avalanche fluc " << returnVal << "  " << x;
0160
0161   return returnVal;
0162 }
0163
0164 double CSCDriftSim::gasGain(const CSCDetId &detId) const {
0165   double result = 130000.;  // 1.30 E05
0166   // if ME1/1, add some extra gas gain to compensate
0167   // for a smaller gas gap
0168   int ring = detId.ring();
0169   if (detId.station() == 1 && (ring == 1 || ring == 4)) {
0170     result = 213000.;  // 2.13 E05 ( to match real world as of Jan-2011)
0171   }
0172   return result;
0173 }