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
CSCCrossGap *	theCrossGap
std::vector< float >	theEnergyBins
std::vector< float >	theGammaBins
const ParticleDataTable *	theParticleDataTable

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

(

const edm::ParameterSet &

pset

)

Definition at line 54 of file CSCGasCollisions.cc.

References clusterExtent, deCut, dumpGasCollisions(), dumpGasCollisions_, eion, ework, gasDensity, edm::ParameterSet::getUntrackedParameter(), me, readCollisionTable(), and saveGasCollisions().

 : me("CSCGasCollisions"),
  gasDensity( 2.1416e-03 ), 
  deCut( 1.e05 ), eion( 14.95 ), ework( 34.0 ), clusterExtent( 0.001 ),
  theGammaBins(N_GAMMA, 0.), theEnergyBins(N_ENERGY, 0.), 
  theCollisionTable(N_ENTRIES, 0.), theCrossGap( nullptr ),
  theParticleDataTable(nullptr),
  saveGasCollisions_ ( false ), dumpGasCollisions_( false )
{

  dumpGasCollisions_ = pset.getUntrackedParameter<bool>("dumpGasCollisions");

  edm::LogInfo(me) << "Constructing a " << me << ":"; 
  edm::LogInfo(me) << "gas density = " << gasDensity << " g/cm3";
  edm::LogInfo(me) << "max eloss per collision allowed = " << deCut/1000. 
       << " keV (for higher elosses, hits should have been simulated.)";
  edm::LogInfo(me) << "ionization threshold = " << eion << " eV";
  edm::LogInfo(me) << "effective work function = " << ework << " eV";
  edm::LogInfo(me) << "cluster extent = " << clusterExtent*1.e04 << " micrometres";
  edm::LogInfo(me) << "dump gas collision and simhit information? " << dumpGasCollisions();
  edm::LogInfo(me) << "save gas collision information? -NOT YET IMPLEMENTED- " << saveGasCollisions();

  readCollisionTable();
}

CSCGasCollisions::~CSCGasCollisions ( )

virtual

Definition at line 79 of file CSCGasCollisions.cc.

References me, and theCrossGap.

                                    {
  edm::LogInfo(me) << "Destructing a " << me;
  delete theCrossGap;
}

Member Function Documentation

bool CSCGasCollisions::dumpGasCollisions ( void  ) const

inline

Definition at line 38 of file CSCGasCollisions.h.

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

{ return dumpGasCollisions_; }

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

private

Definition at line 545 of file CSCGasCollisions.cc.

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

Referenced by simulate().

{
  std::vector<float>::const_iterator bigger = find_if(theGammaBins.begin(),
    theGammaBins.end(), [&logGamma](auto c){return c > logGamma;});

  if ( bigger == theGammaBins.end() ) {
    // use highest bin
    LogTrace(me) << ": using highest gamma bin" 
              << " for logGamma = " << logGamma;
    for (int i=0; i<N_ENERGY; ++i)
      collisions[i] = theCollisionTable[i*N_GAMMA];
  }
  else {
    // use bigger and its lower neighbour
    std::vector<float>::difference_type ihi = bigger - theGammaBins.begin();
    if ( ihi > 0 ) {
      double dlg2 = *bigger--; // and decrement after deref
      //LogTrace(me) << ": using gamma bins " 
      //                   << ihi-1 << " and " << ihi;
      double dlg1 = *bigger;   // now the preceding element
      double dlg = (logGamma-dlg1)/(dlg2-dlg1);
      double omdlg = 1. - dlg;
      for (int i=0; i<N_ENERGY; ++i)
        collisions[i] = theCollisionTable[i*N_GAMMA+ihi-1]*omdlg +
          theCollisionTable[i*N_GAMMA+ihi]*dlg;
    }
    else {
      // bigger has no lower neighbour
       LogTrace(me) << ": using lowest gamma bin" 
                     << " for logGamma = " << logGamma;

      for (int i=0; i<N_ENERGY; ++i)
        collisions[i] = theCollisionTable[i*N_GAMMA];
    }
  }
}

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

