#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().

38 { return dumpGasCollisions_; }

CSCGasCollisions::dumpGasCollisions_

bool dumpGasCollisions_

Definition: CSCGasCollisions.h:78

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

private

Definition at line 545 of file CSCGasCollisions.cc.

References 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(), bind2nd(greater<float>(), 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 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(), bind2nd(less<float>(), 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().

39 { return saveGasCollisions_; }

CSCGasCollisions::saveGasCollisions_

bool saveGasCollisions_

Definition: CSCGasCollisions.h:77

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 eion, CSCCrossGap::electrons(), electrons_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(), bind2nd(greater<float>(), 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"));
   // }
 
 }