#include <SiPixelDigitizerAlgorithm.h>

Classes
class	Amplitude
class	CalParameters
class	EnergyDepositUnit
struct	PixelEfficiencies
class	SignalPoint
Public Member Functions
void	accumulateSimHits (const std::vector< PSimHit >::const_iterator inputBegin, const std::vector< PSimHit >::const_iterator inputEnd, const PixelGeomDetUnit *pixdet, const GlobalVector &bfield)
void	digitize (const PixelGeomDetUnit pixdet, std::vector< PixelDigi > &digis, std::vector< PixelDigiSimLink > &simlinks, const TrackerTopology tTopo)
void	init (const edm::EventSetup &es)
void	initializeEvent ()
	SiPixelDigitizerAlgorithm (const edm::ParameterSet &conf, CLHEP::HepRandomEngine &)
	~SiPixelDigitizerAlgorithm ()
Private Types
typedef GloballyPositioned < double >	Frame
typedef std::vector < edm::ParameterSet >	Parameters
typedef signal_map_type::const_iterator	signal_map_const_iterator
typedef signal_map_type::iterator	signal_map_iterator
typedef std::map< int, Amplitude, std::less< int > >	signal_map_type
typedef std::map< uint32_t, signal_map_type >	signalMaps
typedef std::map< unsigned int, std::vector< float > , std::less< unsigned int > >	simlink_map
Private Member Functions
void	add_noise (const PixelGeomDetUnit *pixdet, float thePixelThreshold)
double	calcQ (float x) const
void	drift (const PSimHit &hit, const PixelGeomDetUnit *pixdet, const GlobalVector &bfield, const std::vector< EnergyDepositUnit > &ionization_points, std::vector< SignalPoint > &collection_points) const
LocalVector	DriftDirection (const PixelGeomDetUnit *pixdet, const GlobalVector &bfield, const DetId &detId) const
void	fluctuateEloss (int particleId, float momentum, float eloss, float length, int NumberOfSegments, float elossVector[]) const
void	induce_signal (const PSimHit &hit, const PixelGeomDetUnit *pixdet, const std::vector< SignalPoint > &collection_points)
std::map< int, CalParameters, std::less< int > >	initCal () const
void	make_digis (float thePixelThresholdInE, uint32_t detID, std::vector< PixelDigi > &digis, std::vector< PixelDigiSimLink > &simlinks, const TrackerTopology *tTopo) const
float	missCalibrate (uint32_t detID, int col, int row, float amp) const
void	module_killing_conf (uint32_t detID)
void	module_killing_DB (uint32_t detID)
void	pixel_inefficiency (const PixelEfficiencies &eff, const PixelGeomDetUnit pixdet, const TrackerTopology tTopo)
void	pixel_inefficiency_db (uint32_t detID)
void	primary_ionization (const PSimHit &hit, std::vector< EnergyDepositUnit > &ionization_points) const
Private Attributes
signalMaps	_signal
const bool	addChargeVCALSmearing
const bool	addNoise
const bool	addNoisyPixels
const bool	AddPixelInefficiency
const bool	addThresholdSmearing
const bool	alpha2Order
const float	BPix_p0
const float	BPix_p1
const float	BPix_p2
const float	BPix_p3
const std::map< int, CalParameters, std::less< int > >	calmap
const float	ClusterWidth
const Parameters	DeadModules
const float	Dist300
const bool	doMissCalibrate
const double	electronsPerVCAL
const double	electronsPerVCAL_Offset
const std::unique_ptr < CLHEP::RandFlat >	flatDistribution_
const std::unique_ptr < SiG4UniversalFluctuation >	fluctuate
const bool	fluctuateCharge
const float	FPix_p0
const float	FPix_p1
const float	FPix_p2
const float	FPix_p3
const std::unique_ptr < CLHEP::RandGaussQ >	gaussDistribution_
const std::unique_ptr < CLHEP::RandGaussQ >	gaussDistributionVCALNoise_
edm::ESHandle< TrackerGeometry >	geom_
const float	GeVperElectron
const bool	makeDigiSimLinks_
edm::ESHandle < SiPixelFedCablingMap >	map_
const int	NumberOfBarrelLayers
const int	NumberOfEndcapDisks
const PixelEfficiencies	pixelEfficiencies_
const double	pseudoRadDamage
const double	pseudoRadDamageRadius
const float	Sigma0
edm::ESHandle< SiPixelQuality >	SiPixelBadModule_
edm::ESHandle < SiPixelLorentzAngle >	SiPixelLorentzAngle_
const std::unique_ptr < CLHEP::RandGaussQ >	smearedThreshold_BPix_
const std::unique_ptr < CLHEP::RandGaussQ >	smearedThreshold_BPix_L1_
const std::unique_ptr < CLHEP::RandGaussQ >	smearedThreshold_FPix_
const float	tanLorentzAnglePerTesla_BPix
const float	tanLorentzAnglePerTesla_FPix
const int	theAdcFullScale
const int	theAdcFullScaleStack
const float	theElectronPerADC
const int	theFirstStackLayer
const float	theGainSmearing
const float	theNoiseInElectrons
const std::unique_ptr < GaussianTailNoiseGenerator >	theNoiser
const float	theOffsetSmearing
const float	theReadoutNoise
const std::unique_ptr < SiPixelGainCalibrationOfflineSimService >	theSiPixelGainCalibrationService_
const float	theThresholdInE_BPix
const float	theThresholdInE_BPix_L1
const float	theThresholdInE_FPix
const double	theThresholdSmearing_BPix
const double	theThresholdSmearing_BPix_L1
const double	theThresholdSmearing_FPix
const float	theTofLowerCut
const float	theTofUpperCut
const double	tMax
const bool	use_deadmodule_DB_
const bool	use_ineff_from_db_
const bool	use_LorentzAngle_DB_
const bool	use_module_killing_

Detailed Description

Definition at line 44 of file SiPixelDigitizerAlgorithm.h.

Member Typedef Documentation

typedef GloballyPositioned<double> SiPixelDigitizerAlgorithm::Frame [private]

Definition at line 246 of file SiPixelDigitizerAlgorithm.h.

typedef std::vector<edm::ParameterSet> SiPixelDigitizerAlgorithm::Parameters [private]

Definition at line 247 of file SiPixelDigitizerAlgorithm.h.

typedef signal_map_type::const_iterator SiPixelDigitizerAlgorithm::signal_map_const_iterator [private]

Definition at line 243 of file SiPixelDigitizerAlgorithm.h.

typedef signal_map_type::iterator SiPixelDigitizerAlgorithm::signal_map_iterator [private]

Definition at line 242 of file SiPixelDigitizerAlgorithm.h.

typedef std::map<int, Amplitude, std::less<int> > SiPixelDigitizerAlgorithm::signal_map_type [private]

Definition at line 241 of file SiPixelDigitizerAlgorithm.h.

typedef std::map<uint32_t, signal_map_type> SiPixelDigitizerAlgorithm::signalMaps [private]

Definition at line 245 of file SiPixelDigitizerAlgorithm.h.

typedef std::map<unsigned int, std::vector<float>,std::less<unsigned int> > SiPixelDigitizerAlgorithm::simlink_map [private]

Definition at line 244 of file SiPixelDigitizerAlgorithm.h.

Constructor & Destructor Documentation

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

Definition at line 125 of file SiPixelDigitizerAlgorithm.cc.

