CMS 3D CMS Logo

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

#include <CSCGasCollisions.h>

Public Member Functions

 CSCGasCollisions (const edm::ParameterSet &pset)
 
bool dumpGasCollisions (void) const
 
bool saveGasCollisions (void) const
 
void setParticleDataTable (const ParticleDataTable *pdt)
 
void simulate (const PSimHit &, std::vector< LocalPoint > &clusters, std::vector< int > &electrons, CLHEP::HepRandomEngine *)
 
virtual ~CSCGasCollisions ()
 

Static Public Attributes

static const int MAX_STEPS = 400
 
static const int N_ENERGY = 63
 
static const int N_ENTRIES = N_GAMMA * N_ENERGY
 
static const int N_GAMMA = 21
 

Private Member Functions

void fillCollisionsForThisGamma (float, std::vector< float > &) const
 
float generateEnergyLoss (double avCollisions, double anmin, double anmax, const std::vector< float > &collisions, CLHEP::HepRandomEngine *) const
 
double generateStep (double avCollisions, CLHEP::HepRandomEngine *) const
 
void ionize (double energyTransferred, LocalPoint startHere) const
 
float lnEnergyLoss (float, const std::vector< float > &) const
 
void readCollisionTable ()
 
void writeSummary (int n_try, int n_steps, double sum_steps, float dedx, const PSimHit &simhit) const
 

Private Attributes

double clusterExtent
 
double deCut
 
bool dumpGasCollisions_
 
double eion
 
double ework
 
double gasDensity
 
const std::string me
 
bool saveGasCollisions_
 
std::vector< float > theCollisionTable
 
CSCCrossGaptheCrossGap
 
std::vector< float > theEnergyBins
 
std::vector< float > theGammaBins
 
const ParticleDataTabletheParticleDataTable
 

Detailed Description

Class to handle CSC gas ionization simulation during digitization.

Author
Tim Cox

This class contains the port of the old CMSIM/ORCA Fortran routine mc_clus.F, and the GEANT3 code it relied on. Of course the version here is much improved :)

Definition at line 27 of file CSCGasCollisions.h.

Constructor & Destructor Documentation

◆ CSCGasCollisions()

CSCGasCollisions::CSCGasCollisions ( const edm::ParameterSet pset)

Definition at line 58 of file CSCGasCollisions.cc.

References clusterExtent, deCut, dumpGasCollisions(), dumpGasCollisions_, eion, ework, gasDensity, me, muonDTDigis_cfi::pset, readCollisionTable(), and saveGasCollisions().

59  : me("CSCGasCollisions"),
60  gasDensity(2.1416e-03),
61  deCut(1.e05),
62  eion(14.95),
63  ework(34.0),
64  clusterExtent(0.001),
65  theGammaBins(N_GAMMA, 0.),
68  theCrossGap(nullptr),
69  theParticleDataTable(nullptr),
70  saveGasCollisions_(false),
71  dumpGasCollisions_(false) {
72  dumpGasCollisions_ = pset.getUntrackedParameter<bool>("dumpGasCollisions");
73 
74  edm::LogInfo(me) << "Constructing a " << me << ":";
75  edm::LogInfo(me) << "gas density = " << gasDensity << " g/cm3";
76  edm::LogInfo(me) << "max eloss per collision allowed = " << deCut / 1000.
77  << " keV (for higher elosses, hits should have been simulated.)";
78  edm::LogInfo(me) << "ionization threshold = " << eion << " eV";
79  edm::LogInfo(me) << "effective work function = " << ework << " eV";
80  edm::LogInfo(me) << "cluster extent = " << clusterExtent * 1.e04 << " micrometres";
81  edm::LogInfo(me) << "dump gas collision and simhit information? " << dumpGasCollisions();
82  edm::LogInfo(me) << "save gas collision information? -NOT YET IMPLEMENTED- " << saveGasCollisions();
83 
85 }
static const int N_ENERGY
std::vector< float > theGammaBins
const std::string me
bool dumpGasCollisions(void) const
const ParticleDataTable * theParticleDataTable
static const int N_GAMMA
static const int N_ENTRIES
Log< level::Info, false > LogInfo
std::vector< float > theCollisionTable
std::vector< float > theEnergyBins
CSCCrossGap * theCrossGap
bool saveGasCollisions(void) const

◆ ~CSCGasCollisions()

CSCGasCollisions::~CSCGasCollisions ( )
virtual

Definition at line 87 of file CSCGasCollisions.cc.

References me, and theCrossGap.

87  {
88  edm::LogInfo(me) << "Destructing a " << me;
89  delete theCrossGap;
90 }
const std::string me
Log< level::Info, false > LogInfo
CSCCrossGap * theCrossGap

Member Function Documentation

◆ dumpGasCollisions()

bool CSCGasCollisions::dumpGasCollisions ( void  ) const
inline

Definition at line 39 of file CSCGasCollisions.h.

References dumpGasCollisions_.

Referenced by CSCGasCollisions(), generateEnergyLoss(), generateStep(), and simulate().