private

Definition at line 291 of file CSCGasCollisions.cc.

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

Referenced by simulate().

{
// Generate a no. of collisions between collisions[0] and [N_ENERGY-1]
   float lnColl = log(CLHEP::RandFlat::shoot(engine, anmin, anmax));

// Without using CLHEP: approx random between anmin and anmax
//    double ra = double(rand())/double(RAND_MAX)*avCollisions;
//    cout << "ra = " << ra << std::endl;
//    float lnColl = static_cast<float>( log( ra ) );

    // Find energy loss for that number
    float lnE    = lnEnergyLoss( lnColl, collisions );
    float eloss  = exp(lnE);
    // Compensate if gamma was actually below 1.1
    if ( theCrossGap->gamma() < 1.1 ) eloss = eloss * 0.173554/theCrossGap->beta2();
    LogTrace(me) << "eloss = " << eloss << " eV";
    // Next line only used to fill container of eloss's for later diagnostic dumps
    if ( dumpGasCollisions() ) theCrossGap->addEloss( eloss );
    return eloss;
}

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

private

Definition at line 276 of file CSCGasCollisions.cc.

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

Referenced by simulate().

{
// Generate a m.f.p.  (1/avCollisions = cm/collision)
  double step = (CLHEP::RandExponential::shoot(engine))/avCollisions;

// Without using CLHEP: approx random exponential by...
//    double da = double(rand())/double(RAND_MAX);
//    double step = -log(1.-da)/avCollisions;

  LogTrace(me)  << " step = " << step << " cm";
    // Next line only used to fill a container of 'step's for later diagnostic dumps
  if ( dumpGasCollisions() ) theCrossGap->addStep( step );
    return step;
}

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

private

Definition at line 313 of file CSCGasCollisions.cc.

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

Referenced by simulate().

{
  while ( energyAvailable > eion ) {
    LogTrace(me) << "     NEW CLUSTER " << theCrossGap->noOfClusters() + 1 <<
           " AT " << startHere;
    LocalPoint newCluster( startHere );
    theCrossGap->addCluster(newCluster);

  //@@ I consider NOT subtracting eion before calculating range to be a bug.
  //@@ But this changes tuning of the algorithm so leave it until after the big rush to 7_5_0
  //@@  energyAvailable -= eion; 

  // Sauli CERN 77-09: delta e range with E in MeV (Sauli references Kobetich & Katz 1968
  // but that has nothing to do with this expression! He seems to have  made a mistake.)
    // I found the Sauli-quoted relationship in R. Glocker, Z. Naturforsch. Ba, 129 (1948):
    // delta e range R = aE^n with a=710, n=1.72 for E in MeV and R in mg/cm^2
    // applicable over the range E = 0.001 to 0.3 MeV.

  // Take HALF that range. //@@ Why? Why not...
    double range = 0.5 * (0.71/gasDensity)*pow( energyAvailable*1.E-6, 1.72);
    LogTrace(me) << " range = " << range << " cm";
    if ( range < clusterExtent ) {

      // short-range delta e
          // How many electrons can we make? Now use *average* energy for ionization (not *minimum*)
      int nelec = static_cast<int>(energyAvailable/ework);
      LogTrace(me) << "short-range delta energy in = " << energyAvailable << " eV";
      //energyAvailable -= nelec*(energyAvailable/ework);
      energyAvailable -= nelec*ework;
        // If still above eion (minimum, not average) add one more e
      if ( energyAvailable > eion ) {
        ++nelec;
        energyAvailable -= eion;
      }
      LogTrace(me) << "short-range delta energy out = " << energyAvailable << " eV, nelec = " << nelec;
      theCrossGap->addElectrons( nelec );
      break;

     }
     else {
      // long-range delta e
         LogTrace(me) << "long-range delta \n"
             << "no. of electrons in cluster now = " << theCrossGap->noOfElectrons();
         theCrossGap->addElectrons( 1 ); // Position is at startHere still

          bool new_range = false;
          while ( !new_range && (energyAvailable>ework) ) {
            energyAvailable -= ework;
            while ( energyAvailable > eion ) {
              double range2 = 0.5 * 0.71/gasDensity*pow( 1.E-6*energyAvailable, 1.72);
              double drange = range - range2;
              LogTrace(me) << "  energy left = " << energyAvailable << 
        " eV, range2 = " << range2 << " cm, drange = " << drange << " cm";
              if ( drange < clusterExtent ) {
                theCrossGap->addElectronToBack(); // increment last element
              }
              else {
                startHere += drange*theCrossGap->unitVector(); // update delta e start position
                range = range2;                       // update range
                new_range = true;                     // Test range again
                LogTrace(me) << "reset range to range2 = " << range 
                             << " from startHere = " << startHere << "  and iterate";
              }
              break; // out of inner while energyAvailable>eion

            } // inner while energyAvailable>eion

          } // while !new_range && energyAvailable>ework

      // energyAvailable now less than ework, but still may be over eion...add an e
          if ( energyAvailable > eion ) {
            energyAvailable -= ework; // yes, it may go negative
            theCrossGap->addElectronToBack(); // add one more e
          }

        } // if range

      } // outer while energyAvailable>eion
}

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