References AddPixelInefficiency, theAdcFullScale, theElectronPerADC, theThresholdInE_BPix, theThresholdInE_BPix_L1, theThresholdInE_FPix, and tMax.

                                                                                                           :

  _signal(),
  makeDigiSimLinks_(conf.getUntrackedParameter<bool>("makeDigiSimLinks", true)),
  use_ineff_from_db_(conf.getParameter<bool>("useDB")),
  use_module_killing_(conf.getParameter<bool>("killModules")), // boolean to kill or not modules
  use_deadmodule_DB_(conf.getParameter<bool>("DeadModules_DB")), // boolean to access dead modules from DB
  use_LorentzAngle_DB_(conf.getParameter<bool>("LorentzAngle_DB")), // boolean to access Lorentz angle from DB

  DeadModules(use_deadmodule_DB_ ? Parameters() : conf.getParameter<Parameters>("DeadModules")), // get dead module from cfg file

  // Common pixel parameters
  // These are parameters which are not likely to be changed
  GeVperElectron(3.61E-09), // 1 electron(3.61eV, 1keV(277e, mod 9/06 d.k.
  Sigma0(0.00037),           // Charge diffusion constant 7->3.7
  Dist300(0.0300),          //   normalized to 300micron Silicon
  alpha2Order(conf.getParameter<bool>("Alpha2Order")),   // switch on/off of E.B effect
  ClusterWidth(3.),     // Charge integration spread on the collection plane

  // get external parameters:
  // To account for upgrade geometries do not assume the number 
  // of layers or disks.
  NumberOfBarrelLayers(conf.exists("NumPixelBarrel")?conf.getParameter<int>("NumPixelBarrel"):3),
  NumberOfEndcapDisks(conf.exists("NumPixelEndcap")?conf.getParameter<int>("NumPixelEndcap"):2),

  // ADC calibration 1adc count(135e.
  // Corresponds to 2adc/kev, 270[e/kev]/135[e/adc](2[adc/kev]
  // Be carefull, this parameter is also used in SiPixelDet.cc to
  // calculate the noise in adc counts from noise in electrons.
  // Both defaults should be the same.
  theElectronPerADC(conf.getParameter<double>("ElectronPerAdc")),

  // ADC saturation value, 255(8bit adc.
  //theAdcFullScale(conf.getUntrackedParameter<int>("AdcFullScale",255)),
  theAdcFullScale(conf.getParameter<int>("AdcFullScale")),
  theAdcFullScaleStack(conf.exists("AdcFullScaleStack")?conf.getParameter<int>("AdcFullScaleStack"):255),
  theFirstStackLayer(conf.exists("FirstStackLayer")?conf.getParameter<int>("FirstStackLayer"):5),

  // Noise in electrons:
  // Pixel cell noise, relevant for generating noisy pixels
  theNoiseInElectrons(conf.getParameter<double>("NoiseInElectrons")),

  // Fill readout noise, including all readout chain, relevant for smearing
  //theReadoutNoise(conf.getUntrackedParameter<double>("ReadoutNoiseInElec",500.)),
  theReadoutNoise(conf.getParameter<double>("ReadoutNoiseInElec")),

  // Pixel threshold in units of noise:
  // thePixelThreshold(conf.getParameter<double>("ThresholdInNoiseUnits")),
  // Pixel threshold in electron units.
  theThresholdInE_FPix(conf.getParameter<double>("ThresholdInElectrons_FPix")),
  theThresholdInE_BPix(conf.getParameter<double>("ThresholdInElectrons_BPix")),
  theThresholdInE_BPix_L1(conf.exists("ThresholdInElectrons_BPix_L1")?conf.getParameter<double>("ThresholdInElectrons_BPix_L1"):theThresholdInE_BPix),

  // Add threshold gaussian smearing:
  theThresholdSmearing_FPix(conf.getParameter<double>("ThresholdSmearing_FPix")),
  theThresholdSmearing_BPix(conf.getParameter<double>("ThresholdSmearing_BPix")),
  theThresholdSmearing_BPix_L1(conf.exists("ThresholdSmearing_BPix_L1")?conf.getParameter<double>("ThresholdSmearing_BPix_L1"):theThresholdSmearing_BPix),

  // electrons to VCAL conversion needed in misscalibrate()
  electronsPerVCAL(conf.getParameter<double>("ElectronsPerVcal")),
  electronsPerVCAL_Offset(conf.getParameter<double>("ElectronsPerVcal_Offset")),

  //theTofCut 12.5, cut in particle TOD +/- 12.5ns
  //theTofCut(conf.getUntrackedParameter<double>("TofCut",12.5)),
  theTofLowerCut(conf.getParameter<double>("TofLowerCut")),
  theTofUpperCut(conf.getParameter<double>("TofUpperCut")),

  // Get the Lorentz angle from the cfg file:
  tanLorentzAnglePerTesla_FPix(use_LorentzAngle_DB_ ? 0.0 : conf.getParameter<double>("TanLorentzAnglePerTesla_FPix")),
  tanLorentzAnglePerTesla_BPix(use_LorentzAngle_DB_ ? 0.0 : conf.getParameter<double>("TanLorentzAnglePerTesla_BPix")),

  // signal response new parameterization: split Fpix and BPix
  FPix_p0(conf.getParameter<double>("FPix_SignalResponse_p0")),
  FPix_p1(conf.getParameter<double>("FPix_SignalResponse_p1")),
  FPix_p2(conf.getParameter<double>("FPix_SignalResponse_p2")),
  FPix_p3(conf.getParameter<double>("FPix_SignalResponse_p3")),

  BPix_p0(conf.getParameter<double>("BPix_SignalResponse_p0")),
  BPix_p1(conf.getParameter<double>("BPix_SignalResponse_p1")),
  BPix_p2(conf.getParameter<double>("BPix_SignalResponse_p2")),
  BPix_p3(conf.getParameter<double>("BPix_SignalResponse_p3")),

  // Add noise
  addNoise(conf.getParameter<bool>("AddNoise")),

  // Smear the pixel charge with a gaussian which RMS is a function of the
  // pixel charge (Danek's study)
  addChargeVCALSmearing(conf.getParameter<bool>("ChargeVCALSmearing")),

  // Add noisy pixels
  addNoisyPixels(conf.getParameter<bool>("AddNoisyPixels")),

  // Fluctuate charge in track subsegments
  fluctuateCharge(conf.getUntrackedParameter<bool>("FluctuateCharge",true)),

  // Control the pixel inefficiency
  AddPixelInefficiency(conf.getParameter<bool>("AddPixelInefficiencyFromPython")),

  // Add threshold gaussian smearing:
  addThresholdSmearing(conf.getParameter<bool>("AddThresholdSmearing")),

  // Get the constants for the miss-calibration studies
  doMissCalibrate(conf.getParameter<bool>("MissCalibrate")), // Enable miss-calibration
  theGainSmearing(conf.getParameter<double>("GainSmearing")), // sigma of the gain smearing
  theOffsetSmearing(conf.getParameter<double>("OffsetSmearing")), //sigma of the offset smearing

  // Add some pseudo-red damage
  pseudoRadDamage(conf.exists("PseudoRadDamage")?conf.getParameter<double>("PseudoRadDamage"):double(0.0)),
  pseudoRadDamageRadius(conf.exists("PseudoRadDamageRadius")?conf.getParameter<double>("PseudoRadDamageRadius"):double(0.0)),

  // delta cutoff in MeV, has to be same as in OSCAR(0.030/cmsim=1.0 MeV
  //tMax(0.030), // In MeV.
  //tMax(conf.getUntrackedParameter<double>("DeltaProductionCut",0.030)),
  tMax(conf.getParameter<double>("DeltaProductionCut")),

  fluctuate(fluctuateCharge ? new SiG4UniversalFluctuation(eng) : 0),
  theNoiser(addNoise ? new GaussianTailNoiseGenerator(eng) : 0),
  calmap(doMissCalibrate ? initCal() : std::map<int,CalParameters,std::less<int> >()),
  theSiPixelGainCalibrationService_(use_ineff_from_db_ ? new SiPixelGainCalibrationOfflineSimService(conf) : 0),
  pixelEfficiencies_(conf, AddPixelInefficiency,NumberOfBarrelLayers,NumberOfEndcapDisks),
  flatDistribution_((addNoise || AddPixelInefficiency || fluctuateCharge || addThresholdSmearing) ? new CLHEP::RandFlat(eng, 0., 1.) : 0),
  gaussDistribution_((addNoise || AddPixelInefficiency || fluctuateCharge || addThresholdSmearing) ? new CLHEP::RandGaussQ(eng, 0., theReadoutNoise) : 0),
  gaussDistributionVCALNoise_((addNoise || AddPixelInefficiency || fluctuateCharge || addThresholdSmearing) ? new CLHEP::RandGaussQ(eng, 0., 1.) : 0),
  // Threshold smearing with gaussian distribution:
  smearedThreshold_FPix_(addThresholdSmearing ? new CLHEP::RandGaussQ(eng, theThresholdInE_FPix , theThresholdSmearing_FPix) : 0),
  smearedThreshold_BPix_(addThresholdSmearing ? new CLHEP::RandGaussQ(eng, theThresholdInE_BPix , theThresholdSmearing_BPix) : 0),
  smearedThreshold_BPix_L1_(addThresholdSmearing ? new CLHEP::RandGaussQ(eng, theThresholdInE_BPix_L1 , theThresholdSmearing_BPix_L1) : 0)
{
  LogInfo ("PixelDigitizer ") <<"SiPixelDigitizerAlgorithm constructed"
                              <<"Configuration parameters:"
                              << "Threshold/Gain = "
                              << "threshold in electron FPix = "
                              << theThresholdInE_FPix
                              << "threshold in electron BPix = "
                              << theThresholdInE_BPix
                              << "threshold in electron BPix Layer1 = "
                              << theThresholdInE_BPix_L1
                              <<" " << theElectronPerADC << " " << theAdcFullScale
                              << " The delta cut-off is set to " << tMax
                              << " pix-inefficiency "<<AddPixelInefficiency;

}

SiPixelDigitizerAlgorithm::~SiPixelDigitizerAlgorithm ( )

Definition at line 361 of file SiPixelDigitizerAlgorithm.cc.

References LogDebug.

                                                      {
  LogDebug ("PixelDigitizer")<<"SiPixelDigitizerAlgorithm deleted";
}

Member Function Documentation

void SiPixelDigitizerAlgorithm::accumulateSimHits	(	const std::vector< PSimHit >::const_iterator	inputBegin,
		const std::vector< PSimHit >::const_iterator	inputEnd,
		const PixelGeomDetUnit *	pixdet,
		const GlobalVector &	bfield
	)

Definition at line 434 of file SiPixelDigitizerAlgorithm.cc.

References drift(), GeomDet::geographicalId(), induce_signal(), LogDebug, mag(), primary_ionization(), DetId::rawId(), GeomDet::surface(), theTofLowerCut, theTofUpperCut, and Surface::toGlobal().

                                                                              {
    // produce SignalPoint's for all SimHit's in detector
    // Loop over hits

    uint32_t detId = pixdet->geographicalId().rawId();
    for (std::vector<PSimHit>::const_iterator ssbegin = inputBegin; ssbegin != inputEnd; ++ssbegin) {
      // skip hits not in this detector.
      if((*ssbegin).detUnitId() != detId) {
        continue;
      }

#ifdef TP_DEBUG
      LogDebug ("Pixel Digitizer")
        << (*ssbegin).particleType() << " " << (*ssbegin).pabs() << " "
        << (*ssbegin).energyLoss() << " " << (*ssbegin).tof() << " "
        << (*ssbegin).trackId() << " " << (*ssbegin).processType() << " "
        << (*ssbegin).detUnitId()
        << (*ssbegin).entryPoint() << " " << (*ssbegin).exitPoint() ;
#endif

      
      std::vector<EnergyDepositUnit> ionization_points;
      std::vector<SignalPoint> collection_points;

      // fill collection_points for this SimHit, indpendent of topology
      // Check the TOF cut
      if (  ((*ssbegin).tof() - pixdet->surface().toGlobal((*ssbegin).localPosition()).mag()/30.)>= theTofLowerCut &&
            ((*ssbegin).tof()- pixdet->surface().toGlobal((*ssbegin).localPosition()).mag()/30.) <= theTofUpperCut ) {
        primary_ionization(*ssbegin, ionization_points); // fills _ionization_points
        drift(*ssbegin, pixdet, bfield, ionization_points, collection_points);  // transforms _ionization_points to collection_points
        // compute induced signal on readout elements and add to _signal
        induce_signal(*ssbegin, pixdet, collection_points); // *ihit needed only for SimHit<-->Digi link
      } //  end if
    } // end for

}

void SiPixelDigitizerAlgorithm::add_noise	(	const PixelGeomDetUnit *	pixdet,
		float	thePixelThreshold
	)		`[private]`

Definition at line 1149 of file SiPixelDigitizerAlgorithm.cc.

References _signal, addChargeVCALSmearing, addNoisyPixels, gaussDistribution_, gaussDistributionVCALNoise_, GeomDet::geographicalId(), i, LogDebug, PixelTopology::ncolumns(), PixelTopology::nrows(), PixelDigi::pixelToChannel(), DetId::rawId(), PixelGeomDetUnit::specificTopology(), theNoiseInElectrons, theNoiser, theThresholdInE_BPix, and theThresholdInE_FPix.

