SiPixelGaussianSmearingRecHitConverterAlgorithm Class Reference

#include <SiPixelGaussianSmearingRecHitConverterAlgorithm.h>

Public Member Functions
LocalError	getError ()
double	getErrorX ()
double	getErrorY ()
double	getErrorZ ()
unsigned int	getPixelMultiplicityAlpha ()
unsigned int	getPixelMultiplicityBeta ()
Local3DPoint	getPosition ()
double	getPositionX ()
double	getPositionY ()
double	getPositionZ ()
	SiPixelGaussianSmearingRecHitConverterAlgorithm (const edm::ParameterSet &pset, GeomDetType::SubDetector pixelPart)
void	smearHit (const PSimHit &simHit, const PixelGeomDetUnit *detUnit, const double boundX, const double boundY, RandomEngineAndDistribution const *)
virtual	~SiPixelGaussianSmearingRecHitConverterAlgorithm ()

Private Member Functions
void	initializeBarrel ()
void	initializeForward ()
bool	isFlipped (const PixelGeomDetUnit *theDet) const

Private Attributes
bool	isForward
edm::ParameterSet	pset_
double	rescotAlpha_binMin
unsigned int	rescotAlpha_binN
double	rescotAlpha_binWidth
double	rescotBeta_binMin
unsigned int	rescotBeta_binN
double	rescotBeta_binWidth
int	resqbin_binMin
unsigned int	resqbin_binN
int	resqbin_binWidth
int	tempId
unsigned int	theClslenx
unsigned int	theClsleny
LocalError	theError
double	theErrorX
double	theErrorY
double	theErrorZ
unsigned int	theLayer
GeomDetType::SubDetector	thePixelPart
TFile *	thePixelResolutionFile1
TFile *	thePixelResolutionFile2
std::string	thePixelResolutionFileName1
std::string	thePixelResolutionFileName2
std::vector< SiPixelTemplateStore >	thePixelTemp_
Local3DPoint	thePosition
double	thePositionX
double	thePositionY
double	thePositionZ
std::map< unsigned, const SimpleHistogramGenerator * >	theXHistos
std::map< unsigned, const SimpleHistogramGenerator * >	theYHistos
bool	useCMSSWPixelParameterization

Detailed Description

Definition at line 36 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

Constructor & Destructor Documentation

SiPixelGaussianSmearingRecHitConverterAlgorithm::SiPixelGaussianSmearingRecHitConverterAlgorithm ( const edm::ParameterSet &  pset, GeomDetType::SubDetector  pixelPart  )

explicit

Definition at line 40 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.cc.

References gather_cfg::cout, edm::hlt::Exception, contentValuesFiles::fullPath, edm::ParameterSet::getParameter(), initializeBarrel(), initializeForward(), isForward, GeomDetEnumerators::PixelBarrel, GeomDetEnumerators::PixelEndcap, pset_, SiPixelTemplate::pushfile(), tempId, thePixelPart, thePixelResolutionFile1, thePixelResolutionFile2, thePixelResolutionFileName1, thePixelResolutionFileName2, thePixelTemp_, theXHistos, theYHistos, and useCMSSWPixelParameterization.

:
  pset_(pset),
  thePixelPart(pixelPart)
{
  // Switch between old (ORCA) and new (CMSSW) pixel parameterization
  useCMSSWPixelParameterization = pset.getParameter<bool>("UseCMSSWPixelParametrization");

  thePixelResolutionFile1=0;
  thePixelResolutionFile2=0;

  if( thePixelPart == GeomDetEnumerators::PixelBarrel ) {
     isForward = false;
     thePixelResolutionFileName1 = pset_.getParameter<string>( "NewPixelBarrelResolutionFile1" );
     thePixelResolutionFile1 = new 
       TFile( edm::FileInPath( thePixelResolutionFileName1 ).fullPath().c_str()  ,"READ");
     thePixelResolutionFileName2 =  pset_.getParameter<string>( "NewPixelBarrelResolutionFile2" );
     thePixelResolutionFile2 = new
       TFile( edm::FileInPath( thePixelResolutionFileName2 ).fullPath().c_str()  ,"READ");
     initializeBarrel();
     tempId = pset_.getParameter<int> ( "templateIdBarrel" );
     if( ! SiPixelTemplate::pushfile(tempId, thePixelTemp_) )
         throw cms::Exception("SiPixelGaussianSmearingRecHitConverterAlgorithm:")
        <<"SiPixel Barrel Template Not Loaded Correctly!"<<endl;
#ifdef  FAMOS_DEBUG
     cout<<"The Barrel map size is "<<theXHistos.size()<<" and "<<theYHistos.size()<<endl;
#endif
  }
  else
  if ( thePixelPart == GeomDetEnumerators::PixelEndcap ) {
     isForward = true;
     thePixelResolutionFileName1 = pset_.getParameter<string>( "NewPixelForwardResolutionFile" );
     thePixelResolutionFile1 = new
       TFile( edm::FileInPath( thePixelResolutionFileName1 ).fullPath().c_str()  ,"READ");
     initializeForward();
     tempId = pset_.getParameter<int> ( "templateIdForward" );
     if( ! SiPixelTemplate::pushfile(tempId, thePixelTemp_) )
         throw cms::Exception("SiPixelGaussianSmearingRecHitConverterAlgorithm:")
            <<"SiPixel Forward Template Not Loaded Correctly!"<<endl;
#ifdef  FAMOS_DEBUG
     cout<<"The Forward map size is "<<theXHistos.size()<<" and "<<theYHistos.size()<<endl;
#endif
  }
  else
     throw cms::Exception("SiPixelGaussianSmearingRecHitConverterAlgorithm :")
       <<"Not a pixel detector"<<endl;

  if ( thePixelResolutionFile2) {
    thePixelResolutionFile2->Close();
    delete thePixelResolutionFile2;
  }
  if ( thePixelResolutionFile1) {
    thePixelResolutionFile1->Close();
    delete thePixelResolutionFile1;
  }
  thePixelResolutionFile1=0;
  thePixelResolutionFile2=0;
}

