CMS 3D CMS Logo

List of all members | Public Member Functions | Private Member Functions | Private Attributes
SiLinearChargeDivider Class Reference

#include <SiLinearChargeDivider.h>

Inheritance diagram for SiLinearChargeDivider:
SiChargeDivider

Public Member Functions

SiChargeDivider::ionization_type divide (const PSimHit *, const LocalVector &, double, const StripGeomDetUnit &det, CLHEP::HepRandomEngine *) override
 
void setParticleDataTable (const ParticleDataTable *pdt) override
 
 SiLinearChargeDivider (const edm::ParameterSet &conf)
 
 ~SiLinearChargeDivider () override
 
- Public Member Functions inherited from SiChargeDivider
virtual ~SiChargeDivider ()
 

Private Member Functions

float driftXPos (const Local3DPoint &pos, const LocalVector &drift, double thickness)
 
void fluctuateEloss (double const particleMass, float momentum, float eloss, float length, int NumberOfSegmentation, float elossVector[], CLHEP::HepRandomEngine *)
 
void readPulseShape (const std::string &pulseShapeFileName)
 
float TimeResponse (const PSimHit *hit, const StripGeomDetUnit &det)
 

Private Attributes

const int chargedivisionsPerStrip
 
const double cosmicShift
 
const double deltaCut
 
std::unique_ptr< SiG4UniversalFluctuationfluctuate
 
const bool fluctuateCharge
 
const bool peakMode
 
double pulseResolution
 
unsigned int pulset0Idx
 
std::vector< double > pulseValues
 
const ParticleDataTabletheParticleDataTable
 

Additional Inherited Members

- Public Types inherited from SiChargeDivider
typedef std::vector< EnergyDepositUnitionization_type
 

Detailed Description

Concrete implementation of SiChargeDivider. It divides the charge on the line connecting entry and exit point of the SimTrack in the Silicon. Effects that are considered here are:

Definition at line 27 of file SiLinearChargeDivider.h.

Constructor & Destructor Documentation

◆ SiLinearChargeDivider()

SiLinearChargeDivider::SiLinearChargeDivider ( const edm::ParameterSet conf)

Definition at line 9 of file SiLinearChargeDivider.cc.

References edm::FileInPath::fullPath(), edm::ParameterSet::getParameter(), peakMode, and readPulseShape().

10  : // Run APV in peak instead of deconvolution mode, which degrades the time resolution.
11  peakMode(conf.getParameter<bool>("APVpeakmode")),
12  // Enable interstrip Landau fluctuations within a cluster.
13  fluctuateCharge(conf.getParameter<bool>("LandauFluctuations")),
14  // Number of segments per strip into which charge is divided during
15  // simulation. If large, precision of simulation improves.
16  chargedivisionsPerStrip(conf.getParameter<int>("chargeDivisionsPerStrip")),
17  // delta cutoff in MeV, has to be same as in Geant (0.120425 MeV corresponding to 100um range for electrons)
18  deltaCut(conf.getParameter<double>("DeltaProductionCut")),
19  //Offset for digitization during the MTCC and in general for taking cosmic particle
20  //The value to be used it must be evaluated and depend on the volume defnition used
21  //for the cosimc generation (Considering only the tracker the value is 11 ns)
22  cosmicShift(conf.getUntrackedParameter<double>("CosmicDelayShift")),
23  theParticleDataTable(nullptr),
24  // Geant4 engine used to fluctuate the charge from segment to segment
26 
27 {
28  readPulseShape(conf.getParameter<edm::FileInPath>(peakMode ? "APVShapePeakFile" : "APVShapeDecoFile").fullPath());
29 }
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
const ParticleDataTable * theParticleDataTable
std::string fullPath() const
Definition: FileInPath.cc:161
T getUntrackedParameter(std::string const &, T const &) const
std::unique_ptr< SiG4UniversalFluctuation > fluctuate
void readPulseShape(const std::string &pulseShapeFileName)

◆ ~SiLinearChargeDivider()

SiLinearChargeDivider::~SiLinearChargeDivider ( )
override

Definition at line 31 of file SiLinearChargeDivider.cc.

31 {}

Member Function Documentation

◆ divide()