39 { return dumpGasCollisions_; }

◆ fillCollisionsForThisGamma()

void CSCGasCollisions::fillCollisionsForThisGamma ( float  logGamma,
std::vector< float > &  collisions 
) const
private

Definition at line 550 of file CSCGasCollisions.cc.

References HltBtagPostValidation_cff::c, mps_fire::i, LogTrace, me, N_ENERGY, N_GAMMA, theCollisionTable, and theGammaBins.

Referenced by simulate().

550  {
551  std::vector<float>::const_iterator bigger =
552  find_if(theGammaBins.begin(), theGammaBins.end(), [&logGamma](auto c) { return c > logGamma; });
553 
554  if (bigger == theGammaBins.end()) {
555  // use highest bin
556  LogTrace(me) << ": using highest gamma bin"
557  << " for logGamma = " << logGamma;
558  for (int i = 0; i < N_ENERGY; ++i)
559  collisions[i] = theCollisionTable[i * N_GAMMA];
560  } else {
561  // use bigger and its lower neighbour
562  std::vector<float>::difference_type ihi = bigger - theGammaBins.begin();
563  if (ihi > 0) {
564  double dlg2 = *bigger--; // and decrement after deref
565  // LogTrace(me) << ": using gamma bins "
566  // << ihi-1 << " and " << ihi;
567  double dlg1 = *bigger; // now the preceding element
568  double dlg = (logGamma - dlg1) / (dlg2 - dlg1);
569  double omdlg = 1. - dlg;
570  for (int i = 0; i < N_ENERGY; ++i)
571  collisions[i] = theCollisionTable[i * N_GAMMA + ihi - 1] * omdlg + theCollisionTable[i * N_GAMMA + ihi] * dlg;
572  } else {
573  // bigger has no lower neighbour
574  LogTrace(me) << ": using lowest gamma bin"
575  << " for logGamma = " << logGamma;
576 
577  for (int i = 0; i < N_ENERGY; ++i)
578  collisions[i] = theCollisionTable[i * N_GAMMA];
579  }
580  }
581 }
static const int N_ENERGY
#define LogTrace(id)
std::vector< float > theGammaBins
const std::string me
static const int N_GAMMA
std::vector< float > theCollisionTable

◆ generateEnergyLoss()

float CSCGasCollisions::generateEnergyLoss ( double  avCollisions,
double  anmin,
double  anmax,
const std::vector< float > &  collisions,
CLHEP::HepRandomEngine *  engine 
) const
private

Definition at line 293 of file CSCGasCollisions.cc.

References CSCCrossGap::addEloss(), CSCCrossGap::beta2(), dumpGasCollisions(), JetChargeProducer_cfi::exp, CSCCrossGap::gamma(), lnEnergyLoss(), dqm-mbProfile::log, LogTrace, me, and theCrossGap.

Referenced by simulate().

297  {
298  // Generate a no. of collisions between collisions[0] and [N_ENERGY-1]
299  float lnColl = log(CLHEP::RandFlat::shoot(engine, anmin, anmax));
300 
301  // Without using CLHEP: approx random between anmin and anmax
302  // double ra = double(rand())/double(RAND_MAX)*avCollisions;
303  // cout << "ra = " << ra << std::endl;
304  // float lnColl = static_cast<float>( log( ra ) );
305 
306  // Find energy loss for that number
307  float lnE = lnEnergyLoss(lnColl, collisions);
308  float eloss = exp(lnE);
309  // Compensate if gamma was actually below 1.1
310  if (theCrossGap->gamma() < 1.1)
311  eloss = eloss * 0.173554 / theCrossGap->beta2();
312  LogTrace(me) << "eloss = " << eloss << " eV";
313  // Next line only used to fill container of eloss's for later diagnostic dumps
314  if (dumpGasCollisions())
315  theCrossGap->addEloss(eloss);
316  return eloss;
317 }
void addEloss(float eloss)
Definition: CSCCrossGap.h:46
float lnEnergyLoss(float, const std::vector< float > &) const
double beta2() const
Definition: CSCCrossGap.h:50
#define LogTrace(id)
const std::string me
bool dumpGasCollisions(void) const
double gamma() const
Definition: CSCCrossGap.h:51
CSCCrossGap * theCrossGap

◆ generateStep()

double CSCGasCollisions::generateStep ( double  avCollisions,
CLHEP::HepRandomEngine *  engine 
) const
private

Definition at line 277 of file CSCGasCollisions.cc.

References CSCCrossGap::addStep(), dumpGasCollisions(), LogTrace, me, and theCrossGap.

Referenced by simulate().

