#include <RPLinearChargeDivider.h>

Public Member Functions
simromanpot::energy_path_distribution	divide (const PSimHit &hit)

	RPLinearChargeDivider (const edm::ParameterSet &params, CLHEP::HepRandomEngine &eng, RPDetId det_id)

	~RPLinearChargeDivider ()

Private Member Functions
void	FluctuateEloss (int pid, double particleMomentum, double eloss, double length, int NumberOfSegs, simromanpot::energy_path_distribution &elossVector)

Private Attributes
int	chargedivisionsPerStrip_

int	chargedivisionsPerThickness_

double	deltaCut_

RPDetId	det_id_

std::unique_ptr< SiG4UniversalFluctuation >	fluctuate_

bool	fluctuateCharge_

const edm::ParameterSet &	params_

double	pitch_

CLHEP::HepRandomEngine &	rndEngine_

simromanpot::energy_path_distribution	the_energy_path_distribution_

double	thickness_

int	verbosity_

Detailed Description

Definition at line 13 of file RPLinearChargeDivider.h.

Constructor & Destructor Documentation

RPLinearChargeDivider::RPLinearChargeDivider	(	const edm::ParameterSet &	params,
		CLHEP::HepRandomEngine &	eng,
		RPDetId	det_id
	)

Definition at line 7 of file RPLinearChargeDivider.cc.

References chargedivisionsPerStrip_, chargedivisionsPerThickness_, deltaCut_, fluctuate_, fluctuateCharge_, edm::ParameterSet::getParameter(), pitch_, thickness_, and verbosity_.

     : params_(params), rndEngine_(eng), det_id_(det_id) {
   verbosity_ = params.getParameter<int>("RPVerbosity");
 
   fluctuate_ = std::make_unique<SiG4UniversalFluctuation>();
 
   // To Run APV in peak instead of deconvolution mode, which degrades the time resolution.
   //use: SimpleConfigurable<bool> SiLinearChargeDivider::peakMode(false,"SiStripDigitizer:APVpeakmode");
 
   // To Enable interstrip Landau fluctuations within a cluster.
   //use: SimpleConfigurable<bool> SiLinearChargeDivider::fluctuateCharge(true,"SiStripDigitizer:LandauFluctuations");
   fluctuateCharge_ = params.getParameter<bool>("RPLandauFluctuations");
 
   // Number of segments per strip into which charge is divided during
   // simulation. If large, precision of simulation improves.
   //to do so: SimpleConfigurable<int> SiLinearChargeDivider::chargeDivisionsPerStrip(10,"SiStripDigitizer:chargeDivisionsPerStrip");
   chargedivisionsPerStrip_ = params.getParameter<int>("RPChargeDivisionsPerStrip");
   chargedivisionsPerThickness_ = params.getParameter<int>("RPChargeDivisionsPerThickness");
 
   // delta cutoff in MeV, has to be same as in OSCAR (0.120425 MeV corresponding // to 100um range for electrons)
   //        SimpleConfigurable<double>  SiLinearChargeDivider::deltaCut(0.120425,
   deltaCut_ = params.getParameter<double>("RPDeltaProductionCut");
 
   RPTopology rp_det_topol;
   pitch_ = rp_det_topol.DetPitch();
   thickness_ = rp_det_topol.DetThickness();
 }

RPLinearChargeDivider::~RPLinearChargeDivider ( )

Definition at line 37 of file RPLinearChargeDivider.cc.

37 {}

Member Function Documentation

simromanpot::energy_path_distribution RPLinearChargeDivider::divide ( const PSimHit & hit )

Definition at line 39 of file RPLinearChargeDivider.cc.

