CSCDigitizer/src/CSCGasCollisions.cc

0001 // This implements the CSC gas ionization simulation.
0002 // It requires the GEANT3-generated tables in the input file
0003 // (default) MuonData/EndcapData/collisions.dat
0004
0005 // Outstanding problems to fix 15-Oct-2003
0006 // - ework must be re-tuned once eion or ework are removed before delta e range
0007 // estimation.
0008 // - logic of long-range delta e's must be revisited.
0009 // - The code here needs to have diagnostic output removed, or reduced.
0010 // PARTICULARLY filling of std::vectors!
0011 // - 'gap' is a misnomer, when simhit entry and exit don't coincide with gap
0012 // edges
0013 //  so the CSCCrossGap might better be named as a simhit-related thing.
0014 // 22-Jan-2004 Corrected position of trap for 'infinite' loop while
0015 //    generating steps. Output files (debugV-flagged only) require SC flag.
0016 //    Increase deCut from 10 keV to 100 keV to accomodate protons!
0017 // Mar-2015: cout to LogVerbatim. Change superseded debugV-flagged code to flag
0018 // from config.
0019
0020 #include "DataFormats/GeometryVector/interface/LocalVector.h"
0021 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0022 #include "FWCore/ParameterSet/interface/FileInPath.h"
0023 #include "FWCore/Utilities/interface/Exception.h"
0024 #include "SimMuon/CSCDigitizer/src/CSCGasCollisions.h"
0025
0026 #include "CLHEP/Random/RandExponential.h"
0027 #include "CLHEP/Random/RandFlat.h"
0028
0029 // Only needed if file writing code is activated again
0030 // #include <iostream>
0031 #include <cassert>
0032 #include <fstream>
0033
0034 // stdlib math functions
0035 // for 'abs':
0036 #include <cstdlib>
0037 // for usual math functions:
0038 #include <cmath>
0039
0040 // stdlib container trickery
0041 // for 'for_each':
0042 #include <algorithm>
0043 // for 'greater' and 'less':
0044 #include <functional>
0045 // for 'accumulate':
0046 #include <iterator>
0047 #include <numeric>
0048
0049 using namespace std;
0050
0051 /* Gas mixture is Ar/CO2/CF4 = 40/50/10
0052 We'll use the ionization energies
0053 Ar    15.8 eV
0054 CO2   13.7 eV
0055 CF4    17.8
0056 to arrive at a weighted average of eion = 14.95 */
0057
0058 CSCGasCollisions::CSCGasCollisions(const edm::ParameterSet &pset)
0059     : me("CSCGasCollisions"),
0060       gasDensity(2.1416e-03),
0061       deCut(1.e05),
0062       eion(14.95),
0063       ework(34.0),
0064       clusterExtent(0.001),
0065       theGammaBins(N_GAMMA, 0.),
0066       theEnergyBins(N_ENERGY, 0.),
0067       theCollisionTable(N_ENTRIES, 0.),
0068       theCrossGap(nullptr),
0069       theParticleDataTable(nullptr),
0070       saveGasCollisions_(false),
0071       dumpGasCollisions_(false) {
0072   dumpGasCollisions_ = pset.getUntrackedParameter<bool>("dumpGasCollisions");
0073
0074   edm::LogInfo(me) << "Constructing a " << me << ":";
0075   edm::LogInfo(me) << "gas density = " << gasDensity << " g/cm3";
0076   edm::LogInfo(me) << "max eloss per collision allowed = " << deCut / 1000.
0077                    << " keV (for higher elosses, hits should have been simulated.)";
0078   edm::LogInfo(me) << "ionization threshold = " << eion << " eV";
0079   edm::LogInfo(me) << "effective work function = " << ework << " eV";
0080   edm::LogInfo(me) << "cluster extent = " << clusterExtent * 1.e04 << " micrometres";
0081   edm::LogInfo(me) << "dump gas collision and simhit information? " << dumpGasCollisions();
0082   edm::LogInfo(me) << "save gas collision information? -NOT YET IMPLEMENTED- " << saveGasCollisions();
0083
0084   readCollisionTable();
0085 }
0086
0087 CSCGasCollisions::~CSCGasCollisions() {
0088   edm::LogInfo(me) << "Destructing a " << me;
0089   delete theCrossGap;
0090 }
0091
0092 void CSCGasCollisions::readCollisionTable() {
0093   // I'd prefer to allow comments in the data file which means
0094   // complications of reading line-by-line and then item-by-item.
0095
0096   // Is float OK? Or do I need double?
0097
0098   // Do I need any sort of error trapping?
0099
0100   // We use the default CMSSW data file path
0101   // We use the default file name 'collisions.dat'
0102
0103   // This can be reset in .orcarc by SimpleConfigurable
0104   //    Muon:Endcap:CollisionsFile
0105
0106   // TODO make configurable
0107   string colliFile = "SimMuon/CSCDigitizer/data/collisions.dat";
0108   edm::FileInPath f1{colliFile};
0109   edm::LogInfo(me) << ": reading " << f1.fullPath();
0110
0111   ifstream fin{f1.fullPath()};
0112
0113   if (fin.fail()) {
0114     string errorMessage = "Cannot open input file " + f1.fullPath();
0115     edm::LogError("CSCGasCollisions") << errorMessage;
0116     throw cms::Exception(errorMessage);
0117   }
0118
0119   fin.clear();             // Clear eof read status
0120   fin.seekg(0, ios::beg);  // Position at start of file
0121
0122   // @@ We had better have the right sizes everywhere or all
0123   // hell will break loose. There's no trapping.
0124   LogTrace(me) << "Reading gamma bins";
0125   for (int i = 0; i < N_GAMMA; ++i) {
0126     fin >> theGammaBins[i];
0127     LogTrace(me) << i + 1 << "   " << theGammaBins[i];
0128   }
0129
0130   LogTrace(me) << "Reading energy bins \n";
0131   for (int i = 0; i < N_ENERGY; ++i) {
0132     fin >> theEnergyBins[i];
0133     LogTrace(me) << i + 1 << "   " << theEnergyBins[i];
0134   }
0135
0136   LogTrace(me) << "Reading collisions table \n";
0137
0138   for (int i = 0; i < N_ENTRIES; ++i) {
0139     fin >> theCollisionTable[i];
0140     LogTrace(me) << i + 1 << "   " << theCollisionTable[i];
0141   }
0142
0143   fin.close();
0144 }
0145
0146 void CSCGasCollisions::setParticleDataTable(const ParticleDataTable *pdt) { theParticleDataTable = pdt; }
0147
0148 void CSCGasCollisions::simulate(const PSimHit &simhit,
0149                                 std::vector<LocalPoint> &positions,
0150                                 std::vector<int> &electrons,
0151                                 CLHEP::HepRandomEngine *engine) {
0152   const float epsilonL = 0.01;  // Shortness of simhit 'length'
0153   //  const float max_gap_z = 1.5;                     // Gas gaps are 0.5
0154   //  or 1.0 cm
0155
0156   // Note that what I call the 'gap' may in fact be the 'length' of a PSimHit
0157   // which does not start and end on the gap edges. This confuses the
0158   // nomenclature at least.
0159
0160   double mom = simhit.pabs();  // in GeV/c - see MuonSensitiveDetector.cc
0161   //  int iam       = simhit.particleType();           // PDG type
0162   delete theCrossGap;  // before building new one
0163   assert(theParticleDataTable != nullptr);
0164   ParticleData const *particle = theParticleDataTable->particle(simhit.particleType());
0165   double mass = 0.105658;  // assume a muon
0166   if (particle == nullptr) {
0167     edm::LogError("CSCGasCollisions") << "Cannot find particle of type " << simhit.particleType() << " in the PDT";
0168   } else {
0169     mass = particle->mass();
0170     LogTrace(me) << "[CSCGasCollisions] Found particle of type " << simhit.particleType()
0171                  << " in the PDT; mass = " << mass;
0172   }
0173
0174   theCrossGap = new CSCCrossGap(mass, mom, simhit.exitPoint() - simhit.entryPoint());
0175   float gapSize = theCrossGap->length();
0176
0177   // Test the simhit 'length' (beware of angular effects)
0178   //  if ( gapSize <= epsilonL || gapSize > max_gap_z ) {
0179   if (gapSize <= epsilonL) {
0180     edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] WARNING! simhit entry and exit are too close - "
0181                                         "skipping simhit:"
0182                                      << "\n entry = " << simhit.entryPoint() << ": exit  = " << simhit.exitPoint()
0183                                      << "\n particle type = " << simhit.particleType()
0184                                      << " : momentum = " << simhit.pabs()
0185                                      << " GeV/c : energy loss = " << simhit.energyLoss() * 1.E06 << " keV"
0186                                      << ", gapSize = " << gapSize << " cm (< epsilonL = " << epsilonL << " cm)";
0187     return;  //@@ Just skip this PSimHit
0188   }
0189
0190   // Interpolate the table for current gamma value
0191   // Extract collisions binned by energy loss values, for this gamma
0192   std::vector<float> collisions(N_ENERGY);
0193   double loggam = theCrossGap->logGamma();
0194   fillCollisionsForThisGamma(static_cast<float>(loggam), collisions);
0195
0196   double anmin = exp(collisions[N_ENERGY - 1]);
0197   double anmax = exp(collisions[0]);
0198   double amu = anmax - anmin;
0199
0200   LogTrace(me) << "collisions extremes = " << collisions[N_ENERGY - 1] << ", " << collisions[0] << "\n"
0201                << "anmin = " << anmin << ", anmax = " << anmax << "\n"
0202                << "amu = " << amu << "\n";
0203
0204   float dedx = 0.;        // total energy loss
0205   double sum_steps = 0.;  // total distance across gap (along simhit direction)
0206   int n_steps = 0;        // no. of steps/primary collisions
0207   int n_try = 0;          // no. of tries to generate steps
0208   double step = -1.;      // Sentinel for start
0209
0210   LocalPoint layerLocalPoint(simhit.entryPoint());
0211
0212   // step/primary collision loop
0213   while (sum_steps < gapSize) {
0214     ++n_try;
0215     if (n_try > MAX_STEPS) {
0216       int maxst = MAX_STEPS;
0217       edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] WARNING! n_try = " << n_try
0218                                        << " is > MAX_STEPS = " << maxst << " - skipping simhit:"
0219                                        << "\n entry = " << simhit.entryPoint() << ": exit  = " << simhit.exitPoint()
0220                                        << "\n particle type = " << simhit.particleType()
0221                                        << " : momentum = " << simhit.pabs()
0222                                        << " GeV/c : energy loss = " << simhit.energyLoss() * 1.E06 << " keV"
0223                                        << "\n gapSize = " << gapSize << " cm, last step = " << step
0224                                        << " cm, sum_steps = " << sum_steps << " cm, n_steps = " << n_steps;
0225       break;
0226     }
0227     step = generateStep(amu, engine);
0228     if (sum_steps + step > gapSize)
0229       break;  // this step goes too far
0230
0231     float eloss = generateEnergyLoss(amu, anmin, anmax, collisions, engine);
0232
0233     // Is the eloss too large? (then GEANT should have produced hits!)
0234     if (eloss > deCut) {
0235       edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] WARNING! eloss = " << eloss << " eV is too large (> "
0236                                        << deCut << " eV) - trying another collision";
0237       continue;  // to generate another collision/step
0238     }
0239
0240     dedx += eloss;      // the energy lost from the ionizing particle
0241     sum_steps += step;  // the position of the ionizing particle
0242     ++n_steps;          // the number of primary collisions
0243
0244     if (n_steps > MAX_STEPS) {  // Extra-careful trap for bizarreness
0245       edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] WARNING! " << n_steps
0246                                        << " is too many steps -  skipping simhit:"
0247                                        << "\n entry = " << simhit.entryPoint() << ": exit  = " << simhit.exitPoint()
0248                                        << "\n particle type = " << simhit.particleType()
0249                                        << " : momentum = " << simhit.pabs()
0250                                        << " GeV/c : energy loss = " << simhit.energyLoss() * 1.E06 << " keV"
0251                                        << "\ngapSize=" << gapSize << " cm, last step=" << step
0252                                        << " cm, sum_steps=" << sum_steps << " cm";
0253       break;
0254     }
0255     LogTrace(me) << "sum_steps = " << sum_steps << " cm , dedx = " << dedx << " eV";
0256
0257     // Generate ionization.
0258     // eion is the minimum energy at which ionization can occur in the gas
0259     if (eloss > eion) {
0260       layerLocalPoint += step * theCrossGap->unitVector();  // local point where the collision occurs
0261       ionize(eloss, layerLocalPoint);
0262     } else {
0263       LogTrace(me) << "Energy available = " << eloss << " eV is too low for ionization (< eion = " << eion << " eV)";
0264     }
0265
0266   }  // step/collision loop
0267
0268   if (dumpGasCollisions())
0269     writeSummary(n_try, n_steps, sum_steps, dedx, simhit);
0270
0271   // Return values in two container arguments
0272   positions = theCrossGap->ionClusters();
0273   electrons = theCrossGap->electrons();
0274   return;
0275 }
0276
0277 double CSCGasCollisions::generateStep(double avCollisions, CLHEP::HepRandomEngine *engine) const {
0278   // Generate a m.f.p.  (1/avCollisions = cm/collision)
0279   double step = (CLHEP::RandExponential::shoot(engine)) / avCollisions;
0280
0281   // Without using CLHEP: approx random exponential by...
0282   //    double da = double(rand())/double(RAND_MAX);
0283   //    double step = -log(1.-da)/avCollisions;
0284
0285   LogTrace(me) << " step = " << step << " cm";
0286   // Next line only used to fill a container of 'step's for later diagnostic
0287   // dumps
0288   if (dumpGasCollisions())
0289     theCrossGap->addStep(step);
0290   return step;
0291 }
0292
0293 float CSCGasCollisions::generateEnergyLoss(double avCollisions,
0294                                            double anmin,
0295                                            double anmax,
0296                                            const std::vector<float> &collisions,
0297                                            CLHEP::HepRandomEngine *engine) const {
0298   // Generate a no. of collisions between collisions[0] and [N_ENERGY-1]
0299   float lnColl = log(CLHEP::RandFlat::shoot(engine, anmin, anmax));
0300
0301   // Without using CLHEP: approx random between anmin and anmax
0302   //    double ra = double(rand())/double(RAND_MAX)*avCollisions;
0303   //    cout << "ra = " << ra << std::endl;
0304   //    float lnColl = static_cast<float>( log( ra ) );
0305
0306   // Find energy loss for that number
0307   float lnE = lnEnergyLoss(lnColl, collisions);
0308   float eloss = exp(lnE);
0309   // Compensate if gamma was actually below 1.1
0310   if (theCrossGap->gamma() < 1.1)
0311     eloss = eloss * 0.173554 / theCrossGap->beta2();
0312   LogTrace(me) << "eloss = " << eloss << " eV";
0313   // Next line only used to fill container of eloss's for later diagnostic dumps
0314   if (dumpGasCollisions())
0315     theCrossGap->addEloss(eloss);
0316   return eloss;
0317 }
0318
0319 void CSCGasCollisions::ionize(double energyAvailable, LocalPoint startHere) const {
0320   while (energyAvailable > eion) {
0321     LogTrace(me) << "     NEW CLUSTER " << theCrossGap->noOfClusters() + 1 << " AT " << startHere;
0322     LocalPoint newCluster(startHere);
0323     theCrossGap->addCluster(newCluster);
0324
0325     //@@ I consider NOT subtracting eion before calculating range to be a bug.
0326     //@@ But this changes tuning of the algorithm so leave it until after the
0327     // big rush to 7_5_0
0328     //@@  energyAvailable -= eion;
0329
0330     // Sauli CERN 77-09: delta e range with E in MeV (Sauli references Kobetich
0331     // & Katz 1968 but that has nothing to do with this expression! He seems to
0332     // have  made a mistake.) I found the Sauli-quoted relationship in R.
0333     // Glocker, Z. Naturforsch. Ba, 129 (1948): delta e range R = aE^n with
0334     // a=710, n=1.72 for E in MeV and R in mg/cm^2 applicable over the range E =
0335     // 0.001 to 0.3 MeV.
0336
0337     // Take HALF that range. //@@ Why? Why not...
0338     double range = 0.5 * (0.71 / gasDensity) * pow(energyAvailable * 1.E-6, 1.72);
0339     LogTrace(me) << " range = " << range << " cm";
0340     if (range < clusterExtent) {
0341       // short-range delta e
0342       // How many electrons can we make? Now use *average* energy for ionization
0343       // (not *minimum*)
0344       int nelec = static_cast<int>(energyAvailable / ework);
0345       LogTrace(me) << "short-range delta energy in = " << energyAvailable << " eV";
0346       // energyAvailable -= nelec*(energyAvailable/ework);
0347       energyAvailable -= nelec * ework;
0348       // If still above eion (minimum, not average) add one more e
0349       if (energyAvailable > eion) {
0350         ++nelec;
0351         energyAvailable -= eion;
0352       }
0353       LogTrace(me) << "short-range delta energy out = " << energyAvailable << " eV, nelec = " << nelec;
0354       theCrossGap->addElectrons(nelec);
0355       break;
0356
0357     } else {
0358       // long-range delta e
0359       LogTrace(me) << "long-range delta \n"
0360                    << "no. of electrons in cluster now = " << theCrossGap->noOfElectrons();
0361       theCrossGap->addElectrons(1);  // Position is at startHere still
0362
0363       bool new_range = false;
0364       while (!new_range && (energyAvailable > ework)) {
0365         energyAvailable -= ework;
0366         while (energyAvailable > eion) {
0367           double range2 = 0.5 * 0.71 / gasDensity * pow(1.E-6 * energyAvailable, 1.72);
0368           double drange = range - range2;
0369           LogTrace(me) << "  energy left = " << energyAvailable << " eV, range2 = " << range2
0370                        << " cm, drange = " << drange << " cm";
0371           if (drange < clusterExtent) {
0372             theCrossGap->addElectronToBack();  // increment last element
0373           } else {
0374             startHere += drange * theCrossGap->unitVector();  // update delta e start position
0375             range = range2;                                   // update range
0376             new_range = true;                                 // Test range again
0377             LogTrace(me) << "reset range to range2 = " << range << " from startHere = " << startHere << "  and iterate";
0378           }
0379           break;  // out of inner while energyAvailable>eion
0380
0381         }  // inner while energyAvailable>eion
0382
0383       }  // while !new_range && energyAvailable>ework
0384
0385       // energyAvailable now less than ework, but still may be over eion...add
0386       // an e
0387       if (energyAvailable > eion) {
0388         energyAvailable -= ework;          // yes, it may go negative
0389         theCrossGap->addElectronToBack();  // add one more e
0390       }
0391
0392     }  // if range
0393
0394   }  // outer while energyAvailable>eion
0395 }
0396
0397 void CSCGasCollisions::writeSummary(int n_try, int n_steps, double sum_steps, float dedx, const PSimHit &simhit) const {
0398   // Switched from std::cout to LogVerbatim, Mar 2015
0399
0400   std::vector<LocalPoint> ion_clusters = theCrossGap->ionClusters();
0401   std::vector<int> electrons = theCrossGap->electrons();
0402   std::vector<float> elosses = theCrossGap->eLossPerStep();
0403   std::vector<double> steps = theCrossGap->stepLengths();
0404
0405   edm::LogVerbatim("CSCDigitizer") << "------------------"
0406                                    << "\nAFTER CROSSING GAP";
0407   /*
0408   edm::LogVerbatim("CSCDigitizer") << "no. of steps tried             = " <<
0409   n_try
0410                                  << "\nno. of steps from theCrossGap  = " <<
0411   theCrossGap->noOfSteps()
0412                                  << "\nsize of steps vector           = " <<
0413   steps.size();
0414
0415   edm::LogVerbatim("CSCDigitizer") << "no. of collisions (steps)      = " <<
0416   n_steps
0417                                  << "\nsize of elosses vector         = " <<
0418   elosses.size()
0419                                  << "\nsize of ion clusters vector    = " <<
0420   ion_clusters.size()
0421                                  << "\nsize of electrons vector       = " <<
0422   electrons.size();
0423   */
0424
0425   size_t nsteps = n_steps;          // force ridiculous type conversion
0426   size_t mstep = steps.size() - 1;  // final step gets filled but is outside gas
0427                                     // gap - unless we reach MAX_STEPS
0428   if ((nsteps != MAX_STEPS) && (nsteps != mstep)) {
0429     edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] WARNING! no. of steps = " << nsteps
0430                                      << " .ne. steps.size()-1 = " << mstep;
0431   }
0432   size_t meloss = elosses.size();
0433
0434   if (nsteps != meloss) {
0435     edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] WARNING! no. of steps = " << nsteps
0436                                      << " .ne. no. of elosses = " << meloss;
0437   } else {
0438     edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] # / length of step / energy loss per collision:";
0439     for (size_t i = 0; i != nsteps; ++i) {
0440       edm::LogVerbatim("CSCDigitizer") << i + 1 << " / S: " << steps[i] << " / E: " << elosses[i];
0441     }
0442   }
0443
0444   size_t mclus = ion_clusters.size();
0445   size_t melec = electrons.size();
0446
0447   if (mclus != melec) {
0448     edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] WARNING! size of cluster vector = " << mclus
0449                                      << " .ne. size of electrons vector = " << melec;
0450   } else {
0451     edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] # / postion of "
0452                                         "cluster / electrons per cluster: ";
0453     for (size_t i = 0; i != mclus; ++i) {
0454       edm::LogVerbatim("CSCDigitizer") << i + 1 << " / I: " << ion_clusters[i] << " / E: " << electrons[i];
0455     }
0456   }
0457
0458   int n_ic = count_if(elosses.begin(), elosses.end(), [&](auto c) { return c > this->eion; });
0459
0460   edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] total no. of collision steps = " << n_steps;
0461   edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] total sum of steps = " << sum_steps << " cm";
0462   if (nsteps > 0)
0463     edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] average step length = " << sum_steps / float(nsteps)
0464                                      << " cm";
0465   edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] total energy loss across gap = " << dedx
0466                                    << " eV = " << dedx / 1000. << " keV";
0467   edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] no. of primary ionizing collisions across gap = "
0468                                       "no. with eloss > eion = "
0469                                    << eion << " eV = " << n_ic;
0470   if (nsteps > 0)
0471     edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] average energy loss/collision = " << dedx / float(nsteps)
0472                                      << " eV";
0473
0474   std::vector<int>::const_iterator bigger = find(electrons.begin(), electrons.end(), 0);
0475   if (bigger != electrons.end()) {
0476     edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] TROUBLE! There is a cluster with 0 electrons.";
0477   }
0478
0479   int n_e = accumulate(electrons.begin(), electrons.end(), 0);
0480
0481   edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] SUMMARY: simhit"
0482                                    << "\n entry = " << simhit.entryPoint() << ": exit  = " << simhit.exitPoint()
0483                                    << "\n particle type = " << simhit.particleType()
0484                                    << " : momentum = " << simhit.pabs()
0485                                    << " GeV/c : energy loss = " << simhit.energyLoss() * 1.E06 << " keV";
0486
0487   if (nsteps > 0) {
0488     edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] SUMMARY: ionization"
0489                                      << " : steps= " << nsteps << " : sum(steps)= " << sum_steps
0490                                      << " cm : <step>= " << sum_steps / float(nsteps) << " cm"
0491                                      << " : ionizing= " << n_ic << " : ionclus= " << mclus << " : total e= " << n_e
0492                                      << " : <dedx/step>= " << dedx / float(nsteps)
0493                                      << " eV : <e/ionclus>= " << float(n_e) / float(mclus)
0494                                      << " : dedx= " << dedx / 1000. << " keV";
0495   } else {
0496     edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisons] ERROR? no collision steps!";
0497   }
0498
0499   // Turn off output file -  used for initial development
0500   //  if ( saveGasCollisions() ) {
0501   //    ofstream of0("osteplen.dat",ios::app);
0502   //    std::copy( steps.begin(), steps.end(),
0503   //    std::ostream_iterator<float>(of0,"\n"));
0504
0505   //    ofstream of1("olperc.dat",ios::app);
0506   //    std::copy( elosses.begin(), elosses.end(),
0507   //    std::ostream_iterator<float>(of1,"\n"));
0508
0509   //   ofstream of2("oclpos.dat",ios::app);
0510   //   std::copy( ion_clusters.begin(), ion_clusters.end(),
0511   //   std::ostream_iterator<LocalPoint>(of2,"\n"));
0512
0513   //   ofstream of3("oepercl.dat",ios::app);
0514   //   std::copy( electrons.begin(), electrons.end(),
0515   //   std::ostream_iterator<int>(of3,"\n"));
0516   // }
0517 }
0518
0519 float CSCGasCollisions::lnEnergyLoss(float lnCollisions, const std::vector<float> &collisions) const {
0520   float lnE = -1.;
0521
0522   // Find collision[] bin in which lnCollisions falls
0523   std::vector<float>::const_iterator it = find(collisions.begin(), collisions.end(), lnCollisions);
0524
0525   if (it != collisions.end()) {
0526     // found the value
0527     std::vector<float>::difference_type ihi = it - collisions.begin();
0528     LogTrace(me) << ": using one energy bin " << ihi << " = " << theEnergyBins[ihi]
0529                  << " for lnCollisions = " << lnCollisions;
0530     lnE = theEnergyBins[ihi];
0531   } else {
0532     // interpolate the value
0533     std::vector<float>::const_iterator loside =
0534         find_if(collisions.begin(), collisions.end(), [&lnCollisions](auto c) { return c < lnCollisions; });
0535     std::vector<float>::difference_type ilo = loside - collisions.begin();
0536     if (ilo > 0) {
0537       LogTrace(me) << ": using energy bin " << ilo - 1 << " and " << ilo;
0538       lnE = theEnergyBins[ilo - 1] + (lnCollisions - collisions[ilo - 1]) *
0539                                          (theEnergyBins[ilo] - theEnergyBins[ilo - 1]) /
0540                                          (collisions[ilo] - collisions[ilo - 1]);
0541     } else {
0542       LogTrace(me) << ": using one energy bin 0 = " << theEnergyBins[0] << " for lnCollisions = " << lnCollisions;
0543       lnE = theEnergyBins[0];  //@@ WHAT ELSE TO DO?
0544     }
0545   }
0546
0547   return lnE;
0548 }
0549
0550 void CSCGasCollisions::fillCollisionsForThisGamma(float logGamma, std::vector<float> &collisions) const {
0551   std::vector<float>::const_iterator bigger =
0552       find_if(theGammaBins.begin(), theGammaBins.end(), [&logGamma](auto c) { return c > logGamma; });
0553
0554   if (bigger == theGammaBins.end()) {
0555     // use highest bin
0556     LogTrace(me) << ": using highest gamma bin"
0557                  << " for logGamma = " << logGamma;
0558     for (int i = 0; i < N_ENERGY; ++i)
0559       collisions[i] = theCollisionTable[i * N_GAMMA];
0560   } else {
0561     // use bigger and its lower neighbour
0562     std::vector<float>::difference_type ihi = bigger - theGammaBins.begin();
0563     if (ihi > 0) {
0564       double dlg2 = *bigger--;  // and decrement after deref
0565       // LogTrace(me) << ": using gamma bins "
0566       //                   << ihi-1 << " and " << ihi;
0567       double dlg1 = *bigger;  // now the preceding element
0568       double dlg = (logGamma - dlg1) / (dlg2 - dlg1);
0569       double omdlg = 1. - dlg;
0570       for (int i = 0; i < N_ENERGY; ++i)
0571         collisions[i] = theCollisionTable[i * N_GAMMA + ihi - 1] * omdlg + theCollisionTable[i * N_GAMMA + ihi] * dlg;
0572     } else {
0573       // bigger has no lower neighbour
0574       LogTrace(me) << ": using lowest gamma bin"
0575                    << " for logGamma = " << logGamma;
0576
0577       for (int i = 0; i < N_ENERGY; ++i)
0578         collisions[i] = theCollisionTable[i * N_GAMMA];
0579     }
0580   }
0581 }