277  {
278  // Generate a m.f.p. (1/avCollisions = cm/collision)
279  double step = (CLHEP::RandExponential::shoot(engine)) / avCollisions;
280 
281  // Without using CLHEP: approx random exponential by...
282  // double da = double(rand())/double(RAND_MAX);
283  // double step = -log(1.-da)/avCollisions;
284 
285  LogTrace(me) << " step = " << step << " cm";
286  // Next line only used to fill a container of 'step's for later diagnostic
287  // dumps
288  if (dumpGasCollisions())
290  return step;
291 }
#define LogTrace(id)
const std::string me
bool dumpGasCollisions(void) const
void addStep(double step)
Definition: CSCCrossGap.h:45
CSCCrossGap * theCrossGap
step
Definition: StallMonitor.cc:83

◆ ionize()

void CSCGasCollisions::ionize ( double  energyTransferred,
LocalPoint  startHere 
) const
private

Definition at line 319 of file CSCGasCollisions.cc.

References CSCCrossGap::addCluster(), CSCCrossGap::addElectrons(), CSCCrossGap::addElectronToBack(), clusterExtent, eion, ework, gasDensity, LogTrace, me, CSCCrossGap::noOfClusters(), CSCCrossGap::noOfElectrons(), funct::pow(), FastTimerService_cff::range, theCrossGap, and CSCCrossGap::unitVector().

Referenced by simulate().

319  {
320  while (energyAvailable > eion) {
321  LogTrace(me) << " NEW CLUSTER " << theCrossGap->noOfClusters() + 1 << " AT " << startHere;
322  LocalPoint newCluster(startHere);
323  theCrossGap->addCluster(newCluster);
324 
325  //@@ I consider NOT subtracting eion before calculating range to be a bug.
326  //@@ But this changes tuning of the algorithm so leave it until after the
327  // big rush to 7_5_0
328  //@@ energyAvailable -= eion;
329 
330  // Sauli CERN 77-09: delta e range with E in MeV (Sauli references Kobetich
331  // & Katz 1968 but that has nothing to do with this expression! He seems to
332  // have made a mistake.) I found the Sauli-quoted relationship in R.
333  // Glocker, Z. Naturforsch. Ba, 129 (1948): delta e range R = aE^n with
334  // a=710, n=1.72 for E in MeV and R in mg/cm^2 applicable over the range E =
335  // 0.001 to 0.3 MeV.
336 
337  // Take HALF that range. //@@ Why? Why not...
338  double range = 0.5 * (0.71 / gasDensity) * pow(energyAvailable * 1.E-6, 1.72);
339  LogTrace(me) << " range = " << range << " cm";
340  if (range < clusterExtent) {
341  // short-range delta e
342  // How many electrons can we make? Now use *average* energy for ionization
343  // (not *minimum*)
344  int nelec = static_cast<int>(energyAvailable / ework);
345  LogTrace(me) << "short-range delta energy in = " << energyAvailable << " eV";
346  // energyAvailable -= nelec*(energyAvailable/ework);
347  energyAvailable -= nelec * ework;
348  // If still above eion (minimum, not average) add one more e
349  if (energyAvailable > eion) {
350  ++nelec;
351  energyAvailable -= eion;
352  }
353  LogTrace(me) << "short-range delta energy out = " << energyAvailable << " eV, nelec = " << nelec;
354  theCrossGap->addElectrons(nelec);
355  break;
356 
357  } else {
358  // long-range delta e
359  LogTrace(me) << "long-range delta \n"
360  << "no. of electrons in cluster now = " << theCrossGap->noOfElectrons();
361  theCrossGap->addElectrons(1); // Position is at startHere still
362 
363  bool new_range = false;
364  while (!new_range && (energyAvailable > ework)) {
365  energyAvailable -= ework;
366  while (energyAvailable > eion) {
367  double range2 = 0.5 * 0.71 / gasDensity * pow(1.E-6 * energyAvailable, 1.72);
368  double drange = range - range2;
369  LogTrace(me) << " energy left = " << energyAvailable << " eV, range2 = " << range2
370  << " cm, drange = " << drange << " cm";
371  if (drange < clusterExtent) {
372  theCrossGap->addElectronToBack(); // increment last element
373  } else {
374  startHere += drange * theCrossGap->unitVector(); // update delta e start position
375  range = range2; // update range
376  new_range = true; // Test range again
377  LogTrace(me) << "reset range to range2 = " << range << " from startHere = " << startHere << " and iterate";
378  }
379  break; // out of inner while energyAvailable>eion
380 
381  } // inner while energyAvailable>eion
382 
383  } // while !new_range && energyAvailable>ework
384 
385  // energyAvailable now less than ework, but still may be over eion...add
386  // an e
387  if (energyAvailable > eion) {
388  energyAvailable -= ework; // yes, it may go negative
389  theCrossGap->addElectronToBack(); // add one more e
390  }
391 
392  } // if range
393 
394  } // outer while energyAvailable>eion
395 }
int noOfClusters() const
Definition: CSCCrossGap.h:33
void addElectronToBack()
Definition: CSCCrossGap.h:43
#define LogTrace(id)
const std::string me
int noOfElectrons() const
Definition: CSCCrossGap.h:35
CSCCrossGap * theCrossGap
void addCluster(LocalPoint here)
Definition: CSCCrossGap.h:41
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
LocalVector unitVector() const
Definition: CSCCrossGap.h:53
void addElectrons(int nelec=1)
Definition: CSCCrossGap.h:42