private

Definition at line 505 of file CSCGasCollisions.cc.

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

Referenced by generateEnergyLoss().

                                              {

  float lnE = -1.;

   // Find collision[] bin in which lnCollisions falls
    std::vector<float>::const_iterator it = find(collisions.begin(),
      collisions.end(), lnCollisions );

    if ( it != collisions.end() ) {
      // found the value
      std::vector<float>::difference_type ihi = it - collisions.begin();
      LogTrace(me) << ": using one energy bin " << ihi << " = " 
                         << theEnergyBins[ihi]
                 << " for lnCollisions = " << lnCollisions;
      lnE = theEnergyBins[ihi];
    }
    else {
      // interpolate the value
      std::vector<float>::const_iterator loside = find_if(collisions.begin(),
        collisions.end(), [&lnCollisions](auto c){return c < lnCollisions;});
      std::vector<float>::difference_type ilo = loside - collisions.begin();
      if ( ilo > 0 ) {
        LogTrace(me) << ": using energy bin " 
                           << ilo-1 << " and " << ilo;
        lnE = theEnergyBins[ilo-1] + (lnCollisions-collisions[ilo-1])*
          (theEnergyBins[ilo]-theEnergyBins[ilo-1]) /
             (collisions[ilo]-collisions[ilo-1]);
      }
      else {
        LogTrace(me) << ": using one energy bin 0 = " 
                           << theEnergyBins[0]
               << " for lnCollisions = " << lnCollisions;
        lnE = theEnergyBins[0]; //@@ WHAT ELSE TO DO?
      }
    }

    return lnE;
}

void CSCGasCollisions::readCollisionTable ( )

private

Definition at line 84 of file CSCGasCollisions.cc.

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

Referenced by CSCGasCollisions().

                                          {

  // I'd prefer to allow comments in the data file which means
  // complications of reading line-by-line and then item-by-item.

  // Is float OK? Or do I need double?

  // Do I need any sort of error trapping?

  // We use the default CMSSW data file path
  // We use the default file name 'collisions.dat'
   
  // This can be reset in .orcarc by SimpleConfigurable
  //    Muon:Endcap:CollisionsFile

  // TODO make configurable
  string colliFile = "SimMuon/CSCDigitizer/data/collisions.dat";
  edm::FileInPath f1{ colliFile };
  edm::LogInfo(me) << ": reading " << f1.fullPath();

  ifstream  fin{ f1.fullPath() };

  if (fin.fail()) {
    string errorMessage = "Cannot open input file " + f1.fullPath();
    edm::LogError("CSCGasCollisions") << errorMessage;
    throw cms::Exception(errorMessage);
  }

  fin.clear( );                  // Clear eof read status
  fin.seekg( 0, ios::beg );      // Position at start of file

  // @@ We had better have the right sizes everywhere or all
  // hell will break loose. There's no trapping.
  LogTrace(me) << "Reading gamma bins";
  int j = 0;
  for(int i = 0; i<N_GAMMA; ++i) {
    fin >> theGammaBins[i];
    LogTrace(me) << ++j << "   " << theGammaBins[i];
  }

  j = 0;
  LogTrace(me) << "Reading energy bins \n";
  for(int i = 0; i<N_ENERGY; ++i) {
     fin >> theEnergyBins[i];
     LogTrace(me) << ++j << "   " << theEnergyBins[i];
  }

  j = 0;
  LogTrace(me) << "Reading collisions table \n";

  for(int i = 0; i<N_ENTRIES; ++i) {
     fin >> theCollisionTable[i];
     LogTrace(me) << ++j << "   " << theCollisionTable[i];
  }

  fin.close();
}

