SiStripDigitizer/plugins/SiLinearChargeDivider.cc

0001 #include "SiLinearChargeDivider.h"
0002 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0003 #include "DataFormats/GeometryVector/interface/LocalVector.h"
0004 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include <fstream>
0007 #include <string>
0008
0009 SiLinearChargeDivider::SiLinearChargeDivider(const edm::ParameterSet& conf)
0010     :  // Run APV in peak instead of deconvolution mode, which degrades the time resolution.
0011       peakMode(conf.getParameter<bool>("APVpeakmode")),
0012       // Enable interstrip Landau fluctuations within a cluster.
0013       fluctuateCharge(conf.getParameter<bool>("LandauFluctuations")),
0014       // Number of segments per strip into which charge is divided during
0015       // simulation. If large, precision of simulation improves.
0016       chargedivisionsPerStrip(conf.getParameter<int>("chargeDivisionsPerStrip")),
0017       // delta cutoff in MeV, has to be same as in Geant (0.120425 MeV corresponding to 100um range for electrons)
0018       deltaCut(conf.getParameter<double>("DeltaProductionCut")),
0019       //Offset for digitization during the MTCC and in general for taking cosmic particle
0020       //The value to be used it must be evaluated and depend on the volume defnition used
0021       //for the cosimc generation (Considering only the tracker the value is 11 ns)
0022       cosmicShift(conf.getUntrackedParameter<double>("CosmicDelayShift")),
0023       theParticleDataTable(nullptr),
0024       // Geant4 engine used to fluctuate the charge from segment to segment
0025       fluctuate(new SiG4UniversalFluctuation())
0026
0027 {
0028   readPulseShape(conf.getParameter<edm::FileInPath>(peakMode ? "APVShapePeakFile" : "APVShapeDecoFile").fullPath());
0029 }
0030
0031 SiLinearChargeDivider::~SiLinearChargeDivider() {}
0032
0033 void SiLinearChargeDivider::readPulseShape(const std::string& pulseShapeFileName) {
0034   // Pulse shape file format: empty lines and comments (lines starting with '#') are ignored
0035   // one line "resolution: value" is interpreted as the resolution
0036   // all other lines are read as consecutive values of the shape
0037   std::ifstream shapeFile(pulseShapeFileName.c_str());
0038   if (!shapeFile.good()) {
0039     throw cms::Exception("FileError") << "Problem opening APV Shape file: " << pulseShapeFileName;
0040   }
0041   pulseResolution = -1.;
0042   std::string line;
0043   const std::string resoPrefix{"resolution: "};
0044   while (std::getline(shapeFile, line)) {
0045     if ((!line.empty()) && (line.substr(1) != "#")) {
0046       std::istringstream lStr{line};
0047       if (line.substr(0, resoPrefix.size()) == resoPrefix) {
0048         lStr.seekg(resoPrefix.size());
0049         lStr >> pulseResolution;
0050       } else {
0051         double value;
0052         while (lStr >> value) {
0053           pulseValues.push_back(value);
0054         }
0055       }
0056     }
0057   }
0058   if (pulseValues.empty() || (pulseResolution == -1.)) {
0059     throw cms::Exception("WrongAPVPulseShape") << "Problem reading from APV pulse shape file " << pulseShapeFileName
0060                                                << ": " << (pulseValues.empty() ? "no values" : "no resolution");
0061   }
0062   const auto maxIt = std::max_element(pulseValues.begin(), pulseValues.end());
0063   if (std::abs((*maxIt) - 1.) > std::numeric_limits<double>::epsilon()) {
0064     throw cms::Exception("WrongAPVPulseShape")
0065         << "The max value of the APV pulse shape stored in the text file used in "
0066            "SimGeneral/MixingModule/python/SiStripSimParameters_cfi.py is not equal to 1. Need to be fixed.";
0067   }
0068   pulset0Idx = std::distance(pulseValues.begin(), maxIt);
0069 }
0070
0071 SiChargeDivider::ionization_type SiLinearChargeDivider::divide(const PSimHit* hit,
0072                                                                const LocalVector& driftdir,
0073                                                                double moduleThickness,
0074                                                                const StripGeomDetUnit& det,
0075                                                                CLHEP::HepRandomEngine* engine) {
0076   // signal after pulse shape correction
0077   float const decSignal = TimeResponse(hit, det);
0078
0079   // if out of time go home!
0080   if (0 == decSignal)
0081     return ionization_type();
0082
0083   // Get the nass if the particle, in MeV.
0084   // Protect from particles with Mass = 0, assuming then the pion mass
0085   assert(theParticleDataTable != nullptr);
0086   ParticleData const* particle = theParticleDataTable->particle(hit->particleType());
0087   double const particleMass = particle ? particle->mass() * 1000 : 139.57;
0088   double const particleCharge = particle ? particle->charge() : 1.;
0089
0090   if (!particle) {
0091     LogDebug("SiLinearChargeDivider") << "Cannot find particle of type " << hit->particleType()
0092                                       << " in the PDT we assign to this particle the mass and charge of the Pion";
0093   }
0094
0095   int NumberOfSegmentation =
0096       // if neutral: just one deposit....
0097       (fabs(particleMass) < 1.e-6 || particleCharge == 0)
0098           ? 1
0099           :
0100           // computes the number of segments from number of segments per strip times number of strips.
0101           (int)(1 + chargedivisionsPerStrip *
0102                         fabs(driftXPos(hit->exitPoint(), driftdir, moduleThickness) -
0103                              driftXPos(hit->entryPoint(), driftdir, moduleThickness)) /
0104                         det.specificTopology().localPitch(hit->localPosition()));
0105
0106   // Eloss in GeV
0107   float eLoss = hit->energyLoss();
0108
0109   // Prepare output
0110   ionization_type _ionization_points;
0111   _ionization_points.resize(NumberOfSegmentation);
0112
0113   // Fluctuate charge in track subsegments
0114   LocalVector direction = hit->exitPoint() - hit->entryPoint();
0115   if (NumberOfSegmentation <= 1) {
0116     // here I need a random... not 0.5
0117     _ionization_points[0] = EnergyDepositUnit(eLoss * decSignal / eLoss, hit->entryPoint() + 0.5f * direction);
0118   } else {
0119     float eLossVector[NumberOfSegmentation];
0120     if (fluctuateCharge) {
0121       fluctuateEloss(particleMass, hit->pabs(), eLoss, direction.mag(), NumberOfSegmentation, eLossVector, engine);
0122       // Save the energy of each segment
0123       for (int i = 0; i != NumberOfSegmentation; i++) {
0124         // take energy value from vector eLossVector,
0125         _ionization_points[i] =
0126             EnergyDepositUnit(eLossVector[i] * decSignal / eLoss,
0127                               hit->entryPoint() + float((i + 0.5) / NumberOfSegmentation) * direction);
0128       }
0129     } else {
0130       // Save the energy of each segment
0131       for (int i = 0; i != NumberOfSegmentation; i++) {
0132         // take energy value from eLoss average over n.segments.
0133         _ionization_points[i] =
0134             EnergyDepositUnit(decSignal / float(NumberOfSegmentation),
0135                               hit->entryPoint() + float((i + 0.5) / NumberOfSegmentation) * direction);
0136       }
0137     }
0138   }
0139   return _ionization_points;
0140 }
0141
0142 void SiLinearChargeDivider::fluctuateEloss(double particleMass,
0143                                            float particleMomentum,
0144                                            float eloss,
0145                                            float length,
0146                                            int NumberOfSegs,
0147                                            float elossVector[],
0148                                            CLHEP::HepRandomEngine* engine) {
0149   // Generate charge fluctuations.
0150   float sum = 0.;
0151   double deltaCutoff;
0152   double mom = particleMomentum * 1000.;
0153   double seglen = length / NumberOfSegs * 10.;
0154   double segeloss = (1000. * eloss) / NumberOfSegs;
0155   for (int i = 0; i < NumberOfSegs; i++) {
0156     // The G4 routine needs momentum in MeV, mass in MeV, delta-cut in MeV,
0157     // track segment length in mm, segment eloss in MeV
0158     // Returns fluctuated eloss in MeV
0159     // the cutoff is sometimes redefined inside, so fix it.
0160     deltaCutoff = deltaCut;
0161     sum += (elossVector[i] =
0162                 fluctuate->SampleFluctuations(mom, particleMass, deltaCutoff, seglen, segeloss, engine) / 1000.);
0163   }
0164
0165   if (sum > 0.) {  // If fluctuations give eloss>0.
0166     // Rescale to the same total eloss
0167     float ratio = eloss / sum;
0168     for (int ii = 0; ii < NumberOfSegs; ii++)
0169       elossVector[ii] = ratio * elossVector[ii];
0170   } else {  // If fluctuations gives 0 eloss
0171     float averageEloss = eloss / NumberOfSegs;
0172     for (int ii = 0; ii < NumberOfSegs; ii++)
0173       elossVector[ii] = averageEloss;
0174   }
0175   return;
0176 }
0177
0178 float SiLinearChargeDivider::TimeResponse(const PSimHit* hit, const StripGeomDetUnit& det) {
0179   // x is difference between the tof and the tof for a photon (reference)
0180   // converted into a bin number
0181   const auto dTOF = det.surface().toGlobal(hit->localPosition()).mag() / 30. + cosmicShift - hit->tof();
0182   const int x = int(dTOF / pulseResolution) + pulset0Idx;
0183   if (x < 0 || x >= int(pulseValues.size()))
0184     return 0;
0185   return hit->energyLoss() * pulseValues[std::size_t(x)];
0186 }