◆ lnEnergyLoss()

float CSCGasCollisions::lnEnergyLoss ( float  lnCollisions,
const std::vector< float > &  collisions 
) const
private

Definition at line 519 of file CSCGasCollisions.cc.

References HltBtagPostValidation_cff::c, spr::find(), LogTrace, me, and theEnergyBins.

Referenced by generateEnergyLoss().

519  {
520  float lnE = -1.;
521 
522  // Find collision[] bin in which lnCollisions falls
523  std::vector<float>::const_iterator it = find(collisions.begin(), collisions.end(), lnCollisions);
524 
525  if (it != collisions.end()) {
526  // found the value
527  std::vector<float>::difference_type ihi = it - collisions.begin();
528  LogTrace(me) << ": using one energy bin " << ihi << " = " << theEnergyBins[ihi]
529  << " for lnCollisions = " << lnCollisions;
530  lnE = theEnergyBins[ihi];
531  } else {
532  // interpolate the value
533  std::vector<float>::const_iterator loside =
534  find_if(collisions.begin(), collisions.end(), [&lnCollisions](auto c) { return c < lnCollisions; });
535  std::vector<float>::difference_type ilo = loside - collisions.begin();
536  if (ilo > 0) {
537  LogTrace(me) << ": using energy bin " << ilo - 1 << " and " << ilo;
538  lnE = theEnergyBins[ilo - 1] + (lnCollisions - collisions[ilo - 1]) *
539  (theEnergyBins[ilo] - theEnergyBins[ilo - 1]) /
540  (collisions[ilo] - collisions[ilo - 1]);
541  } else {
542  LogTrace(me) << ": using one energy bin 0 = " << theEnergyBins[0] << " for lnCollisions = " << lnCollisions;
543  lnE = theEnergyBins[0]; //@@ WHAT ELSE TO DO?
544  }
545  }
546 
547  return lnE;
548 }
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:19
#define LogTrace(id)
const std::string me
std::vector< float > theEnergyBins

◆ readCollisionTable()

void CSCGasCollisions::readCollisionTable ( )
private

Definition at line 92 of file CSCGasCollisions.cc.

References Exception, DeadROC_duringRun::f1, groupFilesInBlocks::fin, mps_fire::i, LogTrace, me, N_ENERGY, N_ENTRIES, N_GAMMA, theCollisionTable, theEnergyBins, and theGammaBins.

Referenced by CSCGasCollisions().

92  {
93  // I'd prefer to allow comments in the data file which means
94  // complications of reading line-by-line and then item-by-item.
95 
96  // Is float OK? Or do I need double?
97 
98  // Do I need any sort of error trapping?
99 
100  // We use the default CMSSW data file path
101  // We use the default file name 'collisions.dat'
102 
103  // This can be reset in .orcarc by SimpleConfigurable
104  // Muon:Endcap:CollisionsFile
105 
106  // TODO make configurable
107  string colliFile = "SimMuon/CSCDigitizer/data/collisions.dat";
108  edm::FileInPath f1{colliFile};
109  edm::LogInfo(me) << ": reading " << f1.fullPath();
110 
111  ifstream fin{f1.fullPath()};
112 
113  if (fin.fail()) {
114  string errorMessage = "Cannot open input file " + f1.fullPath();
115  edm::LogError("CSCGasCollisions") << errorMessage;
116  throw cms::Exception(errorMessage);
117  }
118 
119  fin.clear(); // Clear eof read status
120  fin.seekg(0, ios::beg); // Position at start of file
121 
122  // @@ We had better have the right sizes everywhere or all
123  // hell will break loose. There's no trapping.
124  LogTrace(me) << "Reading gamma bins";
125  for (int i = 0; i < N_GAMMA; ++i) {
126  fin >> theGammaBins[i];
127  LogTrace(me) << i + 1 << " " << theGammaBins[i];
128  }
129 
130  LogTrace(me) << "Reading energy bins \n";
131  for (int i = 0; i < N_ENERGY; ++i) {
132  fin >> theEnergyBins[i];
133  LogTrace(me) << i + 1 << " " << theEnergyBins[i];
134  }
135 
136  LogTrace(me) << "Reading collisions table \n";
137 
138  for (int i = 0; i < N_ENTRIES; ++i) {
139  fin >> theCollisionTable[i];
140  LogTrace(me) << i + 1 << " " << theCollisionTable[i];
141  }
142 
143  fin.close();
144 }
static const int N_ENERGY
Log< level::Error, false > LogError
#define LogTrace(id)
std::vector< float > theGammaBins
const std::string me
static const int N_GAMMA
static const int N_ENTRIES
Log< level::Info, false > LogInfo
std::vector< float > theCollisionTable
std::vector< float > theEnergyBins

◆ saveGasCollisions()

bool CSCGasCollisions::saveGasCollisions ( void  ) const
inline

Definition at line 40 of file CSCGasCollisions.h.