Referenced by digitize().

                                                                   {

#ifdef TP_DEBUG
  LogDebug ("Pixel Digitizer") << " enter add_noise " << theNoiseInElectrons;
#endif

  uint32_t detID= pixdet->geographicalId().rawId();
  signal_map_type& theSignal = _signal[detID];


  // First add noise to hit pixels
  float theSmearedChargeRMS = 0.0;

  for ( signal_map_iterator i = theSignal.begin(); i != theSignal.end(); i++) {

         if(addChargeVCALSmearing)
      {
        if((*i).second < 3000)
          {
            theSmearedChargeRMS = 543.6 - (*i).second * 0.093;
          } else if((*i).second < 6000){
            theSmearedChargeRMS = 307.6 - (*i).second * 0.01;
          } else{
            theSmearedChargeRMS = -432.4 +(*i).second * 0.123;
        }

        // Noise from Vcal smearing:
        float noise_ChargeVCALSmearing = theSmearedChargeRMS * gaussDistributionVCALNoise_->fire() ;
        // Noise from full readout:
        float noise  = gaussDistribution_->fire() ;

                if(((*i).second + Amplitude(noise+noise_ChargeVCALSmearing, -1.)) < 0. ) {
                  (*i).second.set(0);}
                else{
        (*i).second +=Amplitude(noise+noise_ChargeVCALSmearing, -1.);
                }

      } // End if addChargeVCalSmearing
         else
     {
        // Noise: ONLY full READOUT Noise.
        // Use here the FULL readout noise, including TBM,ALT,AOH,OPT-REC.
        float noise  = gaussDistribution_->fire() ;

                if(((*i).second + Amplitude(noise, -1.)) < 0. ) {
                  (*i).second.set(0);}
                else{
        (*i).second +=Amplitude(noise, -1.);
                }
     } // end if only Noise from full readout

  }

  if(!addNoisyPixels)  // Option to skip noise in non-hit pixels
    return;

  const PixelTopology* topol=&pixdet->specificTopology();
  int numColumns = topol->ncolumns();  // det module number of cols&rows
  int numRows = topol->nrows();

  // Add noise on non-hit pixels
  // Use here the pixel noise
  int numberOfPixels = (numRows * numColumns);
  std::map<int,float, std::less<int> > otherPixels;
  std::map<int,float, std::less<int> >::iterator mapI;

  theNoiser->generate(numberOfPixels,
                      thePixelThreshold, //thr. in un. of nois
                      theNoiseInElectrons, // noise in elec.
                      otherPixels );

#ifdef TP_DEBUG
  LogDebug ("Pixel Digitizer")
    <<  " Add noisy pixels " << numRows << " "
    << numColumns << " " << theNoiseInElectrons << " "
    << theThresholdInE_FPix << theThresholdInE_BPix <<" "<< numberOfPixels<<" "
    << otherPixels.size() ;
#endif

  // Add noisy pixels
  for (mapI = otherPixels.begin(); mapI!= otherPixels.end(); mapI++) {
    int iy = ((*mapI).first) / numRows;
    int ix = ((*mapI).first) - (iy*numRows);

    // Keep for a while for testing.
    if( iy < 0 || iy > (numColumns-1) )
      LogWarning ("Pixel Geometry") << " error in iy " << iy ;
    if( ix < 0 || ix > (numRows-1) )
      LogWarning ("Pixel Geometry")  << " error in ix " << ix ;

    int chan = PixelDigi::pixelToChannel(ix, iy);

#ifdef TP_DEBUG
    LogDebug ("Pixel Digitizer")
      <<" Storing noise = " << (*mapI).first << " " << (*mapI).second
      << " " << ix << " " << iy << " " << chan ;
#endif

    if(theSignal[chan] == 0){
      //      float noise = float( (*mapI).second );
      int noise=int( (*mapI).second );
      theSignal[chan] = Amplitude (noise, -1.);
    }
  }
}

double SiPixelDigitizerAlgorithm::calcQ ( float x ) const [inline, private]

Definition at line 398 of file SiPixelDigitizerAlgorithm.h.

References f, min, and mathSSE::sqrt().

Referenced by induce_signal().

                                {
      // need erf(x/sqrt2)
      //float x2=0.5*x*x;
      //float a=0.147;
      //double erf=sqrt(1.0f-exp( -1.0f*x2*( (4/M_PI)+a*x2)/(1.0+a*x2)));
      //if (x<0.) erf*=-1.0;
      //return 0.5*(1.0-erf);

      auto xx=std::min(0.5f*x*x,12.5f);
      return 0.5*(1.0-std::copysign(std::sqrt(1.f- unsafe_expf<4>(-xx*(1.f+0.2733f/(1.f+0.147f*xx)) )),x));
    }

void SiPixelDigitizerAlgorithm::digitize	(	const PixelGeomDetUnit *	pixdet,
		std::vector< PixelDigi > &	digis,
		std::vector< PixelDigiSimLink > &	simlinks,
		const TrackerTopology *	tTopo
	)

Definition at line 475 of file SiPixelDigitizerAlgorithm.cc.

References _signal, add_noise(), addNoise, AddPixelInefficiency, addThresholdSmearing, Surface::bounds(), GeomDet::geographicalId(), LogDebug, make_digis(), module_killing_conf(), module_killing_DB(), PixelTopology::ncolumns(), PixelTopology::nrows(), pixel_inefficiency(), pixel_inefficiency_db(), PixelSubdetector::PixelBarrel, pixelEfficiencies_, TrackerTopology::pxbLayer(), DetId::rawId(), smearedThreshold_BPix_, smearedThreshold_BPix_L1_, smearedThreshold_FPix_, GeomDet::specificSurface(), PixelGeomDetUnit::specificTopology(), DetId::subdetId(), theNoiseInElectrons, theThresholdInE_BPix, theThresholdInE_BPix_L1, theThresholdInE_FPix, Bounds::thickness(), use_deadmodule_DB_, use_ineff_from_db_, and use_module_killing_.

                                                                                                              {

   // Pixel Efficiency moved from the constructor to this method because
   // the information of the det are not available in the constructor
   // Effciency parameters. 0 - no inefficiency, 1-low lumi, 10-high lumi

   uint32_t detID = pixdet->geographicalId().rawId();
   const signal_map_type& theSignal = _signal[detID];

   const PixelTopology* topol=&pixdet->specificTopology();
   int numColumns = topol->ncolumns();  // det module number of cols&rows
   int numRows = topol->nrows();

    // Noise already defined in electrons
    //thePixelThresholdInE = thePixelThreshold * theNoiseInElectrons ;
    // Find the threshold in noise units, needed for the noiser.

  unsigned int Sub_detid=DetId(detID).subdetId();

  float thePixelThresholdInE = 0.;

  if(theNoiseInElectrons>0.){
    if(Sub_detid == PixelSubdetector::PixelBarrel){ // Barrel modules
      int lay = tTopo->pxbLayer(detID);
      if(addThresholdSmearing) {
        if(lay==1) {
          thePixelThresholdInE = smearedThreshold_BPix_L1_->fire(); // gaussian smearing
        } else {
          thePixelThresholdInE = smearedThreshold_BPix_->fire(); // gaussian smearing
        }
      } else {
        if(lay==1) {
          thePixelThresholdInE = theThresholdInE_BPix_L1;
        } else {
          thePixelThresholdInE = theThresholdInE_BPix; // no smearing
        }
      }
    } else { // Forward disks modules
      if(addThresholdSmearing) {
        thePixelThresholdInE = smearedThreshold_FPix_->fire(); // gaussian smearing
      } else {
        thePixelThresholdInE = theThresholdInE_FPix; // no smearing
      }
    }
  }


#ifdef TP_DEBUG
    // full detector thickness
   float moduleThickness = pixdet->specificSurface().bounds().thickness();
    LogDebug ("PixelDigitizer")
      << " PixelDigitizer "
      << numColumns << " " << numRows << " " << moduleThickness;
#endif

    if(addNoise) add_noise(pixdet, thePixelThresholdInE/theNoiseInElectrons);  // generate noise

    // Do only if needed

    if((AddPixelInefficiency) && (theSignal.size()>0))
      pixel_inefficiency(pixelEfficiencies_, pixdet, tTopo); // Kill some pixels

    if(use_ineff_from_db_ && (theSignal.size()>0))
      pixel_inefficiency_db(detID);

    if(use_module_killing_) {
      if (use_deadmodule_DB_) {  // remove dead modules using DB
        module_killing_DB(detID);
      } else { // remove dead modules using the list in cfg file
        module_killing_conf(detID);
      }
    }

    make_digis(thePixelThresholdInE, detID, digis, simlinks, tTopo);

#ifdef TP_DEBUG
  LogDebug ("PixelDigitizer") << "[SiPixelDigitizerAlgorithm] converted " << digis.size() << " PixelDigis in DetUnit" << detID;
#endif
}

void SiPixelDigitizerAlgorithm::drift	(	const PSimHit &	hit,
		const PixelGeomDetUnit *	pixdet,
		const GlobalVector &	bfield,
		const std::vector< EnergyDepositUnit > &	ionization_points,
		std::vector< SignalPoint > &	collection_points
	)		const `[private]`

Definition at line 687 of file SiPixelDigitizerAlgorithm.cc.

