RPDigiProducer/plugins/RPLinearChargeDivider.cc

0001 #include "SimPPS/RPDigiProducer/plugins/RPLinearChargeDivider.h"
0002 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0003 #include "DataFormats/GeometryVector/interface/LocalVector.h"
0004 #include "Geometry/VeryForwardRPTopology/interface/RPTopology.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006
0007 RPLinearChargeDivider::RPLinearChargeDivider(const edm::ParameterSet& params,
0008                                              CLHEP::HepRandomEngine& eng,
0009                                              RPDetId det_id)
0010     : rndEngine_(eng), det_id_(det_id) {
0011   verbosity_ = params.getParameter<int>("RPVerbosity");
0012
0013   fluctuate_ = std::make_unique<SiG4UniversalFluctuation>();
0014
0015   // To Run APV in peak instead of deconvolution mode, which degrades the time resolution.
0016   //use: SimpleConfigurable<bool> SiLinearChargeDivider::peakMode(false,"SiStripDigitizer:APVpeakmode");
0017
0018   // To Enable interstrip Landau fluctuations within a cluster.
0019   //use: SimpleConfigurable<bool> SiLinearChargeDivider::fluctuateCharge(true,"SiStripDigitizer:LandauFluctuations");
0020   fluctuateCharge_ = params.getParameter<bool>("RPLandauFluctuations");
0021
0022   // Number of segments per strip into which charge is divided during
0023   // simulation. If large, precision of simulation improves.
0024   //to do so: SimpleConfigurable<int> SiLinearChargeDivider::chargeDivisionsPerStrip(10,"SiStripDigitizer:chargeDivisionsPerStrip");
0025   chargedivisionsPerStrip_ = params.getParameter<int>("RPChargeDivisionsPerStrip");
0026   chargedivisionsPerThickness_ = params.getParameter<int>("RPChargeDivisionsPerThickness");
0027
0028   // delta cutoff in MeV, has to be same as in OSCAR (0.120425 MeV corresponding // to 100um range for electrons)
0029   //        SimpleConfigurable<double>  SiLinearChargeDivider::deltaCut(0.120425,
0030   deltaCut_ = params.getParameter<double>("RPDeltaProductionCut");
0031
0032   RPTopology rp_det_topol;
0033   pitch_ = rp_det_topol.DetPitch();
0034   thickness_ = rp_det_topol.DetThickness();
0035 }
0036
0037 RPLinearChargeDivider::~RPLinearChargeDivider() {}
0038
0039 simromanpot::energy_path_distribution RPLinearChargeDivider::divide(const PSimHit& hit) {
0040   LocalVector direction = hit.exitPoint() - hit.entryPoint();
0041   if (direction.z() > 10 || direction.x() > 200 || direction.y() > 200) {
0042     the_energy_path_distribution_.clear();
0043     return the_energy_path_distribution_;
0044   }
0045
0046   int NumberOfSegmentation_y = (int)(1 + chargedivisionsPerStrip_ * fabs(direction.y()) / pitch_);
0047   int NumberOfSegmentation_z = (int)(1 + chargedivisionsPerThickness_ * fabs(direction.z()) / thickness_);
0048   int NumberOfSegmentation = std::max(NumberOfSegmentation_y, NumberOfSegmentation_z);
0049
0050   double eLoss = hit.energyLoss();  // Eloss in GeV
0051
0052   the_energy_path_distribution_.resize(NumberOfSegmentation);
0053
0054   if (fluctuateCharge_) {
0055     int pid = hit.particleType();
0056     double momentum = hit.pabs();
0057     double length = direction.mag();  // Track length in Silicon
0058     FluctuateEloss(pid, momentum, eLoss, length, NumberOfSegmentation, the_energy_path_distribution_);
0059     for (int i = 0; i < NumberOfSegmentation; i++) {
0060       the_energy_path_distribution_[i].setPosition(hit.entryPoint() +
0061                                                    double((i + 0.5) / NumberOfSegmentation) * direction);
0062     }
0063   } else {
0064     for (int i = 0; i < NumberOfSegmentation; i++) {
0065       the_energy_path_distribution_[i].setPosition(hit.entryPoint() +
0066                                                    double((i + 0.5) / NumberOfSegmentation) * direction);
0067       the_energy_path_distribution_[i].setEnergy(eLoss / (double)NumberOfSegmentation);
0068     }
0069   }
0070
0071   if (verbosity_) {
0072     edm::LogInfo("RPLinearChargeDivider") << det_id_ << " charge along the track:\n";
0073     double sum = 0;
0074     for (unsigned int i = 0; i < the_energy_path_distribution_.size(); i++) {
0075       edm::LogInfo("RPLinearChargeDivider")
0076           << the_energy_path_distribution_[i].Position().x() << " " << the_energy_path_distribution_[i].Position().y()
0077           << " " << the_energy_path_distribution_[i].Position().z() << " " << the_energy_path_distribution_[i].Energy()
0078           << "\n";
0079       sum += the_energy_path_distribution_[i].Energy();
0080     }
0081     edm::LogInfo("RPLinearChargeDivider") << "energy dep. sum=" << sum << "\n";
0082   }
0083
0084   return the_energy_path_distribution_;
0085 }
0086
0087 void RPLinearChargeDivider::FluctuateEloss(int pid,
0088                                            double particleMomentum,
0089                                            double eloss,
0090                                            double length,
0091                                            int NumberOfSegs,
0092                                            simromanpot::energy_path_distribution& elossVector) {
0093   double particleMass = 139.6;  // Mass in MeV, Assume pion
0094   pid = std::abs(pid);
0095   if (pid != 211) {  // Mass in MeV
0096     if (pid == 11)
0097       particleMass = 0.511;
0098     else if (pid == 13)
0099       particleMass = 105.7;
0100     else if (pid == 321)
0101       particleMass = 493.7;
0102     else if (pid == 2212)
0103       particleMass = 938.3;
0104   }
0105
0106   double segmentLength = length / NumberOfSegs;
0107
0108   // Generate charge fluctuations.
0109   double de = 0.;
0110   double sum = 0.;
0111   double segmentEloss = (eloss * 1000) / NumberOfSegs;  //eloss in MeV
0112   for (int i = 0; i < NumberOfSegs; i++) {
0113     // The G4 routine needs momentum in MeV, mass in Mev, delta-cut in MeV,
0114     // track segment length in mm, segment eloss in MeV
0115     // Returns fluctuated eloss in MeV
0116     double deltaCutoff = deltaCut_;
0117     de = fluctuate_->SampleFluctuations(
0118              particleMomentum * 1000, particleMass, deltaCutoff, segmentLength, segmentEloss, &(rndEngine_)) /
0119          1000;  //convert to GeV
0120     elossVector[i].setEnergy(de);
0121     sum += de;
0122   }
0123
0124   if (sum > 0.) {  // If fluctuations give eloss>0.
0125     // Rescale to the same total eloss
0126     double ratio = eloss / sum;
0127     for (int ii = 0; ii < NumberOfSegs; ii++)
0128       elossVector[ii].setEnergy(ratio * elossVector[ii].Energy());
0129   } else {  // If fluctuations gives 0 eloss
0130     double averageEloss = eloss / NumberOfSegs;
0131     for (int ii = 0; ii < NumberOfSegs; ii++)
0132       elossVector[ii].setEnergy(averageEloss);
0133   }
0134   return;
0135 }