References saveGasCollisions_.

Referenced by CSCGasCollisions().

40 { return saveGasCollisions_; }

◆ setParticleDataTable()

void CSCGasCollisions::setParticleDataTable ( const ParticleDataTable pdt)

Definition at line 146 of file CSCGasCollisions.cc.

References theParticleDataTable.

Referenced by CSCWireHitSim::setParticleDataTable().

146 { theParticleDataTable = pdt; }
const ParticleDataTable * theParticleDataTable

◆ simulate()

void CSCGasCollisions::simulate ( const PSimHit simhit,
std::vector< LocalPoint > &  clusters,
std::vector< int > &  electrons,
CLHEP::HepRandomEngine *  engine 
)

Definition at line 148 of file CSCGasCollisions.cc.

References cms::cuda::assert(), deCut, dumpGasCollisions(), eion, pwdgSkimBPark_cfi::electrons, CSCCrossGap::electrons(), PSimHit::energyLoss(), PSimHit::entryPoint(), PSimHit::exitPoint(), JetChargeProducer_cfi::exp, fillCollisionsForThisGamma(), generateEnergyLoss(), generateStep(), CSCCrossGap::ionClusters(), ionize(), CSCCrossGap::length(), CSCCrossGap::logGamma(), LogTrace, EgHLTOffHistBins_cfi::mass, MAX_STEPS, me, N_ENERGY, PSimHit::pabs(), PSimHit::particleType(), theCrossGap, theParticleDataTable, CSCCrossGap::unitVector(), and writeSummary().

Referenced by CSCWireHitSim::getIonizationClusters().