SiChargeDivider::ionization_type SiLinearChargeDivider::divide ( const PSimHit hit,
const LocalVector driftdir,
double  moduleThickness,
const StripGeomDetUnit det,
CLHEP::HepRandomEngine *  engine 
)
overridevirtual

Implements SiChargeDivider.

Definition at line 71 of file SiLinearChargeDivider.cc.

References cms::cuda::assert(), chargedivisionsPerStrip, driftXPos(), simKBmtfDigis_cfi::eLoss, dqmMemoryStats::float, fluctuateCharge, fluctuateEloss(), mps_fire::i, StripTopology::localPitch(), LogDebug, PV3DBase< T, PVType, FrameType >::mag(), StripGeomDetUnit::specificTopology(), theParticleDataTable, and TimeResponse().

75  {
76  // signal after pulse shape correction
77  float const decSignal = TimeResponse(hit, det);
78 
79  // if out of time go home!
80  if (0 == decSignal)
81  return ionization_type();
82 
83  // Get the nass if the particle, in MeV.
84  // Protect from particles with Mass = 0, assuming then the pion mass
85  assert(theParticleDataTable != nullptr);
86  ParticleData const* particle = theParticleDataTable->particle(hit->particleType());
87  double const particleMass = particle ? particle->mass() * 1000 : 139.57;
88  double const particleCharge = particle ? particle->charge() : 1.;
89 
90  if (!particle) {
91  LogDebug("SiLinearChargeDivider") << "Cannot find particle of type " << hit->particleType()
92  << " in the PDT we assign to this particle the mass and charge of the Pion";
93  }
94 
95  int NumberOfSegmentation =
96  // if neutral: just one deposit....
97  (fabs(particleMass) < 1.e-6 || particleCharge == 0)
98  ? 1
99  :
100  // computes the number of segments from number of segments per strip times number of strips.
101  (int)(1 + chargedivisionsPerStrip *
102  fabs(driftXPos(hit->exitPoint(), driftdir, moduleThickness) -
103  driftXPos(hit->entryPoint(), driftdir, moduleThickness)) /
104  det.specificTopology().localPitch(hit->localPosition()));
105 
106  // Eloss in GeV
107  float eLoss = hit->energyLoss();
108 
109  // Prepare output
110  ionization_type _ionization_points;
111  _ionization_points.resize(NumberOfSegmentation);
112 
113  // Fluctuate charge in track subsegments
114  LocalVector direction = hit->exitPoint() - hit->entryPoint();
115  if (NumberOfSegmentation <= 1) {
116  // here I need a random... not 0.5
117  _ionization_points[0] = EnergyDepositUnit(eLoss * decSignal / eLoss, hit->entryPoint() + 0.5f * direction);
118  } else {
119  float eLossVector[NumberOfSegmentation];
120  if (fluctuateCharge) {
121  fluctuateEloss(particleMass, hit->pabs(), eLoss, direction.mag(), NumberOfSegmentation, eLossVector, engine);
122  // Save the energy of each segment
123  for (int i = 0; i != NumberOfSegmentation; i++) {
124  // take energy value from vector eLossVector,
125  _ionization_points[i] =
126  EnergyDepositUnit(eLossVector[i] * decSignal / eLoss,
127  hit->entryPoint() + float((i + 0.5) / NumberOfSegmentation) * direction);
128  }
129  } else {
130  // Save the energy of each segment
131  for (int i = 0; i != NumberOfSegmentation; i++) {
132  // take energy value from eLoss average over n.segments.
133  _ionization_points[i] =
134  EnergyDepositUnit(decSignal / float(NumberOfSegmentation),
135  hit->entryPoint() + float((i + 0.5) / NumberOfSegmentation) * direction);
136  }
137  }
138  }
139  return _ionization_points;
140 }
const ParticleDataTable * theParticleDataTable
float TimeResponse(const PSimHit *hit, const StripGeomDetUnit &det)
float driftXPos(const Local3DPoint &pos, const LocalVector &drift, double thickness)
void fluctuateEloss(double const particleMass, float momentum, float eloss, float length, int NumberOfSegmentation, float elossVector[], CLHEP::HepRandomEngine *)
assert(be >=bs)
virtual float localPitch(const LocalPoint &) const =0
T mag() const
Definition: PV3DBase.h:64
virtual const StripTopology & specificTopology() const
Returns a reference to the strip proxy topology.
HepPDT::ParticleData ParticleData
std::vector< EnergyDepositUnit > ionization_type
#define LogDebug(id)

