#include <SiLinearChargeDivider.h>

Inheritance diagram for SiLinearChargeDivider:

Public Member Functions
SiChargeDivider::ionization_type	divide (const PSimHit *, const LocalVector &, double, const StripGeomDetUnit &det)
void	setParticleDataTable (const ParticleDataTable *pdt)
	SiLinearChargeDivider (const edm::ParameterSet &conf, CLHEP::HepRandomEngine &)
virtual	~SiLinearChargeDivider ()
Private Member Functions
float	DeconvolutionShape (const PSimHit *, const StripGeomDetUnit &det)
float	driftXPos (const Local3DPoint &pos, const LocalVector &drift, double thickness)
void	fluctuateEloss (int particleId, float momentum, float eloss, float length, int NumberOfSegmentation, float elossVector[])
float	PeakShape (const PSimHit *, const StripGeomDetUnit &det)
float	TimeResponse (const PSimHit *hit, const StripGeomDetUnit &det)
Private Attributes
const int	chargedivisionsPerStrip
const double	cosmicShift
const double	deltaCut
std::unique_ptr < SiG4UniversalFluctuation >	fluctuate
const bool	fluctuateCharge
const bool	peakMode
CLHEP::HepRandomEngine &	rndEngine
const ParticleDataTable *	theParticleDataTable
Static Private Attributes
static float	decoValues [651]
static float	peakValues [921]

Detailed Description

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

fluctuations of the charge deposit fronm segment to segment
pulse shape ( in peak of deconvolution mode)

Definition at line 27 of file SiLinearChargeDivider.h.

Constructor & Destructor Documentation

SiLinearChargeDivider::SiLinearChargeDivider	(	const edm::ParameterSet &	conf,
		CLHEP::HepRandomEngine &	eng
	)

Definition at line 7 of file SiLinearChargeDivider.cc.

                                                                                                   :
  // Run APV in peak instead of deconvolution mode, which degrades the time resolution.
  peakMode(conf.getParameter<bool>("APVpeakmode")),
  // Enable interstrip Landau fluctuations within a cluster.
  fluctuateCharge(conf.getParameter<bool>("LandauFluctuations")),
  // Number of segments per strip into which charge is divided during
  // simulation. If large, precision of simulation improves.
  chargedivisionsPerStrip(conf.getParameter<int>("chargeDivisionsPerStrip")),
  // delta cutoff in MeV, has to be same as in Geant (0.120425 MeV corresponding to 100um range for electrons)
  deltaCut(conf.getParameter<double>("DeltaProductionCut")),
  //Offset for digitization during the MTCC and in general for taking cosmic particle
  //The value to be used it must be evaluated and depend on the volume defnition used
  //for the cosimc generation (Considering only the tracker the value is 11 ns)
  cosmicShift(conf.getUntrackedParameter<double>("CosmicDelayShift")),
  theParticleDataTable(0),
  rndEngine(eng),
  // Geant4 engine used to fluctuate the charge from segment to segment
  fluctuate(new SiG4UniversalFluctuation(rndEngine)) {
}

SiLinearChargeDivider::~SiLinearChargeDivider ( ) [virtual]

Definition at line 27 of file SiLinearChargeDivider.cc.

{
}

Member Function Documentation

float SiLinearChargeDivider::DeconvolutionShape	(	const PSimHit *	hit,
		const StripGeomDetUnit &	det
	)		`[private]`

Definition at line 122 of file SiLinearChargeDivider.cc.

References cosmicShift, decoValues, PSimHit::energyLoss(), PSimHit::localPosition(), mag(), GeomDet::surface(), PSimHit::tof(), Surface::toGlobal(), and x.

Referenced by TimeResponse().

                                                                                              {
  // x is difference between the tof and the tof for a photon (reference)
  // converted into a bin number
  int x = int(((det.surface().toGlobal(hit->localPosition()).mag()/30.) + cosmicShift - hit->tof())*10)+300;
  if(x < 0 || x > 650) return 0;
  return hit->energyLoss()*decoValues[x];
}

SiChargeDivider::ionization_type SiLinearChargeDivider::divide	(	const PSimHit *	hit,
		const LocalVector &	driftdir,
		double	moduleThickness,
		const StripGeomDetUnit &	det
	)		`[virtual]`

Implements SiChargeDivider.

Definition at line 31 of file SiLinearChargeDivider.cc.