151  {
152  const float epsilonL = 0.01; // Shortness of simhit 'length'
153  // const float max_gap_z = 1.5; // Gas gaps are 0.5
154  // or 1.0 cm
155 
156  // Note that what I call the 'gap' may in fact be the 'length' of a PSimHit
157  // which does not start and end on the gap edges. This confuses the
158  // nomenclature at least.
159 
160  double mom = simhit.pabs(); // in GeV/c - see MuonSensitiveDetector.cc
161  // int iam = simhit.particleType(); // PDG type
162  delete theCrossGap; // before building new one
163  assert(theParticleDataTable != nullptr);
164  ParticleData const *particle = theParticleDataTable->particle(simhit.particleType());
165  double mass = 0.105658; // assume a muon
166  if (particle == nullptr) {
167  edm::LogError("CSCGasCollisions") << "Cannot find particle of type " << simhit.particleType() << " in the PDT";
168  } else {
169  mass = particle->mass();
170  LogTrace(me) << "[CSCGasCollisions] Found particle of type " << simhit.particleType()
171  << " in the PDT; mass = " << mass;
172  }
173 
174  theCrossGap = new CSCCrossGap(mass, mom, simhit.exitPoint() - simhit.entryPoint());
175  float gapSize = theCrossGap->length();
176 
177  // Test the simhit 'length' (beware of angular effects)
178  // if ( gapSize <= epsilonL || gapSize > max_gap_z ) {
179  if (gapSize <= epsilonL) {
180  edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] WARNING! simhit entry and exit are too close - "
181  "skipping simhit:"
182  << "\n entry = " << simhit.entryPoint() << ": exit = " << simhit.exitPoint()
183  << "\n particle type = " << simhit.particleType()
184  << " : momentum = " << simhit.pabs()
185  << " GeV/c : energy loss = " << simhit.energyLoss() * 1.E06 << " keV"
186  << ", gapSize = " << gapSize << " cm (< epsilonL = " << epsilonL << " cm)";
187  return; //@@ Just skip this PSimHit
188  }
189 
190  // Interpolate the table for current gamma value
191  // Extract collisions binned by energy loss values, for this gamma
192  std::vector<float> collisions(N_ENERGY);
193  double loggam = theCrossGap->logGamma();
194  fillCollisionsForThisGamma(static_cast<float>(loggam), collisions);
195 
196  double anmin = exp(collisions[N_ENERGY - 1]);
197  double anmax = exp(collisions[0]);
198  double amu = anmax - anmin;
199 
200  LogTrace(me) << "collisions extremes = " << collisions[N_ENERGY - 1] << ", " << collisions[0] << "\n"
201  << "anmin = " << anmin << ", anmax = " << anmax << "\n"
202  << "amu = " << amu << "\n";
203 
204  float dedx = 0.; // total energy loss
205  double sum_steps = 0.; // total distance across gap (along simhit direction)
206  int n_steps = 0; // no. of steps/primary collisions
207  int n_try = 0; // no. of tries to generate steps
208  double step = -1.; // Sentinel for start
209 
210  LocalPoint layerLocalPoint(simhit.entryPoint());
211 
212  // step/primary collision loop
213  while (sum_steps < gapSize) {
214  ++n_try;
215  if (n_try > MAX_STEPS) {
216  int maxst = MAX_STEPS;
217  edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] WARNING! n_try = " << n_try
218  << " is > MAX_STEPS = " << maxst << " - skipping simhit:"
219  << "\n entry = " << simhit.entryPoint() << ": exit = " << simhit.exitPoint()
220  << "\n particle type = " << simhit.particleType()
221  << " : momentum = " << simhit.pabs()
222  << " GeV/c : energy loss = " << simhit.energyLoss() * 1.E06 << " keV"
223  << "\n gapSize = " << gapSize << " cm, last step = " << step
224  << " cm, sum_steps = " << sum_steps << " cm, n_steps = " << n_steps;
225  break;
226  }
227  step = generateStep(amu, engine);
228  if (sum_steps + step > gapSize)
229  break; // this step goes too far
230 
231  float eloss = generateEnergyLoss(amu, anmin, anmax, collisions, engine);
232 
233  // Is the eloss too large? (then GEANT should have produced hits!)
234  if (eloss > deCut) {
235  edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] WARNING! eloss = " << eloss << " eV is too large (> "
236  << deCut << " eV) - trying another collision";
237  continue; // to generate another collision/step
238  }
239 
240  dedx += eloss; // the energy lost from the ionizing particle
241  sum_steps += step; // the position of the ionizing particle
242  ++n_steps; // the number of primary collisions
243 
244  if (n_steps > MAX_STEPS) { // Extra-careful trap for bizarreness
245  edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] WARNING! " << n_steps
246  << " is too many steps - skipping simhit:"
247  << "\n entry = " << simhit.entryPoint() << ": exit = " << simhit.exitPoint()
248  << "\n particle type = " << simhit.particleType()
249  << " : momentum = " << simhit.pabs()
250  << " GeV/c : energy loss = " << simhit.energyLoss() * 1.E06 << " keV"
251  << "\ngapSize=" << gapSize << " cm, last step=" << step
252  << " cm, sum_steps=" << sum_steps << " cm";
253  break;
254  }
255  LogTrace(me) << "sum_steps = " << sum_steps << " cm , dedx = " << dedx << " eV";
256 
257  // Generate ionization.
258  // eion is the minimum energy at which ionization can occur in the gas
259  if (eloss > eion) {
260  layerLocalPoint += step * theCrossGap->unitVector(); // local point where the collision occurs
261  ionize(eloss, layerLocalPoint);
262  } else {
263  LogTrace(me) << "Energy available = " << eloss << " eV is too low for ionization (< eion = " << eion << " eV)";
264  }
265 
266  } // step/collision loop
267 
268  if (dumpGasCollisions())
269  writeSummary(n_try, n_steps, sum_steps, dedx, simhit);
270 
271  // Return values in two container arguments
272  positions = theCrossGap->ionClusters();
274  return;
275 }
Log< level::Info, true > LogVerbatim
float generateEnergyLoss(double avCollisions, double anmin, double anmax, const std::vector< float > &collisions, CLHEP::HepRandomEngine *) const
void writeSummary(int n_try, int n_steps, double sum_steps, float dedx, const PSimHit &simhit) const
static const int N_ENERGY
Local3DPoint exitPoint() const
Exit point in the local Det frame.
Definition: PSimHit.h:48
Log< level::Error, false > LogError
assert(be >=bs)
#define LogTrace(id)
const std::string me
bool dumpGasCollisions(void) const
const ParticleDataTable * theParticleDataTable
double generateStep(double avCollisions, CLHEP::HepRandomEngine *) const
Local3DPoint entryPoint() const
Entry point in the local Det frame.
Definition: PSimHit.h:45
static const int MAX_STEPS
HepPDT::ParticleData ParticleData
double logGamma(double mass, float momentum)
Definition: CSCCrossGap.cc:14
void fillCollisionsForThisGamma(float, std::vector< float > &) const
std::vector< int > electrons() const
Definition: CSCCrossGap.h:34
std::vector< LocalPoint > ionClusters() const
Definition: CSCCrossGap.h:32
float pabs() const
fast and more accurate access to momentumAtEntry().mag()
Definition: PSimHit.h:69
float energyLoss() const
The energy deposit in the PSimHit, in ???.
Definition: PSimHit.h:81
void ionize(double energyTransferred, LocalPoint startHere) const
CSCCrossGap * theCrossGap
step
Definition: StallMonitor.cc:83
int particleType() const
Definition: PSimHit.h:91
LocalVector unitVector() const
Definition: CSCCrossGap.h:53
float length() const
Definition: CSCCrossGap.h:54

◆ writeSummary()

void CSCGasCollisions::writeSummary ( int  n_try,
int  n_steps,
double  sum_steps,
float  dedx,
const PSimHit simhit 
) const
private

Definition at line 397 of file CSCGasCollisions.cc.

References HltBtagPostValidation_cff::c, eion, pwdgSkimBPark_cfi::electrons, CSCCrossGap::electrons(), CSCCrossGap::eLossPerStep(), PSimHit::energyLoss(), PSimHit::entryPoint(), PSimHit::exitPoint(), spr::find(), nano_mu_digi_cff::float, mps_fire::i, CSCCrossGap::ionClusters(), MAX_STEPS, PSimHit::pabs(), PSimHit::particleType(), CSCCrossGap::stepLengths(), customisers::steps, and theCrossGap.

