PPSPixelDigiProducer/src/RPixLinearChargeDivider.cc

0001 #include "SimPPS/PPSPixelDigiProducer/interface/RPixLinearChargeDivider.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0004 #include "DataFormats/GeometryVector/interface/LocalVector.h"
0005
0006 RPixLinearChargeDivider::RPixLinearChargeDivider(const edm::ParameterSet& params,
0007                                                  CLHEP::HepRandomEngine& eng,
0008                                                  uint32_t det_id)
0009     : rndEngine_(eng), det_id_(det_id) {
0010   verbosity_ = params.getParameter<int>("RPixVerbosity");
0011   fluctuate = new SiG4UniversalFluctuation();
0012   fluctuateCharge_ = params.getParameter<bool>("RPixLandauFluctuations");
0013   chargedivisions_ = params.getParameter<int>("RPixChargeDivisions");
0014   deltaCut_ = params.getParameter<double>("RPixDeltaProductionCut");
0015 }
0016
0017 RPixLinearChargeDivider::~RPixLinearChargeDivider() { delete fluctuate; }
0018
0019 std::vector<RPixEnergyDepositUnit> RPixLinearChargeDivider::divide(const PSimHit& hit) {
0020   LocalVector direction = hit.exitPoint() - hit.entryPoint();
0021   if (direction.z() > 10 || direction.x() > 200 || direction.y() > 200) {
0022     the_energy_path_distribution_.clear();
0023     return the_energy_path_distribution_;
0024   }
0025
0026   int NumberOfSegmentation = chargedivisions_;
0027   double eLoss = hit.energyLoss();  // Eloss in GeV
0028   the_energy_path_distribution_.resize(NumberOfSegmentation);
0029
0030   if (fluctuateCharge_) {
0031     int pid = hit.particleType();
0032     double momentum = hit.pabs();
0033     double length = direction.mag();  // Track length in Silicon
0034     FluctuateEloss(pid, momentum, eLoss, length, NumberOfSegmentation, the_energy_path_distribution_);
0035     for (int i = 0; i < NumberOfSegmentation; i++) {
0036       the_energy_path_distribution_[i].setPosition(hit.entryPoint() +
0037                                                    double((i + 0.5) / NumberOfSegmentation) * direction);
0038     }
0039   } else {
0040     for (int i = 0; i < NumberOfSegmentation; i++) {
0041       the_energy_path_distribution_[i].setPosition(hit.entryPoint() +
0042                                                    double((i + 0.5) / NumberOfSegmentation) * direction);
0043       the_energy_path_distribution_[i].setEnergy(eLoss / (double)NumberOfSegmentation);
0044     }
0045   }
0046
0047   if (verbosity_) {
0048     edm::LogInfo("PPS") << "RPixLinearChargeDivider " << det_id_ << " charge along the track:";
0049     double sum = 0;
0050     for (unsigned int i = 0; i < the_energy_path_distribution_.size(); i++) {
0051       edm::LogInfo("RPixLinearChargeDivider")
0052           << the_energy_path_distribution_[i].Position().x() << " " << the_energy_path_distribution_[i].Position().y()
0053           << " " << the_energy_path_distribution_[i].Position().z() << " " << the_energy_path_distribution_[i].Energy();
0054       sum += the_energy_path_distribution_[i].Energy();
0055     }
0056     edm::LogInfo("PPS") << "RPixLinearChargeDivider "
0057                         << "energy dep. sum=" << sum;
0058   }
0059
0060   return the_energy_path_distribution_;
0061 }
0062
0063 void RPixLinearChargeDivider::FluctuateEloss(int pid,
0064                                              double particleMomentum,
0065                                              double eloss,
0066                                              double length,
0067                                              int NumberOfSegs,
0068                                              std::vector<RPixEnergyDepositUnit>& elossVector) {
0069   double particleMass = 139.6;  // Mass in MeV, Assume pion
0070   pid = std::abs(pid);
0071   if (pid != 211) {  // Mass in MeV
0072     if (pid == 11)
0073       particleMass = 0.511;
0074     else if (pid == 13)
0075       particleMass = 105.7;
0076     else if (pid == 321)
0077       particleMass = 493.7;
0078     else if (pid == 2212)
0079       particleMass = 938.3;
0080   }
0081
0082   double segmentLength = length / NumberOfSegs;
0083
0084   // Generate charge fluctuations.
0085   double de = 0.;
0086   double sum = 0.;
0087   double segmentEloss = (eloss * 1000) / NumberOfSegs;  //eloss in MeV
0088   for (int i = 0; i < NumberOfSegs; i++) {
0089     double deltaCutoff = deltaCut_;
0090     de = fluctuate->SampleFluctuations(
0091              particleMomentum * 1000, particleMass, deltaCutoff, segmentLength, segmentEloss, &(rndEngine_)) /
0092          1000;  //convert to GeV
0093     elossVector[i].setEnergy(de);
0094     sum += de;
0095   }
0096
0097   if (sum > 0.) {  // If fluctuations give eloss>0.
0098                    // Rescale to the same total eloss
0099     double ratio = eloss / sum;
0100     for (int ii = 0; ii < NumberOfSegs; ii++)
0101       elossVector[ii].setEnergy(ratio * elossVector[ii].Energy());
0102   } else {  // If fluctuations gives 0 eloss
0103     double averageEloss = eloss / NumberOfSegs;
0104     for (int ii = 0; ii < NumberOfSegs; ii++)
0105       elossVector[ii].setEnergy(averageEloss);
0106   }
0107   return;
0108 }