bool CSCGasCollisions::saveGasCollisions ( void  ) const

inline

Definition at line 39 of file CSCGasCollisions.h.

Referenced by CSCGasCollisions().

{ return saveGasCollisions_; }

void CSCGasCollisions::setParticleDataTable

(

const ParticleDataTable *

pdt

)

Definition at line 143 of file CSCGasCollisions.cc.

References theParticleDataTable.

Referenced by CSCWireHitSim::setParticleDataTable().

{
  theParticleDataTable = pdt;
}

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

Definition at line 149 of file CSCGasCollisions.cc.

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

Referenced by CSCWireHitSim::getIonizationClusters().

                                                                                               {

  const float epsilonL = 0.01;                     // Shortness of simhit 'length'
  //  const float max_gap_z = 1.5;                     // Gas gaps are 0.5 or 1.0 cm

  // Note that what I call the 'gap' may in fact be the 'length' of a PSimHit which
  // does not start and end on the gap edges. This confuses the nomenclature at least.

  double mom    = simhit.pabs();                   // in GeV/c - see MuonSensitiveDetector.cc
  //  int iam       = simhit.particleType();           // PDG type
  delete theCrossGap;                              // before building new one
  assert(theParticleDataTable != nullptr);
  ParticleData const * particle = theParticleDataTable->particle( simhit.particleType() );
  double mass = 0.105658; // assume a muon
  if(particle == nullptr)
  {
     edm::LogError("CSCGasCollisions") << "Cannot find particle of type " << simhit.particleType()
            << " in the PDT";
  }
  else 
  {
    mass = particle->mass();
  }

  theCrossGap   = new CSCCrossGap( mass, mom, simhit.exitPoint() - simhit.entryPoint() );
  float gapSize = theCrossGap->length();

  // Test the simhit 'length' (beware of angular effects)
  //  if ( gapSize <= epsilonL || gapSize > max_gap_z ) {
  if ( gapSize <= epsilonL ) {
    edm::LogVerbatim("CSCDigitizer") 
      << "[CSCGasCollisions] WARNING! simhit entry and exit are too close - skipping simhit:"
      << "\n entry = " << simhit.entryPoint()
      << ": exit  = " << simhit.exitPoint()
      << "\n particle type = " << simhit.particleType() << " : momentum = " << simhit.pabs()
      << " GeV/c : energy loss = " << simhit.energyLoss()*1.E06 << " keV"
      << ", gapSize = " << gapSize << " cm (< epsilonL = " << epsilonL << " cm)";
    return; //@@ Just skip this PSimHit
  }

  // Interpolate the table for current gamma value
  // Extract collisions binned by energy loss values, for this gamma
  std::vector<float> collisions(N_ENERGY);
  double loggam = theCrossGap->logGamma(); 
  fillCollisionsForThisGamma( static_cast<float>(loggam), collisions );

  double anmin = exp(collisions[N_ENERGY-1]);
  double anmax = exp(collisions[0]);
  double amu   = anmax - anmin;

  LogTrace(me) << "collisions extremes = " << collisions[N_ENERGY-1]
     << ", " << collisions[0] << "\n"
     << "anmin = " << anmin << ", anmax = " << anmax << "\n"
     << "amu = " << amu << "\n";

  float dedx       = 0.; // total energy loss
  double sum_steps = 0.; // total distance across gap (along simhit direction)
  int n_steps      = 0;  // no. of steps/primary collisions
  int n_try        = 0;  // no. of tries to generate steps
  double step      = -1.; // Sentinel for start

  LocalPoint layerLocalPoint( simhit.entryPoint() );

  // step/primary collision loop
  while ( sum_steps < gapSize) {
    ++n_try;
    if ( n_try > MAX_STEPS ) {
      int maxst = MAX_STEPS;
      edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] WARNING! n_try = " << n_try 
       << " is > MAX_STEPS = " << maxst << " - skipping simhit:"
           << "\n entry = " << simhit.entryPoint()
           << ": exit  = " << simhit.exitPoint()
           << "\n particle type = " << simhit.particleType() << " : momentum = " << simhit.pabs()
           << " GeV/c : energy loss = " << simhit.energyLoss()*1.E06 << " keV"
           << "\n gapSize = " << gapSize << " cm, last step = " << step 
           << " cm, sum_steps = " << sum_steps << " cm, n_steps = " << n_steps;
      break;
    }
    step = generateStep( amu, engine );
    if ( sum_steps + step > gapSize ) break; // this step goes too far

    float eloss = generateEnergyLoss( amu, anmin, anmax, collisions, engine );

    // Is the eloss too large? (then GEANT should have produced hits!)
    if ( eloss > deCut ) {
      edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] WARNING! eloss = " << eloss 
              << " eV is too large (> " << deCut << " eV) - trying another collision";
      continue; // to generate another collision/step
    }

    dedx      += eloss; // the energy lost from the ionizing particle
    sum_steps += step;  // the position of the ionizing particle
    ++n_steps;          // the number of primary collisions
    
    if (n_steps > MAX_STEPS ) { // Extra-careful trap for bizarreness
        edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] WARNING! " << n_steps 
             << " is too many steps -  skipping simhit:"
             << "\n entry = " << simhit.entryPoint()
             << ": exit  = " << simhit.exitPoint()
             << "\n particle type = " << simhit.particleType() << " : momentum = " << simhit.pabs()
             << " GeV/c : energy loss = " << simhit.energyLoss()*1.E06 << " keV"
         << "\ngapSize=" << gapSize << " cm, last step=" << step << " cm, sum_steps=" << sum_steps << " cm";
        break;
    }
    LogTrace(me) << "sum_steps = " << sum_steps << " cm , dedx = " << dedx << " eV";
    
    // Generate ionization. 
    // eion is the minimum energy at which ionization can occur in the gas
    if ( eloss > eion  ) {
      layerLocalPoint += step * theCrossGap->unitVector(); // local point where the collision occurs
      ionize( eloss, layerLocalPoint );
    }
    else {
      LogTrace(me) << "Energy available = " << eloss << " eV is too low for ionization (< eion = " << eion << " eV)";
    }

  } // step/collision loop

  if ( dumpGasCollisions() ) writeSummary( n_try, n_steps, sum_steps, dedx, simhit );

  // Return values in two container arguments
  positions = theCrossGap->ionClusters();
  electrons = theCrossGap->electrons();
  return;
}

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