Referenced by simulate().

397  {
398  // Switched from std::cout to LogVerbatim, Mar 2015
399 
400  std::vector<LocalPoint> ion_clusters = theCrossGap->ionClusters();
401  std::vector<int> electrons = theCrossGap->electrons();
402  std::vector<float> elosses = theCrossGap->eLossPerStep();
403  std::vector<double> steps = theCrossGap->stepLengths();
404 
405  edm::LogVerbatim("CSCDigitizer") << "------------------"
406  << "\nAFTER CROSSING GAP";
407  /*
408  edm::LogVerbatim("CSCDigitizer") << "no. of steps tried = " <<
409  n_try
410  << "\nno. of steps from theCrossGap = " <<
411  theCrossGap->noOfSteps()
412  << "\nsize of steps vector = " <<
413  steps.size();
414 
415  edm::LogVerbatim("CSCDigitizer") << "no. of collisions (steps) = " <<
416  n_steps
417  << "\nsize of elosses vector = " <<
418  elosses.size()
419  << "\nsize of ion clusters vector = " <<
420  ion_clusters.size()
421  << "\nsize of electrons vector = " <<
422  electrons.size();
423  */
424 
425  size_t nsteps = n_steps; // force ridiculous type conversion
426  size_t mstep = steps.size() - 1; // final step gets filled but is outside gas
427  // gap - unless we reach MAX_STEPS
428  if ((nsteps != MAX_STEPS) && (nsteps != mstep)) {
429  edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] WARNING! no. of steps = " << nsteps
430  << " .ne. steps.size()-1 = " << mstep;
431  }
432  size_t meloss = elosses.size();
433 
434  if (nsteps != meloss) {
435  edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] WARNING! no. of steps = " << nsteps
436  << " .ne. no. of elosses = " << meloss;
437  } else {
438  edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] # / length of step / energy loss per collision:";
439  for (size_t i = 0; i != nsteps; ++i) {
440  edm::LogVerbatim("CSCDigitizer") << i + 1 << " / S: " << steps[i] << " / E: " << elosses[i];
441  }
442  }
443 
444  size_t mclus = ion_clusters.size();
445  size_t melec = electrons.size();
446 
447  if (mclus != melec) {
448  edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] WARNING! size of cluster vector = " << mclus
449  << " .ne. size of electrons vector = " << melec;
450  } else {
451  edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] # / postion of "
452  "cluster / electrons per cluster: ";
453  for (size_t i = 0; i != mclus; ++i) {
454  edm::LogVerbatim("CSCDigitizer") << i + 1 << " / I: " << ion_clusters[i] << " / E: " << electrons[i];
455  }
456  }
457 
458  int n_ic = count_if(elosses.begin(), elosses.end(), [&](auto c) { return c > this->eion; });
459 
460  edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] total no. of collision steps = " << n_steps;
461  edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] total sum of steps = " << sum_steps << " cm";
462  if (nsteps > 0)
463  edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] average step length = " << sum_steps / float(nsteps)
464  << " cm";
465  edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] total energy loss across gap = " << dedx
466  << " eV = " << dedx / 1000. << " keV";
467  edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] no. of primary ionizing collisions across gap = "
468  "no. with eloss > eion = "
469  << eion << " eV = " << n_ic;
470  if (nsteps > 0)
471  edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] average energy loss/collision = " << dedx / float(nsteps)
472  << " eV";
473 
474  std::vector<int>::const_iterator bigger = find(electrons.begin(), electrons.end(), 0);
475  if (bigger != electrons.end()) {
476  edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] TROUBLE! There is a cluster with 0 electrons.";
477  }
478 
479  int n_e = accumulate(electrons.begin(), electrons.end(), 0);
480 
481  edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] SUMMARY: simhit"
482  << "\n entry = " << simhit.entryPoint() << ": exit = " << simhit.exitPoint()
483  << "\n particle type = " << simhit.particleType()
484  << " : momentum = " << simhit.pabs()
485  << " GeV/c : energy loss = " << simhit.energyLoss() * 1.E06 << " keV";
486 
487  if (nsteps > 0) {
488  edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] SUMMARY: ionization"
489  << " : steps= " << nsteps << " : sum(steps)= " << sum_steps
490  << " cm : <step>= " << sum_steps / float(nsteps) << " cm"
491  << " : ionizing= " << n_ic << " : ionclus= " << mclus << " : total e= " << n_e
492  << " : <dedx/step>= " << dedx / float(nsteps)
493  << " eV : <e/ionclus>= " << float(n_e) / float(mclus)
494  << " : dedx= " << dedx / 1000. << " keV";
495  } else {
496  edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisons] ERROR? no collision steps!";
497  }
498 
499  // Turn off output file - used for initial development
500  // if ( saveGasCollisions() ) {
501  // ofstream of0("osteplen.dat",ios::app);
502  // std::copy( steps.begin(), steps.end(),
503  // std::ostream_iterator<float>(of0,"\n"));
504 
505  // ofstream of1("olperc.dat",ios::app);
506  // std::copy( elosses.begin(), elosses.end(),
507  // std::ostream_iterator<float>(of1,"\n"));
508 
509  // ofstream of2("oclpos.dat",ios::app);
510  // std::copy( ion_clusters.begin(), ion_clusters.end(),
511  // std::ostream_iterator<LocalPoint>(of2,"\n"));
512 
513  // ofstream of3("oepercl.dat",ios::app);
514  // std::copy( electrons.begin(), electrons.end(),
515  // std::ostream_iterator<int>(of3,"\n"));
516  // }
517 }
Log< level::Info, true > LogVerbatim
std::vector< double > stepLengths() const
Definition: CSCCrossGap.h:36
Local3DPoint exitPoint() const
Exit point in the local Det frame.
Definition: PSimHit.h:48
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:19
std::vector< float > eLossPerStep() const
Definition: CSCCrossGap.h:38
Local3DPoint entryPoint() const
Entry point in the local Det frame.
Definition: PSimHit.h:45
static const int MAX_STEPS
std::vector< int > electrons() const
Definition: CSCCrossGap.h:34
std::vector< LocalPoint > ionClusters() const
Definition: CSCCrossGap.h:32
float pabs() const
fast and more accurate access to momentumAtEntry().mag()
Definition: PSimHit.h:69
float energyLoss() const
The energy deposit in the PSimHit, in ???.
Definition: PSimHit.h:81
CSCCrossGap * theCrossGap
int particleType() const
Definition: PSimHit.h:91