References chargedivisionsPerStrip_, chargedivisionsPerThickness_, det_id_, simKBmtfDigis_cfi::eLoss, PSimHit::energyLoss(), PSimHit::entryPoint(), PSimHit::exitPoint(), fluctuateCharge_, FluctuateEloss(), mps_fire::i, createfilelist::int, PV3DBase< T, PVType, FrameType >::mag(), SiStripPI::max, PSimHit::pabs(), PSimHit::particleType(), pitch_, the_energy_path_distribution_, thickness_, verbosity_, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

                                                                                     {
   LocalVector direction = hit.exitPoint() - hit.entryPoint();
   if (direction.z() > 10 || direction.x() > 200 || direction.y() > 200) {
     the_energy_path_distribution_.clear();
     return the_energy_path_distribution_;
   }
 
   int NumberOfSegmentation_y = (int)(1 + chargedivisionsPerStrip_ * fabs(direction.y()) / pitch_);
   int NumberOfSegmentation_z = (int)(1 + chargedivisionsPerThickness_ * fabs(direction.z()) / thickness_);
   int NumberOfSegmentation = std::max(NumberOfSegmentation_y, NumberOfSegmentation_z);
 
   double eLoss = hit.energyLoss();  // Eloss in GeV
 
   the_energy_path_distribution_.resize(NumberOfSegmentation);
 
   if (fluctuateCharge_) {
     int pid = hit.particleType();
     double momentum = hit.pabs();
     double length = direction.mag();  // Track length in Silicon
     FluctuateEloss(pid, momentum, eLoss, length, NumberOfSegmentation, the_energy_path_distribution_);
     for (int i = 0; i < NumberOfSegmentation; i++) {
       the_energy_path_distribution_[i].setPosition(hit.entryPoint() +
                                                    double((i + 0.5) / NumberOfSegmentation) * direction);
     }
   } else {
     for (int i = 0; i < NumberOfSegmentation; i++) {
       the_energy_path_distribution_[i].setPosition(hit.entryPoint() +
                                                    double((i + 0.5) / NumberOfSegmentation) * direction);
       the_energy_path_distribution_[i].setEnergy(eLoss / (double)NumberOfSegmentation);
     }
   }
 
   if (verbosity_) {
     edm::LogInfo("RPLinearChargeDivider") << det_id_ << " charge along the track:\n";
     double sum = 0;
     for (unsigned int i = 0; i < the_energy_path_distribution_.size(); i++) {
       edm::LogInfo("RPLinearChargeDivider")
           << the_energy_path_distribution_[i].Position().x() << " " << the_energy_path_distribution_[i].Position().y()
           << " " << the_energy_path_distribution_[i].Position().z() << " " << the_energy_path_distribution_[i].Energy()
           << "\n";
       sum += the_energy_path_distribution_[i].Energy();
     }
     edm::LogInfo("RPLinearChargeDivider") << "energy dep. sum=" << sum << "\n";
   }
 
   return the_energy_path_distribution_;
 }

void RPLinearChargeDivider::FluctuateEloss	(	int	pid,
		double	particleMomentum,
		double	eloss,
		double	length,
		int	NumberOfSegs,
		simromanpot::energy_path_distribution &	elossVector
	)

private

Definition at line 87 of file RPLinearChargeDivider.cc.

References funct::abs(), deltaCut_, EcalCondDBWriter_cfi::Energy, fluctuate_, mps_fire::i, cuy::ii, particleFlowDisplacedVertex_cfi::ratio, and rndEngine_.

Referenced by divide().

                                                                                              {
   double particleMass = 139.6;  // Mass in MeV, Assume pion
   pid = std::abs(pid);
   if (pid != 211) {  // Mass in MeV
     if (pid == 11)
       particleMass = 0.511;
     else if (pid == 13)
       particleMass = 105.7;
     else if (pid == 321)
       particleMass = 493.7;
     else if (pid == 2212)
       particleMass = 938.3;
   }
 
   double segmentLength = length / NumberOfSegs;
 
   // Generate charge fluctuations.
   double de = 0.;
   double sum = 0.;
   double segmentEloss = (eloss * 1000) / NumberOfSegs;  //eloss in MeV
   for (int i = 0; i < NumberOfSegs; i++) {
     // The G4 routine needs momentum in MeV, mass in Mev, delta-cut in MeV,
     // track segment length in mm, segment eloss in MeV
     // Returns fluctuated eloss in MeV
     double deltaCutoff = deltaCut_;
     de = fluctuate_->SampleFluctuations(
              particleMomentum * 1000, particleMass, deltaCutoff, segmentLength, segmentEloss, &(rndEngine_)) /
          1000;  //convert to GeV
     elossVector[i].setEnergy(de);
     sum += de;
   }
 
   if (sum > 0.) {  // If fluctuations give eloss>0.
     // Rescale to the same total eloss
     double ratio = eloss / sum;
     for (int ii = 0; ii < NumberOfSegs; ii++)
       elossVector[ii].setEnergy(ratio * elossVector[ii].Energy());
   } else {  // If fluctuations gives 0 eloss
     double averageEloss = eloss / NumberOfSegs;
     for (int ii = 0; ii < NumberOfSegs; ii++)
       elossVector[ii].setEnergy(averageEloss);
   }
   return;
 }