References alpha2Order, Surface::bounds(), PSimHit::detUnitId(), Dist300, DriftDirection(), create_public_lumi_plots::exp, i, LogDebug, GloballyPositioned< T >::position(), pseudoRadDamage, pseudoRadDamageRadius, Sigma0, GeomDet::specificSurface(), mathSSE::sqrt(), GeomDet::surface(), Bounds::thickness(), PSimHit::tof(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by accumulateSimHits().

                                                                                       {

#ifdef TP_DEBUG
  LogDebug ("Pixel Digitizer") << " enter drift " ;
#endif

  collection_points.resize(ionization_points.size()); // set size

  LocalVector driftDir=DriftDirection(pixdet, bfield, hit.detUnitId());  // get the charge drift direction
  if(driftDir.z() ==0.) {
    LogWarning("Magnetic field") << " pxlx: drift in z is zero ";
    return;
  }
  
  // tangent of Lorentz angle
  //float TanLorenzAngleX = driftDir.x()/driftDir.z();
  //float TanLorenzAngleY = 0.; // force to 0, driftDir.y()/driftDir.z();
  
  float TanLorenzAngleX, TanLorenzAngleY,dir_z, CosLorenzAngleX,
    CosLorenzAngleY;
  if( alpha2Order) {

    TanLorenzAngleX = driftDir.x(); // tangen of Lorentz angle
    TanLorenzAngleY = driftDir.y();
    dir_z = driftDir.z(); // The z drift direction
    CosLorenzAngleX = 1./sqrt(1.+TanLorenzAngleX*TanLorenzAngleX); //cosine
    CosLorenzAngleY = 1./sqrt(1.+TanLorenzAngleY*TanLorenzAngleY); //cosine;

  } else{

    TanLorenzAngleX = driftDir.x();
    TanLorenzAngleY = 0.; // force to 0, driftDir.y()/driftDir.z();
    dir_z = driftDir.z(); // The z drift direction
    CosLorenzAngleX = 1./sqrt(1.+TanLorenzAngleX*TanLorenzAngleX); //cosine to estimate the path length
    CosLorenzAngleY = 1.;
  }

  float moduleThickness = pixdet->specificSurface().bounds().thickness();
#ifdef TP_DEBUG
  LogDebug ("Pixel Digitizer")
    << " Lorentz Tan " << TanLorenzAngleX << " " << TanLorenzAngleY <<" "
    << CosLorenzAngleX << " " << CosLorenzAngleY << " "
    << moduleThickness*TanLorenzAngleX << " " << driftDir;
#endif

  float Sigma_x = 1.;  // Charge spread
  float Sigma_y = 1.;
  float DriftDistance; // Distance between charge generation and collection
  float DriftLength;   // Actual Drift Lentgh
  float Sigma;


  for (unsigned int i = 0; i != ionization_points.size(); i++) {

    float SegX, SegY, SegZ; // position
    SegX = ionization_points[i].x();
    SegY = ionization_points[i].y();
    SegZ = ionization_points[i].z();

    // Distance from the collection plane
    //DriftDistance = (moduleThickness/2. + SegZ); // Drift to -z
    // Include explixitely the E drift direction (for CMS dir_z=-1)
    DriftDistance = moduleThickness/2. - (dir_z * SegZ); // Drift to -z

    //if( DriftDistance <= 0.)
    //cout<<" <=0 "<<DriftDistance<<" "<<i<<" "<<SegZ<<" "<<dir_z<<" "
    //  <<SegX<<" "<<SegY<<" "<<(moduleThickness/2)<<" "
    //  <<ionization_points[i].energy()<<" "
    //  <<hit.particleType()<<" "<<hit.pabs()<<" "<<hit.energyLoss()<<" "
    //  <<hit.entryPoint()<<" "<<hit.exitPoint()
    //  <<std::endl;

    if( DriftDistance < 0.) {
      DriftDistance = 0.;
    } else if( DriftDistance > moduleThickness )
      DriftDistance = moduleThickness;

    // Assume full depletion now, partial depletion will come later.
    float XDriftDueToMagField = DriftDistance * TanLorenzAngleX;
    float YDriftDueToMagField = DriftDistance * TanLorenzAngleY;

    // Shift cloud center
    float CloudCenterX = SegX + XDriftDueToMagField;
    float CloudCenterY = SegY + YDriftDueToMagField;

    // Calculate how long is the charge drift path
    DriftLength = sqrt( DriftDistance*DriftDistance +
                        XDriftDueToMagField*XDriftDueToMagField +
                        YDriftDueToMagField*YDriftDueToMagField );

    // What is the charge diffusion after this path
    Sigma = sqrt(DriftLength/Dist300) * Sigma0;

    // Project the diffusion sigma on the collection plane
    Sigma_x = Sigma / CosLorenzAngleX ;
    Sigma_y = Sigma / CosLorenzAngleY ;

    // Insert a charge loss due to Rad Damage here
    float energyOnCollector = ionization_points[i].energy(); // The energy that reaches the collector
    // pseudoRadDamage
    if (pseudoRadDamage>=0.001){
        float moduleRadius = pixdet->surface().position().perp();
        if (moduleRadius<=pseudoRadDamageRadius){
          float kValue = pseudoRadDamage /(moduleRadius*moduleRadius);
          //std::cout <<"\nmoduleRadius= "<<moduleRadius<<" WITH K="<<kValue <<" wich scales the energy by "<<
          //            exp( -1*kValue*DriftDistance/moduleThickness )<<"\n" ;
          energyOnCollector = energyOnCollector * exp( -1*kValue*DriftDistance/moduleThickness );
          }
        }
#ifdef TP_DEBUG
  LogDebug ("Pixel Digitizer")
                <<" Dift DistanceZ= "<<DriftDistance<<" module thickness= "<<moduleThickness
                <<" Start Energy= "<<ionization_points[i].energy()<<" Energy after loss= "<<energyOnCollector;
#endif
    SignalPoint sp( CloudCenterX, CloudCenterY,
     Sigma_x, Sigma_y, hit.tof(), energyOnCollector );

    // Load the Charge distribution parameters
    collection_points[i] = (sp);

  } // loop over ionization points, i.

} // end drift

LocalVector SiPixelDigitizerAlgorithm::DriftDirection	(	const PixelGeomDetUnit *	pixdet,
		const GlobalVector &	bfield,
		const DetId &	detId
	)		const `[private]`

Definition at line 1473 of file SiPixelDigitizerAlgorithm.cc.

References alpha2Order, GeomDet::geographicalId(), LogDebug, PixelSubdetector::PixelBarrel, GloballyPositioned< T >::position(), DetId::rawId(), GloballyPositioned< T >::rotation(), pileupReCalc_HLTpaths::scale, SiPixelLorentzAngle_, DetId::subdetId(), GeomDet::surface(), tanLorentzAnglePerTesla_BPix, tanLorentzAnglePerTesla_FPix, and use_LorentzAngle_DB_.

Referenced by drift().

                                                                                {
  Frame detFrame(pixdet->surface().position(),pixdet->surface().rotation());
  LocalVector Bfield=detFrame.toLocal(bfield);
  
  float alpha2_FPix;
  float alpha2_BPix;
  float alpha2;
  
  //float dir_x = -tanLorentzAnglePerTesla * Bfield.y();
  //float dir_y = +tanLorentzAnglePerTesla * Bfield.x();
  //float dir_z = -1.; // E field always in z direction, so electrons go to -z
  // The dir_z has to be +/- 1. !
  // LocalVector theDriftDirection = LocalVector(dir_x,dir_y,dir_z);

  float dir_x = 0.0;
  float dir_y = 0.0;
  float dir_z = 0.0;
  float scale = 0.0;

  uint32_t detID= pixdet->geographicalId().rawId();

  unsigned int Sub_detid=DetId(detID).subdetId();

  // Read Lorentz angle from cfg file:**************************************************************

  if(!use_LorentzAngle_DB_){
    
    if( alpha2Order) {
      alpha2_FPix = tanLorentzAnglePerTesla_FPix*tanLorentzAnglePerTesla_FPix;
      alpha2_BPix = tanLorentzAnglePerTesla_BPix*tanLorentzAnglePerTesla_BPix;
    }else {
      alpha2_FPix = 0.0;
      alpha2_BPix = 0.0;
    }
    
    if(Sub_detid == PixelSubdetector::PixelBarrel){// barrel layers
      dir_x = -( tanLorentzAnglePerTesla_BPix * Bfield.y() + alpha2_BPix* Bfield.z()* Bfield.x() );
      dir_y = +( tanLorentzAnglePerTesla_BPix * Bfield.x() - alpha2_BPix* Bfield.z()* Bfield.y() );
      dir_z = -(1 + alpha2_BPix* Bfield.z()*Bfield.z() );
      scale = (1 + alpha2_BPix* Bfield.z()*Bfield.z() );

    } else {// forward disks
      dir_x = -( tanLorentzAnglePerTesla_FPix * Bfield.y() + alpha2_FPix* Bfield.z()* Bfield.x() );
      dir_y = +( tanLorentzAnglePerTesla_FPix * Bfield.x() - alpha2_FPix* Bfield.z()* Bfield.y() );
      dir_z = -(1 + alpha2_FPix* Bfield.z()*Bfield.z() );
      scale = (1 + alpha2_FPix* Bfield.z()*Bfield.z() );
    }
  } // end: Read LA from cfg file.

  //Read Lorentz angle from DB:********************************************************************
  if(use_LorentzAngle_DB_){
    float lorentzAngle = SiPixelLorentzAngle_->getLorentzAngle(detId);
    alpha2 = lorentzAngle * lorentzAngle;
    //std::cout << "detID is: "<< it->first <<"The LA per tesla is: "<< it->second << std::std::endl;
    dir_x = -( lorentzAngle * Bfield.y() + alpha2 * Bfield.z()* Bfield.x() );
    dir_y = +( lorentzAngle * Bfield.x() - alpha2 * Bfield.z()* Bfield.y() );
    dir_z = -(1 + alpha2 * Bfield.z()*Bfield.z() );
    scale = (1 + alpha2 * Bfield.z()*Bfield.z() );
  }// end: Read LA from DataBase.
  
  LocalVector theDriftDirection = LocalVector(dir_x/scale, dir_y/scale, dir_z/scale );
  
#ifdef TP_DEBUG
  LogDebug ("Pixel Digitizer") << " The drift direction in local coordinate is "
                               << theDriftDirection ;
#endif
  
  return theDriftDirection;
}

void SiPixelDigitizerAlgorithm::fluctuateEloss	(	int	particleId,
		float	momentum,
		float	eloss,
		float	length,
		int	NumberOfSegments,
		float	elossVector[]
	)		const `[private]`

Definition at line 632 of file SiPixelDigitizerAlgorithm.cc.

References abs, fluctuate, i, cuy::ii, and tMax.

Referenced by primary_ionization().

                                                                                  {

  // Get dedx for this track
  //float dedx;
  //if( length > 0.) dedx = eloss/length;
  //else dedx = eloss;

  double particleMass = 139.6; // Mass in MeV, Assume pion
  pid = std::abs(pid);
  if(pid!=211) {       // Mass in MeV
    if(pid==11)        particleMass = 0.511;
    else if(pid==13)   particleMass = 105.7;
    else if(pid==321)  particleMass = 493.7;
    else if(pid==2212) particleMass = 938.3;
  }
  // What is the track segment length.
  float segmentLength = length/NumberOfSegs;

  // Generate charge fluctuations.
  float de=0.;
  float sum=0.;
  double segmentEloss = (1000.*eloss)/NumberOfSegs; //eloss in MeV
  for (int i=0;i<NumberOfSegs;i++) {
    //       material,*,   momentum,energy,*, *,  mass
    //myglandz_(14.,segmentLength,2.,2.,dedx,de,0.14);
    // 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
    double deltaCutoff = tMax; // the cutoff is sometimes redefined inside, so fix it.
    de = fluctuate->SampleFluctuations(double(particleMomentum*1000.),
                                      particleMass, deltaCutoff,
                                      double(segmentLength*10.),
                                      segmentEloss )/1000.; //convert to GeV

    elossVector[i]=de;
    sum +=de;
  }

  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;
}

void SiPixelDigitizerAlgorithm::induce_signal	(	const PSimHit &	hit,
		const PixelGeomDetUnit *	pixdet,
		const std::vector< SignalPoint > &	collection_points
	)		`[private]`

Definition at line 817 of file SiPixelDigitizerAlgorithm.cc.

References _signal, calcQ(), Topology::channel(), PixelDigi::channelToPixel(), ClusterWidth, GeomDet::geographicalId(), i, Topology::localPosition(), LogDebug, makeDigiSimLinks_, Topology::measurementPosition(), PixelTopology::ncolumns(), PixelTopology::nrows(), PixelTopology::pitch(), PixelDigi::pixelToChannel(), DetId::rawId(), PixelGeomDetUnit::specificTopology(), PV3DBase< T, PVType, FrameType >::x(), PV2DBase< T, PVType, FrameType >::x(), x, PV3DBase< T, PVType, FrameType >::y(), PV2DBase< T, PVType, FrameType >::y(), and detailsBasic3DVector::y.

Referenced by accumulateSimHits().

                                                                                               {

  // X  - Rows, Left-Right, 160, (1.6cm)   for barrel
  // Y  - Columns, Down-Up, 416, (6.4cm)

   const PixelTopology* topol=&pixdet->specificTopology();
   uint32_t detID= pixdet->geographicalId().rawId();
   signal_map_type& theSignal = _signal[detID];

#ifdef TP_DEBUG
    LogDebug ("Pixel Digitizer")
      << " enter induce_signal, "
      << topol->pitch().first << " " << topol->pitch().second; //OK
#endif

   // local map to store pixels hit by 1 Hit.
   typedef std::map< int, float, std::less<int> > hit_map_type;
   hit_map_type hit_signal;

   // map to store pixel integrals in the x and in the y directions
   std::map<int, float, std::less<int> > x,y;

   // Assign signals to readout channels and store sorted by channel number

   // Iterate over collection points on the collection plane
   for ( std::vector<SignalPoint>::const_iterator i=collection_points.begin();
         i != collection_points.end(); ++i) {

     float CloudCenterX = i->position().x(); // Charge position in x
     float CloudCenterY = i->position().y(); //                 in y
     float SigmaX = i->sigma_x();            // Charge spread in x
     float SigmaY = i->sigma_y();            //               in y
     float Charge = i->amplitude();          // Charge amplitude


     //if(SigmaX==0 || SigmaY==0) {
     //cout<<SigmaX<<" "<<SigmaY
     //   << " cloud " << i->position().x() << " " << i->position().y() << " "
     //   << i->sigma_x() << " " << i->sigma_y() << " " << i->amplitude()<<std::endl;
     //}

#ifdef TP_DEBUG
       LogDebug ("Pixel Digitizer")
         << " cloud " << i->position().x() << " " << i->position().y() << " "
         << i->sigma_x() << " " << i->sigma_y() << " " << i->amplitude();
#endif

     // Find the maximum cloud spread in 2D plane , assume 3*sigma
     float CloudRight = CloudCenterX + ClusterWidth*SigmaX;
     float CloudLeft  = CloudCenterX - ClusterWidth*SigmaX;
     float CloudUp    = CloudCenterY + ClusterWidth*SigmaY;
     float CloudDown  = CloudCenterY - ClusterWidth*SigmaY;

     // Define 2D cloud limit points
     LocalPoint PointRightUp  = LocalPoint(CloudRight,CloudUp);
     LocalPoint PointLeftDown = LocalPoint(CloudLeft,CloudDown);

     // This points can be located outside the sensor area.
     // The conversion to measurement point does not check for that
     // so the returned pixel index might be wrong (outside range).
     // We rely on the limits check below to fix this.
     // But remember whatever we do here THE CHARGE OUTSIDE THE ACTIVE
     // PIXEL ARE IS LOST, it should not be collected.

     // Convert the 2D points to pixel indices
     MeasurementPoint mp = topol->measurementPosition(PointRightUp ); //OK

     int IPixRightUpX = int( floor( mp.x()));
     int IPixRightUpY = int( floor( mp.y()));

#ifdef TP_DEBUG
     LogDebug ("Pixel Digitizer") << " right-up " << PointRightUp << " "
                                  << mp.x() << " " << mp.y() << " "
                                  << IPixRightUpX << " " << IPixRightUpY ;
#endif

     mp = topol->measurementPosition(PointLeftDown ); //OK

     int IPixLeftDownX = int( floor( mp.x()));
     int IPixLeftDownY = int( floor( mp.y()));

#ifdef TP_DEBUG
     LogDebug ("Pixel Digitizer") << " left-down " << PointLeftDown << " "
                                  << mp.x() << " " << mp.y() << " "
                                  << IPixLeftDownX << " " << IPixLeftDownY ;
#endif

     // Check detector limits to correct for pixels outside range.
     int numColumns = topol->ncolumns();  // det module number of cols&rows
     int numRows = topol->nrows();

     IPixRightUpX = numRows>IPixRightUpX ? IPixRightUpX : numRows-1 ;
     IPixRightUpY = numColumns>IPixRightUpY ? IPixRightUpY : numColumns-1 ;
     IPixLeftDownX = 0<IPixLeftDownX ? IPixLeftDownX : 0 ;
     IPixLeftDownY = 0<IPixLeftDownY ? IPixLeftDownY : 0 ;

     x.clear(); // clear temporary integration array
     y.clear();

     // First integrate cahrge strips in x
     int ix; // TT for compatibility
     for (ix=IPixLeftDownX; ix<=IPixRightUpX; ix++) {  // loop over x index
       float xUB, xLB, UpperBound, LowerBound;

       // Why is set to 0 if ix=0, does it meen that we accept charge
       // outside the sensor? CHeck How it was done in ORCA?
       //if(ix == 0) LowerBound = 0.;
       if(ix == 0 || SigmaX==0. )  // skip for surface segemnts
         LowerBound = 0.;
       else {
         mp = MeasurementPoint( float(ix), 0.0);
         xLB = topol->localPosition(mp).x();
         LowerBound = 1-calcQ((xLB-CloudCenterX)/SigmaX);
       }

       if(ix == numRows-1 || SigmaX==0. )
         UpperBound = 1.;
       else {
         mp = MeasurementPoint( float(ix+1), 0.0);
         xUB = topol->localPosition(mp).x();
         UpperBound = 1. - calcQ((xUB-CloudCenterX)/SigmaX);
       }

       float   TotalIntegrationRange = UpperBound - LowerBound; // get strip
       x[ix] = TotalIntegrationRange; // save strip integral
       //if(SigmaX==0 || SigmaY==0)
       //cout<<TotalIntegrationRange<<" "<<ix<<std::endl;

     }

    // Now integarte strips in y
    int iy; // TT for compatibility
    for (iy=IPixLeftDownY; iy<=IPixRightUpY; iy++) { //loope over y ind
      float yUB, yLB, UpperBound, LowerBound;

      if(iy == 0 || SigmaY==0.)
        LowerBound = 0.;
      else {
        mp = MeasurementPoint( 0.0, float(iy) );
        yLB = topol->localPosition(mp).y();
        LowerBound = 1. - calcQ((yLB-CloudCenterY)/SigmaY);
      }

      if(iy == numColumns-1 || SigmaY==0. )
        UpperBound = 1.;
      else {
        mp = MeasurementPoint( 0.0, float(iy+1) );
        yUB = topol->localPosition(mp).y();
        UpperBound = 1. - calcQ((yUB-CloudCenterY)/SigmaY);
      }

      float   TotalIntegrationRange = UpperBound - LowerBound;
      y[iy] = TotalIntegrationRange; // save strip integral
      //if(SigmaX==0 || SigmaY==0)
      //cout<<TotalIntegrationRange<<" "<<iy<<std::endl;
    }

    // Get the 2D charge integrals by folding x and y strips
    int chan;
    for (ix=IPixLeftDownX; ix<=IPixRightUpX; ix++) {  // loop over x index
      for (iy=IPixLeftDownY; iy<=IPixRightUpY; iy++) { //loope over y ind

        float ChargeFraction = Charge*x[ix]*y[iy];

        if( ChargeFraction > 0. ) {
          chan = PixelDigi::pixelToChannel( ix, iy);  // Get index
          // Load the amplitude
          hit_signal[chan] += ChargeFraction;
        } // endif


        mp = MeasurementPoint( float(ix), float(iy) );
        LocalPoint lp = topol->localPosition(mp);
        chan = topol->channel(lp);

#ifdef TP_DEBUG
        LogDebug ("Pixel Digitizer")
          << " pixel " << ix << " " << iy << " - "<<" "
          << chan << " " << ChargeFraction<<" "
          << mp.x() << " " << mp.y() <<" "
          << lp.x() << " " << lp.y() << " "  // givex edge position
          << chan; // edge belongs to previous ?
#endif

      } // endfor iy
    } //endfor ix


    // Test conversions (THIS IS FOR TESTING ONLY) comment-out.
    //     mp = topol->measurementPosition( i->position() ); //OK
    //     LocalPoint lp = topol->localPosition(mp);     //OK
    //     std::pair<float,float> p = topol->pixel( i->position() );  //OK
    //     chan = PixelDigi::pixelToChannel( int(p.first), int(p.second));
    //     std::pair<int,int> ip = PixelDigi::channelToPixel(chan);
    //     MeasurementPoint mp1 = MeasurementPoint( float(ip.first),
    //                                       float(ip.second) );
    //     LogDebug ("Pixel Digitizer") << " Test "<< mp.x() << " " << mp.y()
    //                           << " "<< lp.x() << " " << lp.y() << " "<<" "
    //                           <<p.first <<" "<<p.second<<" "<<chan<< " "
    //                           <<" " << ip.first << " " << ip.second << " "
    //                           << mp1.x() << " " << mp1.y() << " " //OK
    //                           << topol->localPosition(mp1).x() << " "  //OK
    //                           << topol->localPosition(mp1).y() << " "
    //                           << topol->channel( i->position() ); //OK


  } // loop over charge distributions

  // Fill the global map with all hit pixels from this event

  for ( hit_map_type::const_iterator im = hit_signal.begin();
        im != hit_signal.end(); ++im) {
    int chan =  (*im).first;
    theSignal[chan] += (makeDigiSimLinks_ ? Amplitude( (*im).second, &hit, (*im).second) : Amplitude( (*im).second, (*im).second) )  ;

#ifdef TP_DEBUG
    std::pair<int,int> ip = PixelDigi::channelToPixel(chan);
    LogDebug ("Pixel Digitizer")
      << " pixel " << ip.first << " " << ip.second << " "
      << theSignal[chan];
#endif
  }

} // end induce_signal

void SiPixelDigitizerAlgorithm::init ( const edm::EventSetup & es )

Definition at line 107 of file SiPixelDigitizerAlgorithm.cc.

References edm::EventSetup::get().

                                                            {
  if(use_ineff_from_db_){// load gain calibration service fromdb...
    theSiPixelGainCalibrationService_->setESObjects( es );
  }
  if(use_deadmodule_DB_) {
    es.get<SiPixelQualityRcd>().get(SiPixelBadModule_);
  }
  if(use_LorentzAngle_DB_) {
    // Get Lorentz angle from DB record
    es.get<SiPixelLorentzAngleSimRcd>().get(SiPixelLorentzAngle_);
  }
  //gets the map and geometry from the DB (to kill ROCs)
  es.get<SiPixelFedCablingMapRcd>().get(map_);
  es.get<TrackerDigiGeometryRecord>().get(geom_);
}

std::map< int, SiPixelDigitizerAlgorithm::CalParameters, std::less< int > > SiPixelDigitizerAlgorithm::initCal ( ) const [private]

Definition at line 269 of file SiPixelDigitizerAlgorithm.cc.

References calmap, dtNoiseDBValidation_cfg::cerr, PixelIndices::channelToPixelROC(), gather_cfg::cout, lut2db_cfg::filename, i, recoMuon::in, create_public_lumi_plots::in_file, geometryCSVtoXML::line, LogDebug, mergeVDriftHistosByStation::name, AlCaHLTBitMon_ParallelJobs::p, SiPixelDigitizerAlgorithm::CalParameters::p0, SiPixelDigitizerAlgorithm::CalParameters::p1, SiPixelDigitizerAlgorithm::CalParameters::p2, SiPixelDigitizerAlgorithm::CalParameters::p3, fitWZ::par0, PixelIndices::pixelToChannelROC(), and AlCaHLTBitMon_QueryRunRegistry::string.

                                         {

  using std::cerr;
  using std::cout;
  using std::endl;

  std::map<int, SiPixelDigitizerAlgorithm::CalParameters, std::less<int> > calmap;
  // Prepare for the analog amplitude miss-calibration
  LogDebug ("PixelDigitizer ")
    << " miss-calibrate the pixel amplitude ";

  const bool ReadCalParameters = false;
  if(ReadCalParameters) {   // Read the calibration files from file
    // read the calibration constants from a file (testing only)
    std::ifstream in_file;  // data file pointer
    char filename[80] = "phCalibrationFit_C0.dat";

    in_file.open(filename, std::ios::in ); // in C++
    if(in_file.bad()) {
      cout << " File not found " << endl;
      return calmap; // signal error
    }
    cout << " file opened : " << filename << endl;

    char line[500];
    for (int i = 0; i < 3; i++) {
      in_file.getline(line, 500,'\n');
      cout<<line<<endl;
    }

    cout << " test map" << endl;

    float par0,par1,par2,par3;
    int colid,rowid;
    std::string name;
    // Read MC tracks
    for(int i=0;i<(52*80);i++)  { // loop over tracks
      in_file >> par0 >> par1 >> par2 >> par3 >> name >> colid >> rowid;
      if(in_file.bad()) { // check for errors
        cerr << "Cannot read data file" << endl;
        return calmap;
      }
      if( in_file.eof() != 0 ) {
        cerr << in_file.eof() << " " << in_file.gcount() << " "
             << in_file.fail() << " " << in_file.good() << " end of file "
             << endl;
        return calmap;
      }

      //cout << " line " << i << " " <<par0<<" "<<par1<<" "<<par2<<" "<<par3<<" "
      //   <<colid<<" "<<rowid<<endl;

      CalParameters onePix;
      onePix.p0=par0;
      onePix.p1=par1;
      onePix.p2=par2;
      onePix.p3=par3;

      // Convert ROC pixel index to channel
      int chan = PixelIndices::pixelToChannelROC(rowid,colid);
      calmap.insert(std::pair<int,CalParameters>(chan,onePix));

      // Testing the index conversion, can be skipped
      std::pair<int,int> p = PixelIndices::channelToPixelROC(chan);
      if(rowid!=p.first) cout<<" wrong channel row "<<rowid<<" "<<p.first<<endl;
      if(colid!=p.second) cout<<" wrong channel col "<<colid<<" "<<p.second<<endl;

    } // pixel loop in a ROC

    cout << " map size  " << calmap.size() <<" max "<<calmap.max_size() << " "
         <<calmap.empty()<< endl;

//     cout << " map size  " << calmap.size()  << endl;
//     map<int,CalParameters,std::less<int> >::iterator ix,it;
//     map<int,CalParameters,std::less<int> >::const_iterator ip;
//     for (ix = calmap.begin(); ix != calmap.end(); ++ix) {
//       int i = (*ix).first;
//       std::pair<int,int> p = channelToPixelROC(i);
//       it  = calmap.find(i);
//       CalParameters y  = (*it).second;
//       CalParameters z = (*ix).second;
//       cout << i <<" "<<p.first<<" "<<p.second<<" "<<y.p0<<" "<<z.p0<<" "<<calmap[i].p0<<endl;

//       //int dummy=0;
//       //cin>>dummy;
//     }

  } // end if readparameters
  return calmap;
} // end initCal()

void SiPixelDigitizerAlgorithm::initializeEvent ( ) [inline]

Definition at line 52 of file SiPixelDigitizerAlgorithm.h.

References _signal.

                         {
    _signal.clear();
  }

void SiPixelDigitizerAlgorithm::make_digis	(	float	thePixelThresholdInE,
		uint32_t	detID,
		std::vector< PixelDigi > &	digis,
		std::vector< PixelDigiSimLink > &	simlinks,
		const TrackerTopology *	tTopo
	)		const `[private]`

Definition at line 1047 of file SiPixelDigitizerAlgorithm.cc.

References _signal, ecalMGPA::adc(), PixelDigi::channelToPixel(), cuy::col, doMissCalibrate, i, LogDebug, makeDigiSimLinks_, min, missCalibrate(), PixelSubdetector::PixelBarrel, TrackerTopology::pxbLayer(), DetId::subdetId(), theAdcFullScale, theAdcFullScaleStack, theElectronPerADC, theFirstStackLayer, theThresholdInE_BPix, theThresholdInE_BPix_L1, and theThresholdInE_FPix.

Referenced by digitize().

                                                                                {

#ifdef TP_DEBUG
  LogDebug ("Pixel Digitizer") << " make digis "<<" "
                               << " pixel threshold FPix" << theThresholdInE_FPix << " "
                               << " pixel threshold BPix" << theThresholdInE_BPix << " "
                               << " pixel threshold BPix Layer1" << theThresholdInE_BPix_L1 << " "
                               << " List pixels passing threshold ";
#endif

  // Loop over hit pixels

  signalMaps::const_iterator it = _signal.find(detID);
  if (it == _signal.end()) {
    return;
  }

  const signal_map_type& theSignal = (*it).second;

  for (signal_map_const_iterator i = theSignal.begin(); i != theSignal.end(); ++i) {

    float signalInElectrons = (*i).second ;   // signal in electrons

    // Do the miss calibration for calibration studies only.
    //if(doMissCalibrate) signalInElectrons = missCalibrate(signalInElectrons)

    // Do only for pixels above threshold

    if( signalInElectrons >= thePixelThresholdInE) { // check threshold

      int chan =  (*i).first;  // channel number
      std::pair<int,int> ip = PixelDigi::channelToPixel(chan);
      int adc=0;  // ADC count as integer

      // Do the miss calibration for calibration studies only.
      if(doMissCalibrate) {
        int row = ip.first;  // X in row
        int col = ip.second; // Y is in col
        adc = int(missCalibrate(detID, col, row, signalInElectrons)); //full misscalib.
      } else { // Just do a simple electron->adc conversion
        adc = int( signalInElectrons / theElectronPerADC ); // calibrate gain
      }
      adc = std::min(adc, theAdcFullScale); // Check maximum value
// Calculate layerIndex
     if (theAdcFullScale!=theAdcFullScaleStack){
        unsigned int Sub_detid=DetId(detID).subdetId();
        if(Sub_detid == PixelSubdetector::PixelBarrel){ // Barrel modules
          int lay = tTopo->pxbLayer(detID);
          if (lay>=theFirstStackLayer) {
            // Set to 1 if over the threshold
            if (theAdcFullScaleStack==1) {adc=1;}
            // Make it a linear fit to the full scale of the normal adc count.   Start new adc from 1 not zero.
            if (theAdcFullScaleStack!=1&&theAdcFullScaleStack!=theAdcFullScale) {adc = int (1 + adc * (theAdcFullScaleStack-1)/float(theAdcFullScale) );}
          }
        }
     } // Only enter this if the Adc changes for the outer layers
#ifdef TP_DEBUG
      LogDebug ("Pixel Digitizer")
        << (*i).first << " " << (*i).second << " " << signalInElectrons
        << " " << adc << ip.first << " " << ip.second ;
#endif

      // Load digis
      digis.emplace_back(ip.first, ip.second, adc);

      if (makeDigiSimLinks_ && (*i).second.hitInfo()!=0) {
        //digilink
        if((*i).second.trackIds().size()>0){
          simlink_map simi;
          unsigned int il=0;
          for( std::vector<unsigned int>::const_iterator itid = (*i).second.trackIds().begin();
               itid != (*i).second.trackIds().end(); ++itid) {
            simi[*itid].push_back((*i).second.individualampl()[il]);
            il++;
          }

          //sum the contribution of the same trackid
          for( simlink_map::iterator simiiter=simi.begin();
               simiiter!=simi.end();
               simiiter++){

            float sum_samechannel=0;
            for (unsigned int iii=0;iii<(*simiiter).second.size();iii++){
              sum_samechannel+=(*simiiter).second[iii];
            }
            float fraction=sum_samechannel/(*i).second;
            if(fraction>1.) fraction=1.;
            simlinks.emplace_back((*i).first, (*simiiter).first, (*i).second.eventId(), fraction);
          }
        }
      }
    }
  }
}

float SiPixelDigitizerAlgorithm::missCalibrate	(	uint32_t	detID,
		int	col,
		int	row,
		float	amp
	)		const `[private]`

Definition at line 1391 of file SiPixelDigitizerAlgorithm.cc.

References BPix_p0, BPix_p1, BPix_p2, BPix_p3, electronsPerVCAL, electronsPerVCAL_Offset, FPix_p0, FPix_p1, FPix_p2, FPix_p3, p1, p2, p3, PixelSubdetector::PixelBarrel, and DetId::subdetId().

Referenced by make_digis().

                                                                      {
  // Central values
  //const float p0=0.00352, p1=0.868, p2=112., p3=113.; // pix(0,0,0)
  //  const float p0=0.00382, p1=0.886, p2=112.7, p3=113.0; // average roc=0
  //const float p0=0.00492, p1=1.998, p2=90.6, p3=134.1; // average roc=6
  // Smeared (rms)
  //const float s0=0.00020, s1=0.051, s2=5.4, s3=4.4; // average roc=0
  //const float s0=0.00015, s1=0.043, s2=3.2, s3=3.1; // col average roc=0

  // Make 2 sets of parameters for Fpix and BPIx:

  float p0=0.0;
  float p1=0.0;
  float p2=0.0;
  float p3=0.0;

  unsigned int Sub_detid=DetId(detID).subdetId();

    if(Sub_detid == PixelSubdetector::PixelBarrel){// barrel layers
      p0 = BPix_p0;
      p1 = BPix_p1;
      p2 = BPix_p2;
      p3 = BPix_p3;
    } else {// forward disks
      p0 = FPix_p0;
      p1 = FPix_p1;
      p2 = FPix_p2;
      p3 = FPix_p3;
    }

  //  const float electronsPerVCAL = 65.5; // our present VCAL calibration (feb 2009)
  //  const float electronsPerVCAL_Offset = -414.0; // our present VCAL calibration (feb 2009)
  float newAmp = 0.; //Modified signal

  // Convert electrons to VCAL units
  float signal = (signalInElectrons-electronsPerVCAL_Offset)/electronsPerVCAL;

  // Simulate the analog response with fixed parametrization
  newAmp = p3 + p2 * tanh(p0*signal - p1);


  // Use the pixel-by-pixel calibrations
  //transform to ROC index coordinates
  //int chipIndex=0, colROC=0, rowROC=0;
  //std::auto_ptr<PixelIndices> pIndexConverter(new PixelIndices(numColumns,numRows));
  //pIndexConverter->transformToROC(col,row,chipIndex,colROC,rowROC);

  // Use calibration from a file
  //int chanROC = PixelIndices::pixelToChannelROC(rowROC,colROC); // use ROC coordinates
  //float pp0=0, pp1=0,pp2=0,pp3=0;
  //map<int,CalParameters,std::less<int> >::const_iterator it=calmap.find(chanROC);
  //CalParameters y  = (*it).second;
  //pp0 = y.p0;
  //pp1 = y.p1;
  //pp2 = y.p2;
  //pp3 = y.p3;

  //
  // Use random smearing
  // Randomize the pixel response
  //float pp0  = RandGaussQ::shoot(p0,s0);
  //float pp1  = RandGaussQ::shoot(p1,s1);
  //float pp2  = RandGaussQ::shoot(p2,s2);
  //float pp3  = RandGaussQ::shoot(p3,s3);

  //newAmp = pp3 + pp2 * tanh(pp0*signal - pp1); // Final signal

  //cout<<" misscalibrate "<<col<<" "<<row<<" "<<chipIndex<<" "<<colROC<<" "
  //  <<rowROC<<" "<<signalInElectrons<<" "<<signal<<" "<<newAmp<<" "
  //  <<(signalInElectrons/theElectronPerADC)<<std::endl;

  return newAmp;
}

void SiPixelDigitizerAlgorithm::module_killing_conf ( uint32_t detID ) [private]

Definition at line 1572 of file SiPixelDigitizerAlgorithm.cc.

References _signal, PixelDigi::channelToPixel(), DeadModules, cond::rpcobgas::detid, i, Module, AlCaHLTBitMon_QueryRunRegistry::string, and use_module_killing_.

Referenced by digitize().

                                                                  {
  if(!use_module_killing_)
    return;
  
  bool isbad=false;
  
  Parameters::const_iterator itDeadModules=DeadModules.begin();
  
  int detid = detID;
  for(; itDeadModules != DeadModules.end(); ++itDeadModules){
    int Dead_detID = itDeadModules->getParameter<int>("Dead_detID");
    if(detid == Dead_detID){
      isbad=true;
      break;
    }
  }
  
  if(!isbad)
    return;

  signal_map_type& theSignal = _signal[detID];
  
  std::string Module = itDeadModules->getParameter<std::string>("Module");
  
  if(Module=="whole"){
    for(signal_map_iterator i = theSignal.begin();i != theSignal.end(); ++i) {
      i->second.set(0.); // reset amplitude
    }
  }
  
  for(signal_map_iterator i = theSignal.begin();i != theSignal.end(); ++i) {
    std::pair<int,int> ip = PixelDigi::channelToPixel(i->first);//get pixel pos

    if(Module=="tbmA" && ip.first>=80 && ip.first<=159){
      i->second.set(0.);
    }

    if( Module=="tbmB" && ip.first<=79){
      i->second.set(0.);
    }
  }
}

void SiPixelDigitizerAlgorithm::module_killing_DB ( uint32_t detID ) [private]

Definition at line 1615 of file SiPixelDigitizerAlgorithm.cc.

References _signal, PixelDigi::channelToPixel(), SiPixelQuality::disabledModuleType::errorType, i, sipixelobjects::PixelROC::idInDetUnit(), j, map_, getHLTPrescaleColumns::path, edm::ESHandle< T >::product(), SiPixelBadModule_, sipixelobjects::PixelROC::toGlobal(), and use_module_killing_.

Referenced by digitize().

                                                                {
// Not SLHC safe for now
  if(!use_module_killing_)
    return;
  
  bool isbad=false;
  
  std::vector<SiPixelQuality::disabledModuleType>disabledModules = SiPixelBadModule_->getBadComponentList();
  
  SiPixelQuality::disabledModuleType badmodule;
  
  for (size_t id=0;id<disabledModules.size();id++)
    {
      if(detID==disabledModules[id].DetID){
        isbad=true;
        badmodule = disabledModules[id];
        break;
      }
    }
  
  if(!isbad)
    return;

  signal_map_type& theSignal = _signal[detID];
  
  //std::cout<<"Hit in: "<< detID <<" errorType "<< badmodule.errorType<<" BadRocs="<<std::hex<<SiPixelBadModule_->getBadRocs(detID)<<dec<<" "<<std::endl;
  if(badmodule.errorType == 0){ // this is a whole dead module.
    
    for(signal_map_iterator i = theSignal.begin();i != theSignal.end(); ++i) {
      i->second.set(0.); // reset amplitude
    }
  }
  else { // all other module types: half-modules and single ROCs.
    // Get Bad ROC position:
    //follow the example of getBadRocPositions in CondFormats/SiPixelObjects/src/SiPixelQuality.cc
    std::vector<GlobalPixel> badrocpositions (0);
    for(unsigned int j = 0; j < 16; j++){
      if(SiPixelBadModule_->IsRocBad(detID, j) == true){
        
        std::vector<CablingPathToDetUnit> path = map_.product()->pathToDetUnit(detID);
        typedef  std::vector<CablingPathToDetUnit>::const_iterator IT;
        for  (IT it = path.begin(); it != path.end(); ++it) {
          const PixelROC* myroc = map_.product()->findItem(*it);
          if( myroc->idInDetUnit() == j) {
            LocalPixel::RocRowCol  local = { 39, 25};   //corresponding to center of ROC row, col
            GlobalPixel global = myroc->toGlobal( LocalPixel(local) );
            badrocpositions.push_back(global);
            break;
          }
        }
      }
    }// end of getBadRocPositions
    
    
    for(signal_map_iterator i = theSignal.begin();i != theSignal.end(); ++i) {
      std::pair<int,int> ip = PixelDigi::channelToPixel(i->first);//get pixel pos
      
      for(std::vector<GlobalPixel>::const_iterator it = badrocpositions.begin(); it != badrocpositions.end(); ++it){
        if(it->row >= 80 && ip.first >= 80 ){
          if((fabs(ip.second - it->col) < 26) ) {i->second.set(0.);}
          else if(it->row==120 && ip.second-it->col==26){i->second.set(0.);}
          else if(it->row==119 && it->col-ip.second==26){i->second.set(0.);}
        }
        else if(it->row < 80 && ip.first < 80 ){
          if((fabs(ip.second - it->col) < 26) ){i->second.set(0.);}
          else if(it->row==40 && ip.second-it->col==26){i->second.set(0.);}
          else if(it->row==39 && it->col-ip.second==26){i->second.set(0.);}
       }
      }
    }
  }
}

void SiPixelDigitizerAlgorithm::pixel_inefficiency	(	const PixelEfficiencies &	eff,
		const PixelGeomDetUnit *	pixdet,
		const TrackerTopology *	tTopo
	)		`[private]`

Definition at line 1260 of file SiPixelDigitizerAlgorithm.cc.

References _signal, PixelDigi::channelToPixel(), cuy::col, flatDistribution_, SiPixelDigitizerAlgorithm::PixelEfficiencies::FPixIndex, GeomDet::geographicalId(), i, LogDebug, PixelTopology::ncolumns(), PixelTopology::nrows(), NumberOfBarrelLayers, PixelSubdetector::PixelBarrel, TrackerTopology::pxbLayer(), TrackerTopology::pxfDisk(), DetId::rawId(), PixelGeomDetUnit::specificTopology(), DetId::subdetId(), SiPixelDigitizerAlgorithm::PixelEfficiencies::thePixelChipEfficiency, SiPixelDigitizerAlgorithm::PixelEfficiencies::thePixelColEfficiency, and SiPixelDigitizerAlgorithm::PixelEfficiencies::thePixelEfficiency.

Referenced by digitize().

                                                                                 {

  uint32_t detID= pixdet->geographicalId().rawId();

  signal_map_type& theSignal = _signal[detID];

  const PixelTopology* topol=&pixdet->specificTopology();
  int numColumns = topol->ncolumns();  // det module number of cols&rows
  int numRows = topol->nrows();

  // Predefined efficiencies
  float pixelEfficiency  = 1.0;
  float columnEfficiency = 1.0;
  float chipEfficiency   = 1.0;

  // setup the chip indices conversion
  unsigned int Subid=DetId(detID).subdetId();
  if    (Subid==  PixelSubdetector::PixelBarrel){// barrel layers
    int layerIndex=tTopo->pxbLayer(detID);
    pixelEfficiency  = eff.thePixelEfficiency[layerIndex-1];
    columnEfficiency = eff.thePixelColEfficiency[layerIndex-1];
    chipEfficiency   = eff.thePixelChipEfficiency[layerIndex-1];
    //std::cout <<"Using BPix columnEfficiency = "<<columnEfficiency<< " for layer = "<<layerIndex <<"\n";
    // This should never happen, but only check if it is not an upgrade geometry
    if (NumberOfBarrelLayers==3){
       if(numColumns>416)  LogWarning ("Pixel Geometry") <<" wrong columns in barrel "<<numColumns;
       if(numRows>160)  LogWarning ("Pixel Geometry") <<" wrong rows in barrel "<<numRows;
    }
  } else {                // forward disks
    unsigned int diskIndex=tTopo->pxfDisk(detID)+eff.FPixIndex; // Use diskIndex-1 later to stay consistent with BPix
    //if (eff.FPixIndex>diskIndex-1){throw cms::Exception("Configuration") <<"SiPixelDigitizer is using the wrong efficiency value. index = "
    //                                                                       <<diskIndex-1<<" , MinIndex = "<<eff.FPixIndex<<" ... "<<tTopo->pxfDisk(detID);}
    pixelEfficiency  = eff.thePixelEfficiency[diskIndex-1];
    columnEfficiency = eff.thePixelColEfficiency[diskIndex-1];
    chipEfficiency   = eff.thePixelChipEfficiency[diskIndex-1];
    //std::cout <<"Using FPix columnEfficiency = "<<columnEfficiency<<" for Disk = "<< tTopo->pxfDisk(detID)<<"\n";
    // Sometimes the forward pixels have wrong size,
    // this crashes the index conversion, so exit, but only check if it is not an upgrade geometry
    if (NumberOfBarrelLayers==3){
       if(numColumns>260 || numRows>160) {
         if(numColumns>260)  LogWarning ("Pixel Geometry") <<" wrong columns in endcaps "<<numColumns;
         if(numRows>160)  LogWarning ("Pixel Geometry") <<" wrong rows in endcaps "<<numRows;
         return;
       }
    }
  } // if barrel/forward

#ifdef TP_DEBUG
  LogDebug ("Pixel Digitizer") << " enter pixel_inefficiency " << pixelEfficiency << " "
                               << columnEfficiency << " " << chipEfficiency;
#endif

  // Initilize the index converter
  //PixelIndices indexConverter(numColumns,numRows);
  std::auto_ptr<PixelIndices> pIndexConverter(new PixelIndices(numColumns,numRows));

  int chipIndex = 0;
  int rowROC = 0;
  int colROC = 0;
  std::map<int, int, std::less<int> >chips, columns;
  std::map<int, int, std::less<int> >::iterator iter;

  // Find out the number of columns and rocs hits
  // Loop over hit pixels, amplitude in electrons, channel = coded row,col
  for (signal_map_const_iterator i = theSignal.begin(); i != theSignal.end(); ++i) {

    int chan = i->first;
    std::pair<int,int> ip = PixelDigi::channelToPixel(chan);
    int row = ip.first;  // X in row
    int col = ip.second; // Y is in col
    //transform to ROC index coordinates
    pIndexConverter->transformToROC(col,row,chipIndex,colROC,rowROC);
    int dColInChip = pIndexConverter->DColumn(colROC); // get ROC dcol from ROC col
    //dcol in mod
    int dColInDet = pIndexConverter->DColumnInModule(dColInChip,chipIndex);

    chips[chipIndex]++;
    columns[dColInDet]++;
  }

  // Delete some ROC hits.
  for ( iter = chips.begin(); iter != chips.end() ; iter++ ) {
    //float rand  = RandFlat::shoot();
    float rand  = flatDistribution_->fire();
    if( rand > chipEfficiency ) chips[iter->first]=0;
  }

  // Delete some Dcol hits.
  for ( iter = columns.begin(); iter != columns.end() ; iter++ ) {
    //float rand  = RandFlat::shoot();
    float rand  = flatDistribution_->fire();
    if( rand > columnEfficiency ) columns[iter->first]=0;
  }

  // Now loop again over pixels to kill some of them.
  // Loop over hit pixels, amplitude in electrons, channel = coded row,col
  for(signal_map_iterator i = theSignal.begin();i != theSignal.end(); ++i) {

    //    int chan = i->first;
    std::pair<int,int> ip = PixelDigi::channelToPixel(i->first);//get pixel pos
    int row = ip.first;  // X in row
    int col = ip.second; // Y is in col
    //transform to ROC index coordinates
    pIndexConverter->transformToROC(col,row,chipIndex,colROC,rowROC);
    int dColInChip = pIndexConverter->DColumn(colROC); //get ROC dcol from ROC col
    //dcol in mod
    int dColInDet = pIndexConverter->DColumnInModule(dColInChip,chipIndex);


    //float rand  = RandFlat::shoot();
    float rand  = flatDistribution_->fire();
    if( chips[chipIndex]==0 || columns[dColInDet]==0
        || rand>pixelEfficiency ) {
      // make pixel amplitude =0, pixel will be lost at clusterization
      i->second.set(0.); // reset amplitude,
    } // end if

  } // end pixel loop
} // end pixel_indefficiency

void SiPixelDigitizerAlgorithm::pixel_inefficiency_db ( uint32_t detID ) [private]

Definition at line 1547 of file SiPixelDigitizerAlgorithm.cc.

References _signal, PixelDigi::channelToPixel(), cuy::col, i, theSiPixelGainCalibrationService_, and use_ineff_from_db_.

Referenced by digitize().

                                                                    {
  if(!use_ineff_from_db_)
    return;
  
  signal_map_type& theSignal = _signal[detID];

  // Loop over hit pixels, amplitude in electrons, channel = coded row,col
  for(signal_map_iterator i = theSignal.begin();i != theSignal.end(); ++i) {

    //    int chan = i->first;
    std::pair<int,int> ip = PixelDigi::channelToPixel(i->first);//get pixel pos
    int row = ip.first;  // X in row
    int col = ip.second; // Y is in col
    //transform to ROC index coordinates
    if(theSiPixelGainCalibrationService_->isDead(detID, col, row)){
      //      std::cout << "now in isdead check, row " << detID << " " << col << "," << row << std::std::endl;
      // make pixel amplitude =0, pixel will be lost at clusterization
      i->second.set(0.); // reset amplitude,
    } // end if
  } // end pixel loop
} // end pixel_indefficiency

void SiPixelDigitizerAlgorithm::primary_ionization	(	const PSimHit &	hit,
		std::vector< EnergyDepositUnit > &	ionization_points
	)		const `[private]`

Definition at line 560 of file SiPixelDigitizerAlgorithm.cc.

References relval_parameters_module::energy, PSimHit::energyLoss(), PSimHit::entryPoint(), PSimHit::exitPoint(), fluctuateCharge, fluctuateEloss(), GeVperElectron, i, LogDebug, PV3DBase< T, PVType, FrameType >::mag(), PSimHit::pabs(), PSimHit::particleType(), evf::utils::pid, point, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by accumulateSimHits().

                                                                                                                            {

// Straight line approximation for trajectory inside active media

  const float SegmentLength = 0.0010; //10microns in cm
  float energy;

  // Get the 3D segment direction vector
  LocalVector direction = hit.exitPoint() - hit.entryPoint();

  float eLoss = hit.energyLoss();  // Eloss in GeV
  float length = direction.mag();  // Track length in Silicon

  int NumberOfSegments = int ( length / SegmentLength); // Number of segments
  if(NumberOfSegments < 1) NumberOfSegments = 1;

#ifdef TP_DEBUG
  LogDebug ("Pixel Digitizer")
    << " enter primary_ionzation " << NumberOfSegments
    << " shift = "
    << (hit.exitPoint().x()-hit.entryPoint().x()) << " "
    << (hit.exitPoint().y()-hit.entryPoint().y()) << " "
    << (hit.exitPoint().z()-hit.entryPoint().z()) << " "
    << hit.particleType() <<" "<< hit.pabs() ;
#endif

  float* elossVector = new float[NumberOfSegments];  // Eloss vector

  if( fluctuateCharge ) {
    //MP DA RIMUOVERE ASSOLUTAMENTE
    int pid = hit.particleType();
    //int pid=211;  // assume it is a pion

    float momentum = hit.pabs();
    // Generate fluctuated charge points
    fluctuateEloss(pid, momentum, eLoss, length, NumberOfSegments,
                   elossVector);
  }

  ionization_points.resize( NumberOfSegments); // set size

  // loop over segments
  for ( int i = 0; i != NumberOfSegments; i++) {
    // Divide the segment into equal length subsegments
    Local3DPoint point = hit.entryPoint() +
      float((i+0.5)/NumberOfSegments) * direction;

    if( fluctuateCharge )
      energy = elossVector[i]/GeVperElectron; // Convert charge to elec.
    else
      energy = hit.energyLoss()/GeVperElectron/float(NumberOfSegments);

    EnergyDepositUnit edu( energy, point); //define position,energy point
    ionization_points[i] = edu; // save

#ifdef TP_DEBUG
    LogDebug ("Pixel Digitizer")
      << i << " " << ionization_points[i].x() << " "
      << ionization_points[i].y() << " "
      << ionization_points[i].z() << " "
      << ionization_points[i].energy();
#endif

  }  // end for loop

  delete[] elossVector;

}