private

Definition at line 393 of file CSCGasCollisions.cc.

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

Referenced by simulate().

{
  // Switched from std::cout to LogVerbatim, Mar 2015

  std::vector<LocalPoint> ion_clusters   = theCrossGap->ionClusters();
  std::vector<int> electrons             = theCrossGap->electrons();
  std::vector<float> elosses             = theCrossGap->eLossPerStep();
  std::vector<double> steps              = theCrossGap->stepLengths();

  edm::LogVerbatim("CSCDigitizer") << "------------------"
                                 << "\nAFTER CROSSING GAP";
  /*
  edm::LogVerbatim("CSCDigitizer") << "no. of steps tried             = " << n_try
                                 << "\nno. of steps from theCrossGap  = " << theCrossGap->noOfSteps()
                 << "\nsize of steps vector           = " << steps.size();
                                 
  edm::LogVerbatim("CSCDigitizer") << "no. of collisions (steps)      = " << n_steps
                     << "\nsize of elosses vector         = " << elosses.size()
                 << "\nsize of ion clusters vector    = " << ion_clusters.size()
                                 << "\nsize of electrons vector       = " << electrons.size();
  */

  size_t nsteps = n_steps; // force ridiculous type conversion
  size_t mstep = steps.size() - 1; // final step gets filled but is outside gas gap - unless we reach MAX_STEPS
  if ( (nsteps != MAX_STEPS) && (nsteps != mstep) ) {
    edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] WARNING! no. of steps = " << nsteps 
                                     << " .ne. steps.size()-1 = " << mstep;
  }
  size_t meloss = elosses.size();

  if ( nsteps != meloss ){
    edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] WARNING! no. of steps = " << nsteps 
                                     << " .ne. no. of elosses = " << meloss;
  }
  else {
    edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] # / length of step / energy loss per collision:";
    for ( size_t i = 0; i != nsteps; ++i ) {
      edm::LogVerbatim("CSCDigitizer") << i+1 << " / S: " << steps[i] << " / E: " << elosses[i];
    }
  }

  size_t mclus = ion_clusters.size();
  size_t melec = electrons.size();

  if ( mclus != melec ){
    edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] WARNING! size of cluster vector = " << mclus 
            << " .ne. size of electrons vector = " << melec;
  }
  else {
    edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] # / postion of cluster / electrons per cluster: ";
    for ( size_t i = 0; i != mclus; ++i ) {
      edm::LogVerbatim("CSCDigitizer") << i+1 <<  " / I: " << ion_clusters[i] << " / E: " << electrons[i];
    }
  }

  int n_ic = count_if( elosses.begin(), elosses.end(), [&](auto c){return c > this->eion;} );

  edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] total no. of collision steps = " << n_steps;
  edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] total sum of steps = " << sum_steps << " cm";
  if ( nsteps > 0 ) edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] average step length = " << 
      sum_steps/float(nsteps) << " cm";
  edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] total energy loss across gap = " << 
      dedx << " eV = " << dedx/1000. << " keV";
  edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] no. of primary ionizing collisions across gap = no. with eloss > eion = "
      << eion << " eV = " << n_ic;
  if ( nsteps > 0 ) edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] average energy loss/collision = " <<  
      dedx/float(nsteps) << " eV";

  std::vector<int>::const_iterator bigger = find(electrons.begin(), electrons.end(), 0 );
  if ( bigger != electrons.end() ) {
    edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] TROUBLE! There is a cluster with 0 electrons.";
  }

  int n_e = accumulate(electrons.begin(), electrons.end(), 0 );

  edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] SUMMARY: simhit"
      << "\n entry = " << simhit.entryPoint()
      << ": exit  = " << simhit.exitPoint()
      << "\n particle type = " << simhit.particleType() << " : momentum = " << simhit.pabs()
      << " GeV/c : energy loss = " << simhit.energyLoss()*1.E06 << " keV";

  if ( nsteps > 0 ) {
    edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisions] SUMMARY: ionization" 
      << " : steps= " << nsteps << " : sum(steps)= " << sum_steps << " cm : <step>= " 
      << sum_steps/float(nsteps) << " cm"
      << " : ionizing= " << n_ic << " : ionclus= " << mclus << " : total e= " << n_e
      << " : <dedx/step>= " << dedx/float(nsteps)
      << " eV : <e/ionclus>= " << float(n_e)/float(mclus)
      << " : dedx= " << dedx/1000. << " keV";
  }
  else {
    edm::LogVerbatim("CSCDigitizer") << "[CSCGasCollisons] ERROR? no collision steps!";
  }

  // Turn off output file -  used for initial development
  //  if ( saveGasCollisions() ) {
  //    ofstream of0("osteplen.dat",ios::app);
  //    std::copy( steps.begin(), steps.end(), std::ostream_iterator<float>(of0,"\n"));

  //    ofstream of1("olperc.dat",ios::app);
  //    std::copy( elosses.begin(), elosses.end(), std::ostream_iterator<float>(of1,"\n"));

  //   ofstream of2("oclpos.dat",ios::app);
  //   std::copy( ion_clusters.begin(), ion_clusters.end(), std::ostream_iterator<LocalPoint>(of2,"\n"));

  //   ofstream of3("oepercl.dat",ios::app);
  //   std::copy( electrons.begin(), electrons.end(), std::ostream_iterator<int>(of3,"\n"));
  // }

}