Member Data Documentation

◆ clusterExtent

double CSCGasCollisions::clusterExtent
private

Definition at line 72 of file CSCGasCollisions.h.

Referenced by CSCGasCollisions(), and ionize().

◆ deCut

double CSCGasCollisions::deCut
private

Definition at line 69 of file CSCGasCollisions.h.

Referenced by CSCGasCollisions(), and simulate().

◆ dumpGasCollisions_

bool CSCGasCollisions::dumpGasCollisions_
private

Definition at line 83 of file CSCGasCollisions.h.

Referenced by CSCGasCollisions(), and dumpGasCollisions().

◆ eion

double CSCGasCollisions::eion
private

Definition at line 70 of file CSCGasCollisions.h.

Referenced by CSCGasCollisions(), ionize(), simulate(), and writeSummary().

◆ ework

double CSCGasCollisions::ework
private

Definition at line 71 of file CSCGasCollisions.h.

Referenced by CSCGasCollisions(), and ionize().

◆ gasDensity

double CSCGasCollisions::gasDensity
private

Definition at line 63 of file CSCGasCollisions.h.

Referenced by CSCGasCollisions(), and ionize().

◆ MAX_STEPS

const int CSCGasCollisions::MAX_STEPS = 400
static

Definition at line 45 of file CSCGasCollisions.h.

Referenced by simulate(), and writeSummary().

◆ me

const std::string CSCGasCollisions::me
private

◆ N_ENERGY

const int CSCGasCollisions::N_ENERGY = 63
static

Definition at line 43 of file CSCGasCollisions.h.

Referenced by fillCollisionsForThisGamma(), readCollisionTable(), and simulate().

◆ N_ENTRIES

const int CSCGasCollisions::N_ENTRIES = N_GAMMA * N_ENERGY
static

Definition at line 44 of file CSCGasCollisions.h.

Referenced by readCollisionTable().

◆ N_GAMMA

const int CSCGasCollisions::N_GAMMA = 21
static

Definition at line 42 of file CSCGasCollisions.h.

Referenced by fillCollisionsForThisGamma(), and readCollisionTable().

◆ saveGasCollisions_

bool CSCGasCollisions::saveGasCollisions_
private

Definition at line 81 of file CSCGasCollisions.h.

Referenced by saveGasCollisions().

◆ theCollisionTable

std::vector<float> CSCGasCollisions::theCollisionTable
private

Definition at line 77 of file CSCGasCollisions.h.

Referenced by fillCollisionsForThisGamma(), and readCollisionTable().

◆ theCrossGap

CSCCrossGap* CSCGasCollisions::theCrossGap
private

◆ theEnergyBins

std::vector<float> CSCGasCollisions::theEnergyBins
private

Definition at line 76 of file CSCGasCollisions.h.

Referenced by lnEnergyLoss(), and readCollisionTable().

◆ theGammaBins

std::vector<float> CSCGasCollisions::theGammaBins
private

Definition at line 75 of file CSCGasCollisions.h.

Referenced by fillCollisionsForThisGamma(), and readCollisionTable().

◆ theParticleDataTable

const ParticleDataTable* CSCGasCollisions::theParticleDataTable
private

Definition at line 80 of file CSCGasCollisions.h.

Referenced by setParticleDataTable(), and simulate().