SiPixelGaussianSmearingRecHitConverterAlgorithm::~SiPixelGaussianSmearingRecHitConverterAlgorithm ( )

virtual

Definition at line 101 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.cc.

References citk::for(), theXHistos, and theYHistos.

{
  
  std::map<unsigned,const SimpleHistogramGenerator*>::const_iterator it;
  for ( it=theXHistos.begin(); it!=theXHistos.end(); ++it )
    delete it->second;
  for ( it=theYHistos.begin(); it!=theYHistos.end(); ++it )
    delete it->second;

  theXHistos.clear();
  theYHistos.clear();

}

Member Function Documentation

LocalError SiPixelGaussianSmearingRecHitConverterAlgorithm::getError ( )

inline

Definition at line 51 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

References theError.

Referenced by SiTrackerGaussianSmearingRecHitConverter::gaussianSmearing().

{return theError;}

double SiPixelGaussianSmearingRecHitConverterAlgorithm::getErrorX ( )

inline

Definition at line 52 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

References theErrorX.

{return theErrorX;}

double SiPixelGaussianSmearingRecHitConverterAlgorithm::getErrorY ( )

inline

Definition at line 53 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

References theErrorY.

{return theErrorY;}

double SiPixelGaussianSmearingRecHitConverterAlgorithm::getErrorZ ( )

inline

Definition at line 54 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

References theErrorZ.

{return theErrorZ;}

unsigned int SiPixelGaussianSmearingRecHitConverterAlgorithm::getPixelMultiplicityAlpha ( )

inline

Definition at line 55 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

References theClslenx.

Referenced by SiTrackerGaussianSmearingRecHitConverter::gaussianSmearing().

{return theClslenx;}

unsigned int SiPixelGaussianSmearingRecHitConverterAlgorithm::getPixelMultiplicityBeta ( )

inline

Definition at line 56 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

References theClsleny.

Referenced by SiTrackerGaussianSmearingRecHitConverter::gaussianSmearing().

{return theClsleny;}

Local3DPoint SiPixelGaussianSmearingRecHitConverterAlgorithm::getPosition ( )

inline

Definition at line 47 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

References thePosition.

Referenced by SiTrackerGaussianSmearingRecHitConverter::gaussianSmearing().

{return thePosition;}

double SiPixelGaussianSmearingRecHitConverterAlgorithm::getPositionX ( )

inline

Definition at line 48 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

References thePositionX.

{return thePositionX;}

double SiPixelGaussianSmearingRecHitConverterAlgorithm::getPositionY ( )

inline

Definition at line 49 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

References thePositionY.

{return thePositionY;}

double SiPixelGaussianSmearingRecHitConverterAlgorithm::getPositionZ ( )

inline

Definition at line 50 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

References thePositionZ.

{return thePositionZ;}

void SiPixelGaussianSmearingRecHitConverterAlgorithm::initializeBarrel ( )

private

Definition at line 498 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.cc.

References rescotAlpha_binMin, rescotAlpha_binN, rescotAlpha_binWidth, rescotBeta_binMin, rescotBeta_binN, rescotBeta_binWidth, resqbin_binMin, resqbin_binN, resqbin_binWidth, thePixelResolutionFile1, thePixelResolutionFile2, theXHistos, and theYHistos.

Referenced by SiPixelGaussianSmearingRecHitConverterAlgorithm().