◆ driftXPos()

float SiLinearChargeDivider::driftXPos ( const Local3DPoint pos,
const LocalVector drift,
double  thickness 
)
inlineprivate

Definition at line 57 of file SiLinearChargeDivider.h.

References shallow::drift(), Calorimetry_cff::thickness, PV3DBase< T, PVType, FrameType >::x(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by divide().

57  {
58  return pos.x() + (thickness / 2. - pos.z()) * drift.x() / drift.z();
59  }
T z() const
Definition: PV3DBase.h:61
LocalVector drift(const StripGeomDetUnit *, const MagneticField &, const SiStripLorentzAngle &)
Definition: ShallowTools.cc:36
T x() const
Definition: PV3DBase.h:59

◆ fluctuateEloss()

void SiLinearChargeDivider::fluctuateEloss ( double const  particleMass,
float  momentum,
float  eloss,
float  length,
int  NumberOfSegmentation,
float  elossVector[],
CLHEP::HepRandomEngine *  engine 
)
private

Definition at line 142 of file SiLinearChargeDivider.cc.

References deltaCut, fluctuate, mps_fire::i, cuy::ii, and particleFlowDisplacedVertex_cfi::ratio.

Referenced by divide().

148  {
149  // Generate charge fluctuations.
150  float sum = 0.;
151  double deltaCutoff;
152  double mom = particleMomentum * 1000.;
153  double seglen = length / NumberOfSegs * 10.;
154  double segeloss = (1000. * eloss) / NumberOfSegs;
155  for (int i = 0; i < NumberOfSegs; i++) {
156  // The G4 routine needs momentum in MeV, mass in MeV, delta-cut in MeV,
157  // track segment length in mm, segment eloss in MeV
158  // Returns fluctuated eloss in MeV
159  // the cutoff is sometimes redefined inside, so fix it.
160  deltaCutoff = deltaCut;
161  sum += (elossVector[i] =
162  fluctuate->SampleFluctuations(mom, particleMass, deltaCutoff, seglen, segeloss, engine) / 1000.);
163  }
164 
165  if (sum > 0.) { // If fluctuations give eloss>0.
166  // Rescale to the same total eloss
167  float ratio = eloss / sum;
168  for (int ii = 0; ii < NumberOfSegs; ii++)
169  elossVector[ii] = ratio * elossVector[ii];
170  } else { // If fluctuations gives 0 eloss
171  float averageEloss = eloss / NumberOfSegs;
172  for (int ii = 0; ii < NumberOfSegs; ii++)
173  elossVector[ii] = averageEloss;
174  }
175  return;
176 }
ii
Definition: cuy.py:589
std::unique_ptr< SiG4UniversalFluctuation > fluctuate

◆ readPulseShape()

void SiLinearChargeDivider::readPulseShape ( const std::string &  pulseShapeFileName)
private

Definition at line 33 of file SiLinearChargeDivider.cc.

References funct::abs(), HLT_2022v14_cff::distance, geometryDiff::epsilon, Exception, mps_splice::line, pulseResolution, pulset0Idx, pulseValues, AlCaHLTBitMon_QueryRunRegistry::string, and relativeConstraints::value.

Referenced by SiLinearChargeDivider().

33  {
34  // Pulse shape file format: empty lines and comments (lines starting with '#') are ignored
35  // one line "resolution: value" is interpreted as the resolution
36  // all other lines are read as consecutive values of the shape
37  std::ifstream shapeFile(pulseShapeFileName.c_str());
38  if (!shapeFile.good()) {
39  throw cms::Exception("FileError") << "Problem opening APV Shape file: " << pulseShapeFileName;
40  }
41  pulseResolution = -1.;
43  const std::string resoPrefix{"resolution: "};
44  while (std::getline(shapeFile, line)) {
45  if ((!line.empty()) && (line.substr(1) != "#")) {
46  std::istringstream lStr{line};
47  if (line.substr(0, resoPrefix.size()) == resoPrefix) {
48  lStr.seekg(resoPrefix.size());
49  lStr >> pulseResolution;
50  } else {
51  double value;
52  while (lStr >> value) {
53  pulseValues.push_back(value);
54  }
55  }
56  }
57  }
58  if (pulseValues.empty() || (pulseResolution == -1.)) {
59  throw cms::Exception("WrongAPVPulseShape") << "Problem reading from APV pulse shape file " << pulseShapeFileName
60  << ": " << (pulseValues.empty() ? "no values" : "no resolution");
61  }
62  const auto maxIt = std::max_element(pulseValues.begin(), pulseValues.end());
63  if (std::abs((*maxIt) - 1.) > std::numeric_limits<double>::epsilon()) {
64  throw cms::Exception("WrongAPVPulseShape")
65  << "The max value of the APV pulse shape stored in the text file used in "
66  "SimGeneral/MixingModule/python/SiStripSimParameters_cfi.py is not equal to 1. Need to be fixed.";
67  }
68  pulset0Idx = std::distance(pulseValues.begin(), maxIt);
69 }
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
Definition: value.py:1
std::vector< double > pulseValues

◆ setParticleDataTable()

void SiLinearChargeDivider::setParticleDataTable ( const ParticleDataTable pdt)
inlineoverridevirtual

Implements SiChargeDivider.

Definition at line 40 of file SiLinearChargeDivider.h.

References theParticleDataTable.

40 { theParticleDataTable = pdt; }
const ParticleDataTable * theParticleDataTable

◆ TimeResponse()

float SiLinearChargeDivider::TimeResponse ( const PSimHit hit,
const StripGeomDetUnit det 
)
private

Definition at line 178 of file SiLinearChargeDivider.cc.

References cosmicShift, createfilelist::int, mag(), pulseResolution, pulset0Idx, pulseValues, GeomDet::surface(), Surface::toGlobal(), and x.

Referenced by divide().

178  {
179  // x is difference between the tof and the tof for a photon (reference)
180  // converted into a bin number
181  const auto dTOF = det.surface().toGlobal(hit->localPosition()).mag() / 30. + cosmicShift - hit->tof();
182  const int x = int(dTOF / pulseResolution) + pulset0Idx;
183  if (x < 0 || x >= int(pulseValues.size()))
184  return 0;
185  return hit->energyLoss() * pulseValues[std::size_t(x)];
186 }
GlobalPoint toGlobal(const Point2DBase< Scalar, LocalTag > lp) const
Definition: Surface.h:79
const Plane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:37
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
std::vector< double > pulseValues

Member Data Documentation

◆ chargedivisionsPerStrip

const int SiLinearChargeDivider::chargedivisionsPerStrip
private

Definition at line 46 of file SiLinearChargeDivider.h.

Referenced by divide().

◆ cosmicShift

const double SiLinearChargeDivider::cosmicShift
private

Definition at line 48 of file SiLinearChargeDivider.h.

Referenced by TimeResponse().

◆ deltaCut

const double SiLinearChargeDivider::deltaCut
private

Definition at line 47 of file SiLinearChargeDivider.h.

Referenced by fluctuateEloss().

◆ fluctuate

std::unique_ptr<SiG4UniversalFluctuation> SiLinearChargeDivider::fluctuate
private

Definition at line 55 of file SiLinearChargeDivider.h.

Referenced by fluctuateEloss().

◆ fluctuateCharge

const bool SiLinearChargeDivider::fluctuateCharge
private

Definition at line 45 of file SiLinearChargeDivider.h.

Referenced by divide().

◆ peakMode

const bool SiLinearChargeDivider::peakMode
private

Definition at line 44 of file SiLinearChargeDivider.h.

Referenced by SiLinearChargeDivider().

◆ pulseResolution

double SiLinearChargeDivider::pulseResolution
private

Definition at line 50 of file SiLinearChargeDivider.h.

Referenced by readPulseShape(), and TimeResponse().

◆ pulset0Idx

unsigned int SiLinearChargeDivider::pulset0Idx
private

Definition at line 51 of file SiLinearChargeDivider.h.

Referenced by readPulseShape(), and TimeResponse().

◆ pulseValues

std::vector<double> SiLinearChargeDivider::pulseValues
private

Definition at line 52 of file SiLinearChargeDivider.h.

Referenced by readPulseShape(), and TimeResponse().

◆ theParticleDataTable

const ParticleDataTable* SiLinearChargeDivider::theParticleDataTable
private

Definition at line 49 of file SiLinearChargeDivider.h.

Referenced by divide(), and setParticleDataTable().