References chargedivisionsPerStrip, driftXPos(), PSimHit::energyLoss(), PSimHit::entryPoint(), PSimHit::exitPoint(), fluctuateCharge, fluctuateEloss(), i, StripTopology::localPitch(), PSimHit::localPosition(), PV3DBase< T, PVType, FrameType >::mag(), PSimHit::pabs(), PSimHit::particleType(), StripGeomDetUnit::specificTopology(), and TimeResponse().

                                                                                                                                  {

  // computes the number of segments from number of segments per strip times number of strips.
  int NumberOfSegmentation =  
    (int)(1 + chargedivisionsPerStrip*fabs(driftXPos(hit->exitPoint(), driftdir, moduleThickness)-
                                           driftXPos(hit->entryPoint(), driftdir, moduleThickness) )
                                     /det.specificTopology().localPitch(hit->localPosition())         ); 
 
  // Eloss in GeV
  float eLoss = hit->energyLoss();

  // signal after pulse shape correction
  float decSignal = TimeResponse(hit, det);

  // Prepare output
  ionization_type _ionization_points;
  _ionization_points.resize(NumberOfSegmentation);

  // Fluctuate charge in track subsegments
  LocalVector direction = hit->exitPoint() - hit->entryPoint();  
  float* eLossVector = new float[NumberOfSegmentation];
  if( fluctuateCharge ) {
    fluctuateEloss(hit->particleType(), hit->pabs(), eLoss, direction.mag(), NumberOfSegmentation, eLossVector);   
    // Save the energy of each segment
    for ( int i = 0; i != NumberOfSegmentation; i++) {
      // take energy value from vector eLossVector, 
      _ionization_points[i] = EnergyDepositUnit(eLossVector[i]*decSignal/eLoss,
                                                hit->entryPoint()+float((i+0.5)/NumberOfSegmentation)*direction);
    }
  } else {
    // Save the energy of each segment
    for ( int i = 0; i != NumberOfSegmentation; i++) {
      // take energy value from eLoss average over n.segments.
      _ionization_points[i] = EnergyDepositUnit(decSignal/float(NumberOfSegmentation),
                                                hit->entryPoint()+float((i+0.5)/NumberOfSegmentation)*direction);
    }
  }
 
  delete[] eLossVector;
  return _ionization_points;
}

float SiLinearChargeDivider::driftXPos	(	const Local3DPoint &	pos,
		const LocalVector &	drift,
		double	thickness
	)		`[inline, private]`

Definition at line 56 of file SiLinearChargeDivider.h.

References PV3DBase< T, PVType, FrameType >::x(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by divide().

                                                                                              { 
    return pos.x()+(thickness/2.-pos.z())*drift.x()/drift.z();
  }

void SiLinearChargeDivider::fluctuateEloss	(	int	particleId,
		float	momentum,
		float	eloss,
		float	length,
		int	NumberOfSegmentation,
		float	elossVector[]
	)		`[private]`

Definition at line 73 of file SiLinearChargeDivider.cc.

References deltaCut, alignCSCRings::e, fluctuate, i, LogDebug, and theParticleDataTable.

Referenced by divide().

                                                                                 {

  // Get the nass if the particle, in MeV.
  // Protect from particles with Mass = 0, assuming then the pion mass
  assert(theParticleDataTable != 0);
  ParticleData const * particle = theParticleDataTable->particle( pid );
  double particleMass = particle ? particle->mass()*1000 : 139.57;
  if(!particle) {
    LogDebug("SiLinearChargeDivider") << "Cannot find particle of type "<<pid
                                      << " in the PDT we assign to this particle the mass of the Pion";
  }
  if(fabs(particleMass)<1.e-6 || pid == 22) particleMass = 139.57;

  // Generate charge fluctuations.
  float sum=0.;
  double deltaCutoff;
  double mom = particleMomentum*1000.;
  double seglen = length/NumberOfSegs*10.;
  double segeloss = (1000.*eloss)/NumberOfSegs;
  for (int i=0;i<NumberOfSegs;i++) {
    // The G4 routine needs momentum in MeV, mass in MeV, delta-cut in MeV,
    // track segment length in mm, segment eloss in MeV 
    // Returns fluctuated eloss in MeV
    // the cutoff is sometimes redefined inside, so fix it.
    deltaCutoff = deltaCut;
    sum += (elossVector[i] = fluctuate->SampleFluctuations(mom,particleMass,deltaCutoff,seglen,segeloss)/1000.);
  }
  
  if(sum>0.) {  // If fluctuations give eloss>0.
    // Rescale to the same total eloss
    float ratio = eloss/sum;
    for (int ii=0;ii<NumberOfSegs;ii++) elossVector[ii]= ratio*elossVector[ii];
  } else {  // If fluctuations gives 0 eloss
    float averageEloss = eloss/NumberOfSegs;
    for (int ii=0;ii<NumberOfSegs;ii++) elossVector[ii]= averageEloss; 
  }
  return;
}

float SiLinearChargeDivider::PeakShape	(	const PSimHit *	hit,
		const StripGeomDetUnit &	det
	)		`[private]`

Definition at line 114 of file SiLinearChargeDivider.cc.

References cosmicShift, PSimHit::energyLoss(), PSimHit::localPosition(), mag(), peakValues, GeomDet::surface(), PSimHit::tof(), Surface::toGlobal(), and x.

Referenced by TimeResponse().

                                                                                     {
  // x is difference between the tof and the tof for a photon (reference)
  // converted into a bin number
  int x = int(((det.surface().toGlobal(hit->localPosition()).mag()/30.) + cosmicShift - hit->tof())*2)+120;
  if(x < 0 || x > 920) return 0;
  return hit->energyLoss()*peakValues[x];
}

void SiLinearChargeDivider::setParticleDataTable ( const ParticleDataTable * pdt ) [inline, virtual]

Implements SiChargeDivider.

Definition at line 40 of file SiLinearChargeDivider.h.

References theParticleDataTable.

{ theParticleDataTable = pdt; }

float SiLinearChargeDivider::TimeResponse	(	const PSimHit *	hit,
		const StripGeomDetUnit &	det
	)		`[inline, private]`

Definition at line 62 of file SiLinearChargeDivider.h.

References DeconvolutionShape(), peakMode, and PeakShape().

Referenced by divide().

                                                                              {
    return (peakMode ? PeakShape(hit,det) : DeconvolutionShape(hit,det));
  }