{
  //Hard coded at the moment, can easily be changed to be configurable
  rescotAlpha_binMin = -0.2;
  rescotAlpha_binWidth = 0.08 ;
  rescotAlpha_binN = 5;
  rescotBeta_binMin = -5.5;
  rescotBeta_binWidth = 1.0;
  rescotBeta_binN = 11;
  resqbin_binMin = 0;
  resqbin_binWidth = 1;
  resqbin_binN = 4;

  // Initialize the barrel histos once and for all, and prepare the random generation
  for ( unsigned cotalphaHistBin=1; cotalphaHistBin<=rescotAlpha_binN; ++cotalphaHistBin )
     for ( unsigned cotbetaHistBin=1; cotbetaHistBin<=rescotBeta_binN; ++cotbetaHistBin )  {
         unsigned int singleBigPixelHistN = 1 * 100000
                                        + cotalphaHistBin * 100
                                        + cotbetaHistBin ;
         theXHistos[singleBigPixelHistN] = new SimpleHistogramGenerator(
              (TH1F*) thePixelResolutionFile1->Get(  Form( "DQMData/clustBPIX/hx%u" , singleBigPixelHistN ) ) );
         theYHistos[singleBigPixelHistN] = new SimpleHistogramGenerator(
              (TH1F*) thePixelResolutionFile1->Get(  Form( "DQMData/clustBPIX/hy%u" , singleBigPixelHistN ) ) );
         unsigned int singlePixelHistN = 1 * 100000 + 1 * 1000
                                        + cotalphaHistBin * 100
                                        + cotbetaHistBin ;
         theXHistos[singlePixelHistN] = new SimpleHistogramGenerator(
              (TH1F*) thePixelResolutionFile1->Get(  Form( "DQMData/clustBPIX/hx%u" , singlePixelHistN ) ) );
         theYHistos[singlePixelHistN] = new SimpleHistogramGenerator(
              (TH1F*) thePixelResolutionFile1->Get(  Form( "DQMData/clustBPIX/hy%u" , singlePixelHistN ) ) );
         for( unsigned qbinBin=1;  qbinBin<=resqbin_binN; ++qbinBin )  {
             unsigned int edgePixelHistN = cotalphaHistBin * 1000
                                        +  cotbetaHistBin * 10
                                        +  qbinBin;
             theXHistos[edgePixelHistN] = new SimpleHistogramGenerator(
              (TH1F*) thePixelResolutionFile2->Get(  Form( "DQMData/clustBPIX/hx0%u" ,edgePixelHistN ) ) );
             theYHistos[edgePixelHistN] = new SimpleHistogramGenerator(
              (TH1F*) thePixelResolutionFile2->Get(  Form( "DQMData/clustBPIX/hy0%u" ,edgePixelHistN ) ) );
             unsigned int multiPixelBigHistN = 1 * 1000000 + 1 * 100000
                                           + cotalphaHistBin * 1000
                                           + cotbetaHistBin * 10
                                           + qbinBin;
             theXHistos[multiPixelBigHistN] = new SimpleHistogramGenerator(
              (TH1F*) thePixelResolutionFile1->Get(  Form( "DQMData/clustBPIX/hx%u" ,multiPixelBigHistN ) ) );
             theYHistos[multiPixelBigHistN] = new SimpleHistogramGenerator(
              (TH1F*) thePixelResolutionFile1->Get(  Form( "DQMData/clustBPIX/hy%u" ,multiPixelBigHistN ) ) );
             unsigned int multiPixelHistN = 1 * 1000000 + 1 * 100000 + 1 * 10000
                                           + cotalphaHistBin * 1000
                                           + cotbetaHistBin * 10
                                           + qbinBin;
             theXHistos[multiPixelHistN] = new SimpleHistogramGenerator(
             (TH1F*) thePixelResolutionFile1->Get(  Form( "DQMData/clustBPIX/hx%u" , multiPixelHistN ) ) );
             theYHistos[multiPixelHistN] = new SimpleHistogramGenerator(
             (TH1F*) thePixelResolutionFile1->Get(  Form( "DQMData/clustBPIX/hy%u" , multiPixelHistN ) ) );
          } //end for qbinBin
     }//end for cotalphaHistBin, cotbetaHistBin
}

void SiPixelGaussianSmearingRecHitConverterAlgorithm::initializeForward ( )

private

Definition at line 556 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.cc.

References rescotAlpha_binMin, rescotAlpha_binN, rescotAlpha_binWidth, rescotBeta_binMin, rescotBeta_binN, rescotBeta_binWidth, resqbin_binMin, resqbin_binN, resqbin_binWidth, thePixelResolutionFile1, theXHistos, and theYHistos.

Referenced by SiPixelGaussianSmearingRecHitConverterAlgorithm().

{
  //Hard coded at the moment, can easily be changed to be configurable
  rescotAlpha_binMin = 0.1;
  rescotAlpha_binWidth = 0.1 ;
  rescotAlpha_binN = 4;
  rescotBeta_binMin = 0.;
  rescotBeta_binWidth = 0.15;
  rescotBeta_binN = 4;
  resqbin_binMin = 0;
  resqbin_binWidth = 1;
  resqbin_binN = 4;

  // Initialize the forward histos once and for all, and prepare the random generation
  for ( unsigned cotalphaHistBin=1; cotalphaHistBin<=rescotAlpha_binN; ++cotalphaHistBin )
     for ( unsigned cotbetaHistBin=1; cotbetaHistBin<=rescotBeta_binN; ++cotbetaHistBin )  
         for( unsigned qbinBin=1;  qbinBin<=resqbin_binN; ++qbinBin )  {
        unsigned int edgePixelHistN = cotalphaHistBin * 1000 +  cotbetaHistBin * 10 +  qbinBin;
        theXHistos[edgePixelHistN] = new SimpleHistogramGenerator(
        (TH1F*) thePixelResolutionFile1->Get(  Form( "DQMData/clustFPIX/fhx0%u" ,edgePixelHistN ) ) );
        theYHistos[edgePixelHistN] = new SimpleHistogramGenerator(
        (TH1F*) thePixelResolutionFile1->Get(  Form( "DQMData/clustFPIX/fhy0%u" ,edgePixelHistN ) ) );
        unsigned int PixelHistN = 100000 + 10000 + cotalphaHistBin * 1000 +  cotbetaHistBin * 10 +  qbinBin;
        theXHistos[PixelHistN] = new SimpleHistogramGenerator(
        (TH1F*) thePixelResolutionFile1->Get(  Form( "DQMData/clustFPIX/fhx%u" ,PixelHistN ) ) );
        theYHistos[PixelHistN] = new SimpleHistogramGenerator(
        (TH1F*) thePixelResolutionFile1->Get(  Form( "DQMData/clustFPIX/fhy%u" ,PixelHistN ) ) );
     }//end cotalphaHistBin, cotbetaHistBin, qbinBin

  for ( unsigned cotalphaHistBin=1; cotalphaHistBin<=rescotAlpha_binN; ++cotalphaHistBin )
    for ( unsigned cotbetaHistBin=1; cotbetaHistBin<=rescotBeta_binN; ++cotbetaHistBin )
    {
      unsigned int SingleBigPixelHistN = 100000 + cotalphaHistBin * 100 + cotbetaHistBin;
      theXHistos[SingleBigPixelHistN] = new SimpleHistogramGenerator(
      (TH1F*) thePixelResolutionFile1->Get(  Form( "DQMData/clustFPIX/fhx%u" ,SingleBigPixelHistN ) ) );
      theYHistos[SingleBigPixelHistN] = new SimpleHistogramGenerator(
      (TH1F*) thePixelResolutionFile1->Get(  Form( "DQMData/clustFPIX/fhy%u" ,SingleBigPixelHistN ) ) );
      unsigned int SinglePixelHistN = 100000 + 1000 + cotalphaHistBin * 100 + cotbetaHistBin;
      theXHistos[SinglePixelHistN]  = new SimpleHistogramGenerator(
      (TH1F*) thePixelResolutionFile1->Get(  Form( "DQMData/clustFPIX/fhx%u" ,SinglePixelHistN ) ) );
      theYHistos[SinglePixelHistN]  = new SimpleHistogramGenerator(
      (TH1F*) thePixelResolutionFile1->Get(  Form( "DQMData/clustFPIX/fhy%u" ,SinglePixelHistN ) ) );

    }
}