Member Data Documentation

double CSCGasCollisions::clusterExtent

private

Definition at line 69 of file CSCGasCollisions.h.

Referenced by CSCGasCollisions(), and ionize().

double CSCGasCollisions::deCut

private

Definition at line 66 of file CSCGasCollisions.h.

Referenced by CSCGasCollisions(), and simulate().

bool CSCGasCollisions::dumpGasCollisions_

private

Definition at line 78 of file CSCGasCollisions.h.

Referenced by CSCGasCollisions().

double CSCGasCollisions::eion

private

Definition at line 67 of file CSCGasCollisions.h.

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

double CSCGasCollisions::ework

private

Definition at line 68 of file CSCGasCollisions.h.

Referenced by CSCGasCollisions(), and ionize().

double CSCGasCollisions::gasDensity

private

Definition at line 60 of file CSCGasCollisions.h.

Referenced by CSCGasCollisions(), and ionize().

const int CSCGasCollisions::MAX_STEPS = 400

static

Definition at line 44 of file CSCGasCollisions.h.

Referenced by simulate(), and writeSummary().

const std::string CSCGasCollisions::me

private

Definition at line 59 of file CSCGasCollisions.h.

Referenced by CSCGasCollisions(), fillCollisionsForThisGamma(), generateEnergyLoss(), generateStep(), ionize(), lnEnergyLoss(), readCollisionTable(), simulate(), and ~CSCGasCollisions().

const int CSCGasCollisions::N_ENERGY = 63

static

Definition at line 42 of file CSCGasCollisions.h.

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

const int CSCGasCollisions::N_ENTRIES = N_GAMMA*N_ENERGY

static

Definition at line 43 of file CSCGasCollisions.h.

Referenced by readCollisionTable().

const int CSCGasCollisions::N_GAMMA = 21

static

Definition at line 41 of file CSCGasCollisions.h.

Referenced by fillCollisionsForThisGamma(), and readCollisionTable().

bool CSCGasCollisions::saveGasCollisions_

private

Definition at line 77 of file CSCGasCollisions.h.

std::vector<float> CSCGasCollisions::theCollisionTable

private

Definition at line 73 of file CSCGasCollisions.h.

Referenced by fillCollisionsForThisGamma(), and readCollisionTable().

CSCCrossGap* CSCGasCollisions::theCrossGap

private

Definition at line 75 of file CSCGasCollisions.h.

Referenced by generateEnergyLoss(), generateStep(), ionize(), simulate(), writeSummary(), and ~CSCGasCollisions().

std::vector<float> CSCGasCollisions::theEnergyBins

private

Definition at line 72 of file CSCGasCollisions.h.

Referenced by lnEnergyLoss(), and readCollisionTable().

std::vector<float> CSCGasCollisions::theGammaBins

private

Definition at line 71 of file CSCGasCollisions.h.

Referenced by fillCollisionsForThisGamma(), and readCollisionTable().

const ParticleDataTable* CSCGasCollisions::theParticleDataTable

private

Definition at line 76 of file CSCGasCollisions.h.

Referenced by setParticleDataTable(), and simulate().