bool SiPixelGaussianSmearingRecHitConverterAlgorithm::isFlipped ( const PixelGeomDetUnit *  theDet ) const

private

Definition at line 489 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.cc.

References PV3DBase< T, PVType, FrameType >::perp(), GeomDet::surface(), and Surface::toGlobal().

Referenced by smearHit().

                                                                                                    {
  // Check the relative position of the local +/- z in global coordinates.
  float tmp1 = theDet->surface().toGlobal(Local3DPoint(0.,0.,0.)).perp();
  float tmp2 = theDet->surface().toGlobal(Local3DPoint(0.,0.,1.)).perp();
  //  std::cout << " 1: " << tmp1 << " 2: " << tmp2 << std::endl;
  if ( tmp2<tmp1 ) return true;
  else return false;    
}

void SiPixelGaussianSmearingRecHitConverterAlgorithm::smearHit	(	const PSimHit &	simHit,
		const PixelGeomDetUnit *	detUnit,
		const double	boundX,
		const double	boundY,
		RandomEngineAndDistribution const *	random
	)

Definition at line 115 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.cc.

References BHX, BHY, BXSIZE, BYSIZE, RectangularPixelTopology::containsBigPixelInX(), RectangularPixelTopology::containsBigPixelInY(), counter, gather_cfg::cout, alignCSCRings::e, prof2calltree::edge, RandomEngineAndDistribution::flatShoot(), i, SiPixelTemplate::interpolate(), isFlipped(), isForward, RectangularPixelTopology::isItBigPixelInX(), RectangularPixelTopology::isItBigPixelInY(), RectangularPixelTopology::isItEdgePixelInX(), RectangularPixelTopology::isItEdgePixelInY(), PSimHit::localDirection(), PSimHit::localPosition(), RectangularPixelTopology::localPosition(), RectangularPixelTopology::measurementPosition(), microntocm, PSimHit::momentumAtEntry(), PixelTopology::ncolumns(), PixelTopology::nrows(), SiPixelTemplate::qbin_dist(), rescotAlpha_binMin, rescotAlpha_binN, rescotAlpha_binWidth, rescotBeta_binMin, rescotBeta_binN, rescotBeta_binWidth, SiPixelTemplate::s50(), jetcorrextractor::sign(), PixelGeomDetType::specificTopology(), PixelGeomDetUnit::specificType(), SiPixelTemplate::temperrors(), tempId, theClslenx, theClsleny, theError, theErrorX, theErrorY, theErrorZ, thePixelPart, thePixelTemp_, thePosition, thePositionX, thePositionY, thePositionZ, theXHistos, theYHistos, Vector3DBase< T, FrameTag >::unit(), PV2DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::x(), xsize, SiPixelTemplate::xtemp(), PV2DBase< T, PVType, FrameType >::y(), PV3DBase< T, PVType, FrameType >::y(), ysize, SiPixelTemplate::ytemp(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by SiTrackerGaussianSmearingRecHitConverter::gaussianSmearing().

{

#ifdef FAMOS_DEBUG
  std::cout << " Pixel smearing in " << thePixelPart 
        << std::endl;
#endif
  //
  // at the beginning the position is the Local Point in the local pixel module reference frame
  // same code as in PixelCPEBase
  LocalVector localDir = simHit.momentumAtEntry().unit();
  float locx = localDir.x();
  float locy = localDir.y();
  float locz = localDir.z();

  // alpha: angle with respect to local x axis in local (x,z) plane
  float cotalpha = locx/locz;
  if ( isFlipped( detUnit ) ) { // &&& check for FPIX !!!
#ifdef FAMOS_DEBUG
    std::cout << " isFlipped " << std::endl;
#endif
  }
  // beta: angle with respect to local y axis in local (y,z) plane
  float cotbeta = locy/locz;
  float sign=1.;
  if( isForward ) {
    if( cotbeta < 0 ) sign=-1.;
    cotbeta = sign*cotbeta;
  }
  
  
  //
#ifdef FAMOS_DEBUG
  std::cout << " Local Direction " << simHit.localDirection()
        << " cotalpha(x) = " << cotalpha
        << " cotbeta(y) = "  << cotbeta
        << std::endl;
#endif
  
  const PixelTopology* theSpecificTopology = &(detUnit->specificType().specificTopology());
  const RectangularPixelTopology *rectPixelTopology = static_cast<const RectangularPixelTopology*>(theSpecificTopology);

  const int nrows = theSpecificTopology->nrows();
  const int ncolumns = theSpecificTopology->ncolumns();

  const Local3DPoint lp = simHit.localPosition();
  //Transform local position to measurement position
  const MeasurementPoint mp = rectPixelTopology->measurementPosition( lp );
  float mpy = mp.y();
  float mpx = mp.x();
  //Get the center of the struck pixel in measurement position
  float pixelCenterY = 0.5 + (int)mpy;
  float pixelCenterX = 0.5 + (int)mpx;
#ifdef FAMOS_DEBUG
  cout<<"Struck pixel center at pitch units x: "<<pixelCenterX<<" y: "<<pixelCenterY<<endl;
#endif

  const MeasurementPoint mpCenter(pixelCenterX, pixelCenterY);
  //Transform the center of the struck pixel back into local position
  const Local3DPoint lpCenter = rectPixelTopology->localPosition( mpCenter );
#ifdef FAMOS_DEBUG
  cout<<"Struck point at cm x: "<<lp.x()<<" y: "<<lp.y()<<endl;
  cout<<"Struck pixel center at cm x: "<<lpCenter.x()<<" y: "<<lpCenter.y()<<endl;
  cout<<"The boundX is "<<boundX<<" boundY is "<<boundY<<endl;
#endif

  //Get the relative position of struck point to the center of the struck pixel
  float xtrk = lp.x() - lpCenter.x();
  float ytrk = lp.y() - lpCenter.y();
  //Pixel Y, X pitch
  const float ysize={0.015}, xsize={0.01};
  //Variables for SiPixelTemplate input, see SiPixelTemplate reco
  float yhit = 20. + 8.*(ytrk/ysize);
  float xhit = 20. + 8.*(xtrk/xsize);
  int   ybin = (int)yhit;
  int   xbin = (int)xhit;
  float yfrac= yhit - (float)ybin;
  float xfrac= xhit - (float)xbin;
  //Protect againt ybin, xbin being outside of range [0-39]
  if( ybin < 0 )    ybin = 0;
  if( ybin > 39 )   ybin = 39;
  if( xbin < 0 )    xbin = 0;
  if( xbin > 39 )   xbin = 39; 

  //Variables for SiPixelTemplate output
  //qBin -- normalized pixel charge deposition
  float qbin_frac[4];
  //Single pixel cluster projection possibility
  float ny1_frac, ny2_frac, nx1_frac, nx2_frac;
  bool singlex = false, singley = false;
  SiPixelTemplate templ(thePixelTemp_);
  templ.interpolate(tempId, cotalpha, cotbeta);
  templ.qbin_dist(tempId, cotalpha, cotbeta, qbin_frac, ny1_frac, ny2_frac, nx1_frac, nx2_frac );
  int  nqbin;

  double xsizeProbability = random->flatShoot();
  double ysizeProbability = random->flatShoot();
  bool hitbigx = rectPixelTopology->isItBigPixelInX( (int)mpx );
  bool hitbigy = rectPixelTopology->isItBigPixelInY( (int)mpy );
  
  if( hitbigx ) 
    if( xsizeProbability < nx2_frac )  singlex = true;
    else singlex = false;
  else
    if( xsizeProbability < nx1_frac )  singlex = true;
    else singlex = false;

  if( hitbigy )
    if( ysizeProbability < ny2_frac )  singley = true;
    else singley = false;
  else
    if( ysizeProbability < ny1_frac )  singley = true;
    else singley = false;
  


  // random multiplicity for alpha and beta
  double qbinProbability = random->flatShoot();
  for(int i = 0; i<4; ++i) {
     nqbin = i;
     if(qbinProbability < qbin_frac[i]) break;
  }

  //Store interpolated pixel cluster profile
  //BYSIZE, BXSIZE, const definition from SiPixelTemplate
  float ytempl[41][BYSIZE] = {{0}}, xtempl[41][BXSIZE] = {{0}} ;
  templ.ytemp(0, 40, ytempl);
  templ.xtemp(0, 40, xtempl);

  std::vector<double> ytemp(BYSIZE);
  for( int i=0; i<BYSIZE; ++i) {
     ytemp[i]=(1.-yfrac)*ytempl[ybin][i]+yfrac*ytempl[ybin+1][i];
  }

  std::vector<double> xtemp(BXSIZE);
  for(int i=0; i<BXSIZE; ++i) {
     xtemp[i]=(1.-xfrac)*xtempl[xbin][i]+xfrac*xtempl[xbin+1][i];
  }

  //Pixel readout threshold
  const float qThreshold = templ.s50()*2.0;

  //Cut away pixels below readout threshold
  //For cluster lengths calculation
  int offsetX1=0, offsetX2=0, offsetY1=0, offsetY2=0;
  int firstY, lastY, firstX, lastX;
  for( firstY = 0; firstY < BYSIZE; ++firstY ) {
    bool yCluster = ytemp[firstY] > qThreshold ;
    if( yCluster )
    {
      offsetY1 = BHY -firstY;
      break;
    }
  }
  for( lastY = firstY; lastY < BYSIZE; ++lastY )
  {
    bool yCluster = ytemp[lastY] > qThreshold ;
    if( !yCluster )
    {
      lastY = lastY - 1;
      offsetY2 = lastY - BHY;
      break;
    }
  }

  for( firstX = 0; firstX < BXSIZE; ++firstX )  {
    bool xCluster = xtemp[firstX] > qThreshold ;
    if( xCluster ) {
      offsetX1 = BHX - firstX;
      break;
    }
  }
  for( lastX = firstX; lastX < BXSIZE; ++ lastX ) {
    bool xCluster = xtemp[lastX] > qThreshold ;
    if( !xCluster ) {
      lastX = lastX - 1;
      offsetX2 = lastX - BHX;
      break;
    }
  }

  bool edge, edgex, edgey;
  //  bool bigx, bigy;
  unsigned int clslenx = offsetX1 + offsetX2 + 1;
  unsigned int clsleny = offsetY1 + offsetY2 + 1;

  theClslenx = clslenx;
  theClsleny = clsleny;

  int firstPixelInX = (int)mpx - offsetX1 ;
  int firstPixelInY = (int)mpy - offsetY1 ;
  int lastPixelInX  = (int)mpx + offsetX2 ;
  int lastPixelInY  = (int)mpy + offsetY2 ;
  firstPixelInX = (firstPixelInX >= 0) ? firstPixelInX : 0 ;
  firstPixelInY = (firstPixelInY >= 0) ? firstPixelInY : 0 ;
  lastPixelInX  = (lastPixelInX < nrows ) ? lastPixelInX : nrows-1 ;
  lastPixelInY  = (lastPixelInY < ncolumns ) ? lastPixelInY : ncolumns-1;

  edgex = rectPixelTopology->isItEdgePixelInX( firstPixelInX ) || rectPixelTopology->isItEdgePixelInX( lastPixelInX );
  edgey = rectPixelTopology->isItEdgePixelInY( firstPixelInY ) || rectPixelTopology->isItEdgePixelInY( lastPixelInY );
  edge = edgex || edgey;

  //  bigx = rectPixelTopology->isItBigPixelInX( firstPixelInX ) || rectPixelTopology->isItBigPixelInX( lastPixelInX );
  //  bigy = rectPixelTopology->isItBigPixelInY( firstPixelInY ) || rectPixelTopology->isItBigPixelInY( lastPixelInY );
  bool hasBigPixelInX = rectPixelTopology->containsBigPixelInX( firstPixelInX, lastPixelInX );
  bool hasBigPixelInY = rectPixelTopology->containsBigPixelInY( firstPixelInY, lastPixelInY );

  //Variables for SiPixelTemplate pixel hit error output
  float sigmay, sigmax, sy1, sy2, sx1, sx2;
  templ.temperrors(tempId, cotalpha, cotbeta, nqbin, sigmay, sigmax, sy1, sy2, sx1, sx2 );

  // define private mebers --> Errors
  if( edge ) {
    if( edgex && !edgey ) {
      theErrorX = 23.0*microntocm;
      theErrorY = 39.0*microntocm;
    }
    else if( !edgex && edgey ) {
      theErrorX = 24.0*microntocm;
      theErrorY = 96.0*microntocm;
    }
    else
    {
      theErrorX = 31.0*microntocm;
      theErrorY = 90.0*microntocm;
    }
    
  }  
  else {
    if( singlex )
      if ( hitbigx )
        theErrorX = sx2*microntocm;
      else
        theErrorX = sx1*microntocm;
    else  theErrorX = sigmax*microntocm;
    if( singley )
      if( hitbigy )
        theErrorY = sy2*microntocm;
      else
        theErrorY = sy1*microntocm;
    else  theErrorY = sigmay*microntocm;
  }
  theErrorZ = 1e-8; // 1 um means zero
  theError = LocalError( theErrorX*theErrorX, 0., theErrorY*theErrorY);
  // Local Error is 2D: (xx,xy,yy), square of sigma in first an third position 
  // as for resolution matrix
  //
#ifdef FAMOS_DEBUG
  std::cout << " Pixel Errors "
        << "\talpha(x) = " << theErrorX
        << "\tbeta(y) = "  << theErrorY
        << std::endl;   
#endif
  // Generate position
  // get resolution histograms
  int cotalphaHistBin = (int)( ( cotalpha - rescotAlpha_binMin ) / rescotAlpha_binWidth + 1 );
  int cotbetaHistBin  = (int)( ( cotbeta  - rescotBeta_binMin )  / rescotBeta_binWidth + 1 );
  // protection against out-of-range (undeflows and overflows)
  if( cotalphaHistBin < 1 ) cotalphaHistBin = 1; 
  if( cotbetaHistBin  < 1 ) cotbetaHistBin  = 1; 
  if( cotalphaHistBin > (int)rescotAlpha_binN ) cotalphaHistBin = (int)rescotAlpha_binN; 
  if( cotbetaHistBin  > (int)rescotBeta_binN  ) cotbetaHistBin  = (int)rescotBeta_binN; 
  //
  unsigned int theXHistN;
  unsigned int theYHistN;
  if( !isForward ) {
     if(edge)
     {
       theXHistN = cotalphaHistBin * 1000 + cotbetaHistBin * 10  +  (nqbin+1);
       theYHistN = theXHistN;         
     }
     else
     {
       if(singlex)
       {
     if(hitbigx)  theXHistN = 1 * 100000 + cotalphaHistBin * 100 + cotbetaHistBin ;
     else   theXHistN = 1 * 100000 + 1 * 1000 + cotalphaHistBin * 100 + cotbetaHistBin ;
       }
       else
       {
     if(hasBigPixelInX)  theXHistN = 1 * 1000000 + 1 * 100000 + cotalphaHistBin * 1000 + cotbetaHistBin * 10 + (nqbin+1);
     else   theXHistN = 1 * 1000000 + 1 * 100000 + 1 * 10000 + cotalphaHistBin * 1000 + cotbetaHistBin * 10 + (nqbin+1);
       }
       if(singley)
       {
     if(hitbigy)  theYHistN = 1 * 100000 + cotalphaHistBin * 100 + cotbetaHistBin ;
     else   theYHistN = 1 * 100000 + 1 * 1000 + cotalphaHistBin * 100 + cotbetaHistBin ;
       }
       else
       {
     if(hasBigPixelInY)  theYHistN = 1 * 1000000 + 1 * 100000 + cotalphaHistBin * 1000 + cotbetaHistBin * 10 + (nqbin+1);
     else   theYHistN = 1 * 1000000 + 1 * 100000 + 1 * 10000 + cotalphaHistBin * 1000 + cotbetaHistBin * 10 + (nqbin+1);
       }
     }
  }
  else
  {
     if(edge)
     {
       theXHistN = cotalphaHistBin * 1000 +  cotbetaHistBin * 10 +  (nqbin+1);
       theYHistN = theXHistN;
     }
     else
     {
       if( singlex )
         if( hitbigx )
       theXHistN = 100000 + cotalphaHistBin * 100 + cotbetaHistBin;
     else
       theXHistN = 100000 + 1000 + cotalphaHistBin * 100 + cotbetaHistBin;
       else
         theXHistN = 100000 + 10000 + cotalphaHistBin * 1000 +  cotbetaHistBin * 10 +  (nqbin+1);
       if( singley )
         if( hitbigy )
       theYHistN = 100000 + cotalphaHistBin * 100 + cotbetaHistBin;
     else
       theYHistN = 100000 + 1000 + cotalphaHistBin * 100 + cotbetaHistBin;
       else
         theYHistN = 100000 + 10000 + cotalphaHistBin * 1000 +  cotbetaHistBin * 10 +  (nqbin+1);
     }
  }
  unsigned int counter = 0;
  do {
    //
    // Smear the hit Position
    thePositionX = theXHistos[theXHistN]->generate(random);
    thePositionY = theYHistos[theYHistN]->generate(random);
    if( isForward ) thePositionY *= sign;
    thePositionZ = 0.0; // set at the centre of the active area
    //protect from empty resolution histograms
    //if( thePositionX > 0.0799 )  thePositionX = 0;
    //if( thePositionY > 0.0799 )  thePositionY = 0;
    thePosition = 
      Local3DPoint(simHit.localPosition().x() + thePositionX , 
                   simHit.localPosition().y() + thePositionY , 
                   simHit.localPosition().z() + thePositionZ );
#ifdef FAMOS_DEBUG
    std::cout << " Detector bounds: "
              << "\t\tx = " << boundX
              << "\ty = " << boundY
              << std::endl;
    std::cout << " Generated local position "
              << "\tx = " << thePosition.x()
              << "\ty = " << thePosition.y()
              << std::endl;       
#endif  
    counter++;
    if(counter > 20) {
      thePosition = Local3DPoint(simHit.localPosition().x(), 
                                 simHit.localPosition().y(), 
                                 simHit.localPosition().z());
      break;
    }
  } while(fabs(thePosition.x()) > boundX  || fabs(thePosition.y()) > boundY);

}

Member Data Documentation

bool SiPixelGaussianSmearingRecHitConverterAlgorithm::isForward

private

Definition at line 74 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

Referenced by SiPixelGaussianSmearingRecHitConverterAlgorithm(), and smearHit().

edm::ParameterSet SiPixelGaussianSmearingRecHitConverterAlgorithm::pset_

private

Definition at line 85 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

Referenced by SiPixelGaussianSmearingRecHitConverterAlgorithm().

double SiPixelGaussianSmearingRecHitConverterAlgorithm::rescotAlpha_binMin

private

Definition at line 78 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

Referenced by initializeBarrel(), initializeForward(), and smearHit().

unsigned int SiPixelGaussianSmearingRecHitConverterAlgorithm::rescotAlpha_binN

private

Definition at line 79 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

Referenced by initializeBarrel(), initializeForward(), and smearHit().

double SiPixelGaussianSmearingRecHitConverterAlgorithm::rescotAlpha_binWidth

private

Definition at line 78 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

Referenced by initializeBarrel(), initializeForward(), and smearHit().

double SiPixelGaussianSmearingRecHitConverterAlgorithm::rescotBeta_binMin

private

Definition at line 80 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

Referenced by initializeBarrel(), initializeForward(), and smearHit().

unsigned int SiPixelGaussianSmearingRecHitConverterAlgorithm::rescotBeta_binN

private

Definition at line 81 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

Referenced by initializeBarrel(), initializeForward(), and smearHit().

double SiPixelGaussianSmearingRecHitConverterAlgorithm::rescotBeta_binWidth

private

Definition at line 80 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

Referenced by initializeBarrel(), initializeForward(), and smearHit().

int SiPixelGaussianSmearingRecHitConverterAlgorithm::resqbin_binMin

private

Definition at line 82 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

Referenced by initializeBarrel(), and initializeForward().

unsigned int SiPixelGaussianSmearingRecHitConverterAlgorithm::resqbin_binN

private

Definition at line 83 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

Referenced by initializeBarrel(), and initializeForward().

int SiPixelGaussianSmearingRecHitConverterAlgorithm::resqbin_binWidth

private

Definition at line 82 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

Referenced by initializeBarrel(), and initializeForward().

int SiPixelGaussianSmearingRecHitConverterAlgorithm::tempId

private

Definition at line 68 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

Referenced by SiPixelGaussianSmearingRecHitConverterAlgorithm(), and smearHit().

unsigned int SiPixelGaussianSmearingRecHitConverterAlgorithm::theClslenx

private

Definition at line 108 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

Referenced by getPixelMultiplicityAlpha(), and smearHit().

unsigned int SiPixelGaussianSmearingRecHitConverterAlgorithm::theClsleny

private

Definition at line 109 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

Referenced by getPixelMultiplicityBeta(), and smearHit().

LocalError SiPixelGaussianSmearingRecHitConverterAlgorithm::theError

private

Definition at line 104 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

Referenced by getError(), and smearHit().

double SiPixelGaussianSmearingRecHitConverterAlgorithm::theErrorX

private

Definition at line 105 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

Referenced by getErrorX(), and smearHit().

double SiPixelGaussianSmearingRecHitConverterAlgorithm::theErrorY

private

Definition at line 106 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

Referenced by getErrorY(), and smearHit().

double SiPixelGaussianSmearingRecHitConverterAlgorithm::theErrorZ

private

Definition at line 107 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

Referenced by getErrorZ(), and smearHit().

unsigned int SiPixelGaussianSmearingRecHitConverterAlgorithm::theLayer

private

Definition at line 98 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

GeomDetType::SubDetector SiPixelGaussianSmearingRecHitConverterAlgorithm::thePixelPart

private

Definition at line 87 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

Referenced by SiPixelGaussianSmearingRecHitConverterAlgorithm(), and smearHit().

TFile* SiPixelGaussianSmearingRecHitConverterAlgorithm::thePixelResolutionFile1

private

Definition at line 92 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

Referenced by initializeBarrel(), initializeForward(), and SiPixelGaussianSmearingRecHitConverterAlgorithm().

TFile* SiPixelGaussianSmearingRecHitConverterAlgorithm::thePixelResolutionFile2

private

Definition at line 95 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

Referenced by initializeBarrel(), and SiPixelGaussianSmearingRecHitConverterAlgorithm().

std::string SiPixelGaussianSmearingRecHitConverterAlgorithm::thePixelResolutionFileName1

private

Definition at line 93 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

Referenced by SiPixelGaussianSmearingRecHitConverterAlgorithm().

std::string SiPixelGaussianSmearingRecHitConverterAlgorithm::thePixelResolutionFileName2

private

Definition at line 96 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

Referenced by SiPixelGaussianSmearingRecHitConverterAlgorithm().

std::vector< SiPixelTemplateStore > SiPixelGaussianSmearingRecHitConverterAlgorithm::thePixelTemp_

private

Definition at line 67 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

Referenced by SiPixelGaussianSmearingRecHitConverterAlgorithm(), and smearHit().

Local3DPoint SiPixelGaussianSmearingRecHitConverterAlgorithm::thePosition

private

Definition at line 100 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

Referenced by getPosition(), and smearHit().

double SiPixelGaussianSmearingRecHitConverterAlgorithm::thePositionX

private

Definition at line 101 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

Referenced by getPositionX(), and smearHit().

double SiPixelGaussianSmearingRecHitConverterAlgorithm::thePositionY

private

Definition at line 102 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

Referenced by getPositionY(), and smearHit().

double SiPixelGaussianSmearingRecHitConverterAlgorithm::thePositionZ

private

Definition at line 103 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

Referenced by getPositionZ(), and smearHit().

std::map<unsigned,const SimpleHistogramGenerator*> SiPixelGaussianSmearingRecHitConverterAlgorithm::theXHistos

private

Definition at line 89 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

Referenced by initializeBarrel(), initializeForward(), SiPixelGaussianSmearingRecHitConverterAlgorithm(), smearHit(), and ~SiPixelGaussianSmearingRecHitConverterAlgorithm().

std::map<unsigned,const SimpleHistogramGenerator*> SiPixelGaussianSmearingRecHitConverterAlgorithm::theYHistos

private

Definition at line 90 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

Referenced by initializeBarrel(), initializeForward(), SiPixelGaussianSmearingRecHitConverterAlgorithm(), smearHit(), and ~SiPixelGaussianSmearingRecHitConverterAlgorithm().

bool SiPixelGaussianSmearingRecHitConverterAlgorithm::useCMSSWPixelParameterization

private

Definition at line 65 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.h.

Referenced by SiPixelGaussianSmearingRecHitConverterAlgorithm().