CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
List of all members | Public Member Functions | Private Member Functions | Private Attributes
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, const RandomEngine *engine)
 
void smearHit (const PSimHit &simHit, const PixelGeomDetUnit *detUnit, const double boundX, const double boundY)
 
virtual ~SiPixelGaussianSmearingRecHitConverterAlgorithm ()
 

Private Member Functions

void initializeBarrel ()
 
void initializeForward ()
 
bool isFlipped (const PixelGeomDetUnit *theDet) const
 

Private Attributes

bool isForward
 
edm::ParameterSet pset_
 
const RandomEnginerandom
 
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
 
SiPixelTemplate templ
 
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
 
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,
const RandomEngine engine 
)
explicit

Definition at line 41 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.cc.

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

45 :
46  pset_(pset),
47  thePixelPart(pixelPart),
48  random(engine)
49 {
50  // Switch between old (ORCA) and new (CMSSW) pixel parameterization
51  useCMSSWPixelParameterization = pset.getParameter<bool>("UseCMSSWPixelParametrization");
52 
54  isForward = false;
55  thePixelResolutionFileName1 = pset_.getParameter<string>( "NewPixelBarrelResolutionFile1" );
57  TFile( edm::FileInPath( thePixelResolutionFileName1 ).fullPath().c_str() ,"READ");
58  thePixelResolutionFileName2 = pset_.getParameter<string>( "NewPixelBarrelResolutionFile2" );
60  TFile( edm::FileInPath( thePixelResolutionFileName2 ).fullPath().c_str() ,"READ");
62  tempId = pset_.getParameter<int> ( "templateIdBarrel" );
63  if( ! templ.pushfile(tempId) )
64  throw cms::Exception("SiPixelGaussianSmearingRecHitConverterAlgorithm:")
65  <<"SiPixel Barrel Template Not Loaded Correctly!"<<endl;
66 #ifdef FAMOS_DEBUG
67  cout<<"The Barrel map size is "<<theXHistos.size()<<" and "<<theYHistos.size()<<endl;
68 #endif
69  }
70  else
72  isForward = true;
73  thePixelResolutionFileName1 = pset_.getParameter<string>( "NewPixelForwardResolutionFile" );
74  thePixelResolutionFile1 = new
75  TFile( edm::FileInPath( thePixelResolutionFileName1 ).fullPath().c_str() ,"READ");
77  tempId = pset_.getParameter<int> ( "templateIdForward" );
78  if( ! templ.pushfile(tempId) )
79  throw cms::Exception("SiPixelGaussianSmearingRecHitConverterAlgorithm:")
80  <<"SiPixel Forward Template Not Loaded Correctly!"<<endl;
81 #ifdef FAMOS_DEBUG
82  cout<<"The Forward map size is "<<theXHistos.size()<<" and "<<theYHistos.size()<<endl;
83 #endif
84  }
85  else
86  throw cms::Exception("SiPixelGaussianSmearingRecHitConverterAlgorithm :")
87  <<"Not a pixel detector"<<endl;
88 }
T getParameter(std::string const &) const
std::map< unsigned, const SimpleHistogramGenerator * > theXHistos
bool pushfile(int filenum)
std::map< unsigned, const SimpleHistogramGenerator * > theYHistos
tuple cout
Definition: gather_cfg.py:121
SiPixelGaussianSmearingRecHitConverterAlgorithm::~SiPixelGaussianSmearingRecHitConverterAlgorithm ( )
virtual

Definition at line 90 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.cc.

References theXHistos, and theYHistos.

91 {
92 
93  std::map<unsigned,const SimpleHistogramGenerator*>::const_iterator it;
94  for ( it=theXHistos.begin(); it!=theXHistos.end(); ++it )
95  delete it->second;
96  for ( it=theYHistos.begin(); it!=theYHistos.end(); ++it )
97  delete it->second;
98 
99  theXHistos.clear();
100  theYHistos.clear();
101 
102 }
std::map< unsigned, const SimpleHistogramGenerator * > theXHistos
std::map< unsigned, const SimpleHistogramGenerator * > theYHistos

Member Function Documentation

LocalError SiPixelGaussianSmearingRecHitConverterAlgorithm::getError ( )
inline
double SiPixelGaussianSmearingRecHitConverterAlgorithm::getErrorX ( )
inline
double SiPixelGaussianSmearingRecHitConverterAlgorithm::getErrorY ( )
inline
double SiPixelGaussianSmearingRecHitConverterAlgorithm::getErrorZ ( )
inline
unsigned int SiPixelGaussianSmearingRecHitConverterAlgorithm::getPixelMultiplicityAlpha ( )
inline
unsigned int SiPixelGaussianSmearingRecHitConverterAlgorithm::getPixelMultiplicityBeta ( )
inline
Local3DPoint SiPixelGaussianSmearingRecHitConverterAlgorithm::getPosition ( )
inline
double SiPixelGaussianSmearingRecHitConverterAlgorithm::getPositionX ( )
inline
double SiPixelGaussianSmearingRecHitConverterAlgorithm::getPositionY ( )
inline
double SiPixelGaussianSmearingRecHitConverterAlgorithm::getPositionZ ( )
inline
void SiPixelGaussianSmearingRecHitConverterAlgorithm::initializeBarrel ( )
private

Definition at line 485 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.cc.

References random, 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().

486 {
487  //Hard coded at the moment, can easily be changed to be configurable
488  rescotAlpha_binMin = -0.2;
489  rescotAlpha_binWidth = 0.08 ;
490  rescotAlpha_binN = 5;
491  rescotBeta_binMin = -5.5;
492  rescotBeta_binWidth = 1.0;
493  rescotBeta_binN = 11;
494  resqbin_binMin = 0;
495  resqbin_binWidth = 1;
496  resqbin_binN = 4;
497 
498  // Initialize the barrel histos once and for all, and prepare the random generation
499  for ( unsigned cotalphaHistBin=1; cotalphaHistBin<=rescotAlpha_binN; ++cotalphaHistBin )
500  for ( unsigned cotbetaHistBin=1; cotbetaHistBin<=rescotBeta_binN; ++cotbetaHistBin ) {
501  unsigned int singleBigPixelHistN = 1 * 100000
502  + cotalphaHistBin * 100
503  + cotbetaHistBin ;
504  theXHistos[singleBigPixelHistN] = new SimpleHistogramGenerator(
505  (TH1F*) thePixelResolutionFile1->Get( Form( "DQMData/clustBPIX/hx%u" , singleBigPixelHistN ) ),
506  random);
507  theYHistos[singleBigPixelHistN] = new SimpleHistogramGenerator(
508  (TH1F*) thePixelResolutionFile1->Get( Form( "DQMData/clustBPIX/hy%u" , singleBigPixelHistN ) ),
509  random);
510  unsigned int singlePixelHistN = 1 * 100000 + 1 * 1000
511  + cotalphaHistBin * 100
512  + cotbetaHistBin ;
513  theXHistos[singlePixelHistN] = new SimpleHistogramGenerator(
514  (TH1F*) thePixelResolutionFile1->Get( Form( "DQMData/clustBPIX/hx%u" , singlePixelHistN ) ),
515  random);
516  theYHistos[singlePixelHistN] = new SimpleHistogramGenerator(
517  (TH1F*) thePixelResolutionFile1->Get( Form( "DQMData/clustBPIX/hy%u" , singlePixelHistN ) ),
518  random);
519  for( unsigned qbinBin=1; qbinBin<=resqbin_binN; ++qbinBin ) {
520  unsigned int edgePixelHistN = cotalphaHistBin * 1000
521  + cotbetaHistBin * 10
522  + qbinBin;
523  theXHistos[edgePixelHistN] = new SimpleHistogramGenerator(
524  (TH1F*) thePixelResolutionFile2->Get( Form( "DQMData/clustBPIX/hx0%u" ,edgePixelHistN ) ), random );
525  theYHistos[edgePixelHistN] = new SimpleHistogramGenerator(
526  (TH1F*) thePixelResolutionFile2->Get( Form( "DQMData/clustBPIX/hy0%u" ,edgePixelHistN ) ), random );
527  unsigned int multiPixelBigHistN = 1 * 1000000 + 1 * 100000
528  + cotalphaHistBin * 1000
529  + cotbetaHistBin * 10
530  + qbinBin;
531  theXHistos[multiPixelBigHistN] = new SimpleHistogramGenerator(
532  (TH1F*) thePixelResolutionFile1->Get( Form( "DQMData/clustBPIX/hx%u" ,multiPixelBigHistN ) ), random );
533  theYHistos[multiPixelBigHistN] = new SimpleHistogramGenerator(
534  (TH1F*) thePixelResolutionFile1->Get( Form( "DQMData/clustBPIX/hy%u" ,multiPixelBigHistN ) ), random );
535  unsigned int multiPixelHistN = 1 * 1000000 + 1 * 100000 + 1 * 10000
536  + cotalphaHistBin * 1000
537  + cotbetaHistBin * 10
538  + qbinBin;
539  theXHistos[multiPixelHistN] = new SimpleHistogramGenerator(
540  (TH1F*) thePixelResolutionFile1->Get( Form( "DQMData/clustBPIX/hx%u" , multiPixelHistN ) ),
541  random);
542  theYHistos[multiPixelHistN] = new SimpleHistogramGenerator(
543  (TH1F*) thePixelResolutionFile1->Get( Form( "DQMData/clustBPIX/hy%u" , multiPixelHistN ) ),
544  random);
545  } //end for qbinBin
546  }//end for cotalphaHistBin, cotbetaHistBin
547 }
std::map< unsigned, const SimpleHistogramGenerator * > theXHistos
std::map< unsigned, const SimpleHistogramGenerator * > theYHistos
void SiPixelGaussianSmearingRecHitConverterAlgorithm::initializeForward ( )
private

Definition at line 549 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.cc.

References random, 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().

550 {
551  //Hard coded at the moment, can easily be changed to be configurable
552  rescotAlpha_binMin = 0.1;
553  rescotAlpha_binWidth = 0.1 ;
554  rescotAlpha_binN = 4;
555  rescotBeta_binMin = 0.;
556  rescotBeta_binWidth = 0.15;
557  rescotBeta_binN = 4;
558  resqbin_binMin = 0;
559  resqbin_binWidth = 1;
560  resqbin_binN = 4;
561 
562  // Initialize the forward histos once and for all, and prepare the random generation
563  for ( unsigned cotalphaHistBin=1; cotalphaHistBin<=rescotAlpha_binN; ++cotalphaHistBin )
564  for ( unsigned cotbetaHistBin=1; cotbetaHistBin<=rescotBeta_binN; ++cotbetaHistBin )
565  for( unsigned qbinBin=1; qbinBin<=resqbin_binN; ++qbinBin ) {
566  unsigned int edgePixelHistN = cotalphaHistBin * 1000 + cotbetaHistBin * 10 + qbinBin;
567  theXHistos[edgePixelHistN] = new SimpleHistogramGenerator(
568  (TH1F*) thePixelResolutionFile1->Get( Form( "DQMData/clustFPIX/fhx0%u" ,edgePixelHistN ) ), random );
569  theYHistos[edgePixelHistN] = new SimpleHistogramGenerator(
570  (TH1F*) thePixelResolutionFile1->Get( Form( "DQMData/clustFPIX/fhy0%u" ,edgePixelHistN ) ), random );
571  unsigned int PixelHistN = 100000 + 10000 + cotalphaHistBin * 1000 + cotbetaHistBin * 10 + qbinBin;
572  theXHistos[PixelHistN] = new SimpleHistogramGenerator(
573  (TH1F*) thePixelResolutionFile1->Get( Form( "DQMData/clustFPIX/fhx%u" ,PixelHistN ) ), random );
574  theYHistos[PixelHistN] = new SimpleHistogramGenerator(
575  (TH1F*) thePixelResolutionFile1->Get( Form( "DQMData/clustFPIX/fhy%u" ,PixelHistN ) ), random );
576  }//end cotalphaHistBin, cotbetaHistBin, qbinBin
577 
578  for ( unsigned cotalphaHistBin=1; cotalphaHistBin<=rescotAlpha_binN; ++cotalphaHistBin )
579  for ( unsigned cotbetaHistBin=1; cotbetaHistBin<=rescotBeta_binN; ++cotbetaHistBin )
580  {
581  unsigned int SingleBigPixelHistN = 100000 + cotalphaHistBin * 100 + cotbetaHistBin;
582  theXHistos[SingleBigPixelHistN] = new SimpleHistogramGenerator(
583  (TH1F*) thePixelResolutionFile1->Get( Form( "DQMData/clustFPIX/fhx%u" ,SingleBigPixelHistN ) ), random );
584  theYHistos[SingleBigPixelHistN] = new SimpleHistogramGenerator(
585  (TH1F*) thePixelResolutionFile1->Get( Form( "DQMData/clustFPIX/fhy%u" ,SingleBigPixelHistN ) ), random );
586  unsigned int SinglePixelHistN = 100000 + 1000 + cotalphaHistBin * 100 + cotbetaHistBin;
587  theXHistos[SinglePixelHistN] = new SimpleHistogramGenerator(
588  (TH1F*) thePixelResolutionFile1->Get( Form( "DQMData/clustFPIX/fhx%u" ,SinglePixelHistN ) ), random );
589  theYHistos[SinglePixelHistN] = new SimpleHistogramGenerator(
590  (TH1F*) thePixelResolutionFile1->Get( Form( "DQMData/clustFPIX/fhy%u" ,SinglePixelHistN ) ), random );
591 
592  }
593 }
std::map< unsigned, const SimpleHistogramGenerator * > theXHistos
std::map< unsigned, const SimpleHistogramGenerator * > theYHistos
bool SiPixelGaussianSmearingRecHitConverterAlgorithm::isFlipped ( const PixelGeomDetUnit theDet) const
private

Definition at line 476 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.cc.

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

Referenced by smearHit().

476  {
477  // Check the relative position of the local +/- z in global coordinates.
478  float tmp1 = theDet->surface().toGlobal(Local3DPoint(0.,0.,0.)).perp();
479  float tmp2 = theDet->surface().toGlobal(Local3DPoint(0.,0.,1.)).perp();
480  // std::cout << " 1: " << tmp1 << " 2: " << tmp2 << std::endl;
481  if ( tmp2<tmp1 ) return true;
482  else return false;
483 }
GlobalPoint toGlobal(const Point2DBase< Scalar, LocalTag > lp) const
Definition: Surface.h:114
T perp() const
Definition: PV3DBase.h:72
const Plane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:35
Point3DBase< float, LocalTag > Local3DPoint
Definition: LocalPoint.h:9
void SiPixelGaussianSmearingRecHitConverterAlgorithm::smearHit ( const PSimHit simHit,
const PixelGeomDetUnit detUnit,
const double  boundX,
const double  boundY 
)

Definition at line 104 of file SiPixelGaussianSmearingRecHitConverterAlgorithm.cc.

References BHX, BHY, BXSIZE, BYSIZE, gather_cfg::cout, alignCSCRings::e, prof2calltree::edge, RandomEngine::flatShoot(), i, SiPixelTemplate::interpolate(), isFlipped(), isForward, PSimHit::localDirection(), PSimHit::localPosition(), microntocm, PSimHit::momentumAtEntry(), PixelTopology::ncolumns(), PixelTopology::nrows(), SiPixelTemplate::qbin_dist(), random, rescotAlpha_binMin, rescotAlpha_binN, rescotAlpha_binWidth, rescotBeta_binMin, rescotBeta_binN, rescotBeta_binWidth, SiPixelTemplate::s50(), PixelGeomDetUnit::specificType(), SiPixelTemplate::temperrors(), tempId, templ, theClslenx, theClsleny, theError, theErrorX, theErrorY, theErrorZ, thePixelPart, thePosition, thePositionX, thePositionY, thePositionZ, theXHistos, theYHistos, Vector3DBase< T, FrameTag >::unit(), PV2DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::x(), SiPixelTemplate::xtemp(), PV2DBase< T, PVType, FrameType >::y(), PV3DBase< T, PVType, FrameType >::y(), SiPixelTemplate::ytemp(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by SiTrackerGaussianSmearingRecHitConverter::gaussianSmearing().

109 {
110 
111 #ifdef FAMOS_DEBUG
112  std::cout << " Pixel smearing in " << thePixelPart
113  << std::endl;
114 #endif
115  //
116  // at the beginning the position is the Local Point in the local pixel module reference frame
117  // same code as in PixelCPEBase
118  LocalVector localDir = simHit.momentumAtEntry().unit();
119  float locx = localDir.x();
120  float locy = localDir.y();
121  float locz = localDir.z();
122 
123  // alpha: angle with respect to local x axis in local (x,z) plane
124  float cotalpha = locx/locz;
125  if ( isFlipped( detUnit ) ) { // &&& check for FPIX !!!
126 #ifdef FAMOS_DEBUG
127  std::cout << " isFlipped " << std::endl;
128 #endif
129  }
130  // beta: angle with respect to local y axis in local (y,z) plane
131  float cotbeta = locy/locz;
132  float sign=1.;
133  if( isForward ) {
134  if( cotbeta < 0 ) sign=-1.;
135  cotbeta = sign*cotbeta;
136  }
137 
138 
139  //
140 #ifdef FAMOS_DEBUG
141  std::cout << " Local Direction " << simHit.localDirection()
142  << " cotalpha(x) = " << cotalpha
143  << " cotbeta(y) = " << cotbeta
144  << std::endl;
145 #endif
146 
147  const PixelTopology* theSpecificTopology = &(detUnit->specificType().specificTopology());
148  const RectangularPixelTopology *rectPixelTopology = static_cast<const RectangularPixelTopology*>(theSpecificTopology);
149 
150  const int nrows = theSpecificTopology->nrows();
151  const int ncolumns = theSpecificTopology->ncolumns();
152 
153  const Local3DPoint lp = simHit.localPosition();
154  //Transform local position to measurement position
155  const MeasurementPoint mp = rectPixelTopology->measurementPosition( lp );
156  float mpy = mp.y();
157  float mpx = mp.x();
158  //Get the center of the struck pixel in measurement position
159  float pixelCenterY = 0.5 + (int)mpy;
160  float pixelCenterX = 0.5 + (int)mpx;
161 #ifdef FAMOS_DEBUG
162  cout<<"Struck pixel center at pitch units x: "<<pixelCenterX<<" y: "<<pixelCenterY<<endl;
163 #endif
164 
165  const MeasurementPoint mpCenter(pixelCenterX, pixelCenterY);
166  //Transform the center of the struck pixel back into local position
167  const Local3DPoint lpCenter = rectPixelTopology->localPosition( mpCenter );
168 #ifdef FAMOS_DEBUG
169  cout<<"Struck point at cm x: "<<lp.x()<<" y: "<<lp.y()<<endl;
170  cout<<"Struck pixel center at cm x: "<<lpCenter.x()<<" y: "<<lpCenter.y()<<endl;
171  cout<<"The boundX is "<<boundX<<" boundY is "<<boundY<<endl;
172 #endif
173 
174  //Get the relative position of struck point to the center of the struck pixel
175  float xtrk = lp.x() - lpCenter.x();
176  float ytrk = lp.y() - lpCenter.y();
177  //Pixel Y, X pitch
178  const float ysize={0.015}, xsize={0.01};
179  //Variables for SiPixelTemplate input, see SiPixelTemplate reco
180  float yhit = 20. + 8.*(ytrk/ysize);
181  float xhit = 20. + 8.*(xtrk/xsize);
182  int ybin = (int)yhit;
183  int xbin = (int)xhit;
184  float yfrac= yhit - (float)ybin;
185  float xfrac= xhit - (float)xbin;
186  //Protect againt ybin, xbin being outside of range [0-39]
187  if( ybin < 0 ) ybin = 0;
188  if( ybin > 39 ) ybin = 39;
189  if( xbin < 0 ) xbin = 0;
190  if( xbin > 39 ) xbin = 39;
191 
192  //Variables for SiPixelTemplate output
193  //qBin -- normalized pixel charge deposition
194  float qbin_frac[4];
195  //Single pixel cluster projection possibility
196  float ny1_frac, ny2_frac, nx1_frac, nx2_frac;
197  bool singlex = false, singley = false;
198  templ.interpolate(tempId, cotalpha, cotbeta);
199  templ.qbin_dist(tempId, cotalpha, cotbeta, qbin_frac, ny1_frac, ny2_frac, nx1_frac, nx2_frac );
200  int nqbin;
201 
202  double xsizeProbability = random->flatShoot();
203  double ysizeProbability = random->flatShoot();
204  bool hitbigx = rectPixelTopology->isItBigPixelInX( (int)mpx );
205  bool hitbigy = rectPixelTopology->isItBigPixelInY( (int)mpy );
206 
207  if( hitbigx )
208  if( xsizeProbability < nx2_frac ) singlex = true;
209  else singlex = false;
210  else
211  if( xsizeProbability < nx1_frac ) singlex = true;
212  else singlex = false;
213 
214  if( hitbigy )
215  if( ysizeProbability < ny2_frac ) singley = true;
216  else singley = false;
217  else
218  if( ysizeProbability < ny1_frac ) singley = true;
219  else singley = false;
220 
221 
222 
223  // random multiplicity for alpha and beta
224  double qbinProbability = random->flatShoot();
225  for(int i = 0; i<4; ++i) {
226  nqbin = i;
227  if(qbinProbability < qbin_frac[i]) break;
228  }
229 
230  //Store interpolated pixel cluster profile
231  //BYSIZE, BXSIZE, const definition from SiPixelTemplate
232  float ytempl[41][BYSIZE] = {{0}}, xtempl[41][BXSIZE] = {{0}} ;
233  templ.ytemp(0, 40, ytempl);
234  templ.xtemp(0, 40, xtempl);
235 
236  std::vector<double> ytemp(BYSIZE);
237  for( int i=0; i<BYSIZE; ++i) {
238  ytemp[i]=(1.-yfrac)*ytempl[ybin][i]+yfrac*ytempl[ybin+1][i];
239  }
240 
241  std::vector<double> xtemp(BXSIZE);
242  for(int i=0; i<BXSIZE; ++i) {
243  xtemp[i]=(1.-xfrac)*xtempl[xbin][i]+xfrac*xtempl[xbin+1][i];
244  }
245 
246  //Pixel readout threshold
247  const float qThreshold = templ.s50()*2.0;
248 
249  //Cut away pixels below readout threshold
250  //For cluster lengths calculation
251  int offsetX1=0, offsetX2=0, offsetY1=0, offsetY2=0;
252  int firstY, lastY, firstX, lastX;
253  for( firstY = 0; firstY < BYSIZE; ++firstY ) {
254  bool yCluster = ytemp[firstY] > qThreshold ;
255  if( yCluster )
256  {
257  offsetY1 = BHY -firstY;
258  break;
259  }
260  }
261  for( lastY = firstY; lastY < BYSIZE; ++lastY )
262  {
263  bool yCluster = ytemp[lastY] > qThreshold ;
264  if( !yCluster )
265  {
266  lastY = lastY - 1;
267  offsetY2 = lastY - BHY;
268  break;
269  }
270  }
271 
272  for( firstX = 0; firstX < BXSIZE; ++firstX ) {
273  bool xCluster = xtemp[firstX] > qThreshold ;
274  if( xCluster ) {
275  offsetX1 = BHX - firstX;
276  break;
277  }
278  }
279  for( lastX = firstX; lastX < BXSIZE; ++ lastX ) {
280  bool xCluster = xtemp[lastX] > qThreshold ;
281  if( !xCluster ) {
282  lastX = lastX - 1;
283  offsetX2 = lastX - BHX;
284  break;
285  }
286  }
287 
288  bool edge, edgex, edgey;
289  // bool bigx, bigy;
290  unsigned int clslenx = offsetX1 + offsetX2 + 1;
291  unsigned int clsleny = offsetY1 + offsetY2 + 1;
292 
293  theClslenx = clslenx;
294  theClsleny = clsleny;
295 
296  int firstPixelInX = (int)mpx - offsetX1 ;
297  int firstPixelInY = (int)mpy - offsetY1 ;
298  int lastPixelInX = (int)mpx + offsetX2 ;
299  int lastPixelInY = (int)mpy + offsetY2 ;
300  firstPixelInX = (firstPixelInX >= 0) ? firstPixelInX : 0 ;
301  firstPixelInY = (firstPixelInY >= 0) ? firstPixelInY : 0 ;
302  lastPixelInX = (lastPixelInX < nrows ) ? lastPixelInX : nrows-1 ;
303  lastPixelInY = (lastPixelInY < ncolumns ) ? lastPixelInY : ncolumns-1;
304 
305  edgex = rectPixelTopology->isItEdgePixelInX( firstPixelInX ) || rectPixelTopology->isItEdgePixelInX( lastPixelInX );
306  edgey = rectPixelTopology->isItEdgePixelInY( firstPixelInY ) || rectPixelTopology->isItEdgePixelInY( lastPixelInY );
307  edge = edgex || edgey;
308 
309  // bigx = rectPixelTopology->isItBigPixelInX( firstPixelInX ) || rectPixelTopology->isItBigPixelInX( lastPixelInX );
310  // bigy = rectPixelTopology->isItBigPixelInY( firstPixelInY ) || rectPixelTopology->isItBigPixelInY( lastPixelInY );
311  bool hasBigPixelInX = rectPixelTopology->containsBigPixelInX( firstPixelInX, lastPixelInX );
312  bool hasBigPixelInY = rectPixelTopology->containsBigPixelInY( firstPixelInY, lastPixelInY );
313 
314  //Variables for SiPixelTemplate pixel hit error output
315  float sigmay, sigmax, sy1, sy2, sx1, sx2;
316  templ.temperrors(tempId, cotalpha, cotbeta, nqbin, sigmay, sigmax, sy1, sy2, sx1, sx2 );
317 
318  // define private mebers --> Errors
319  if( edge ) {
320  if( edgex && !edgey ) {
321  theErrorX = 23.0*microntocm;
322  theErrorY = 39.0*microntocm;
323  }
324  else if( !edgex && edgey ) {
325  theErrorX = 24.0*microntocm;
326  theErrorY = 96.0*microntocm;
327  }
328  else
329  {
330  theErrorX = 31.0*microntocm;
331  theErrorY = 90.0*microntocm;
332  }
333 
334  }
335  else {
336  if( singlex )
337  if ( hitbigx )
338  theErrorX = sx2*microntocm;
339  else
340  theErrorX = sx1*microntocm;
341  else theErrorX = sigmax*microntocm;
342  if( singley )
343  if( hitbigy )
344  theErrorY = sy2*microntocm;
345  else
346  theErrorY = sy1*microntocm;
347  else theErrorY = sigmay*microntocm;
348  }
349  theErrorZ = 1e-8; // 1 um means zero
351  // Local Error is 2D: (xx,xy,yy), square of sigma in first an third position
352  // as for resolution matrix
353  //
354 #ifdef FAMOS_DEBUG
355  std::cout << " Pixel Errors "
356  << "\talpha(x) = " << theErrorX
357  << "\tbeta(y) = " << theErrorY
358  << std::endl;
359 #endif
360  // Generate position
361  // get resolution histograms
362  int cotalphaHistBin = (int)( ( cotalpha - rescotAlpha_binMin ) / rescotAlpha_binWidth + 1 );
363  int cotbetaHistBin = (int)( ( cotbeta - rescotBeta_binMin ) / rescotBeta_binWidth + 1 );
364  // protection against out-of-range (undeflows and overflows)
365  if( cotalphaHistBin < 1 ) cotalphaHistBin = 1;
366  if( cotbetaHistBin < 1 ) cotbetaHistBin = 1;
367  if( cotalphaHistBin > (int)rescotAlpha_binN ) cotalphaHistBin = (int)rescotAlpha_binN;
368  if( cotbetaHistBin > (int)rescotBeta_binN ) cotbetaHistBin = (int)rescotBeta_binN;
369  //
370  unsigned int theXHistN;
371  unsigned int theYHistN;
372  if( !isForward ) {
373  if(edge)
374  {
375  theXHistN = cotalphaHistBin * 1000 + cotbetaHistBin * 10 + (nqbin+1);
376  theYHistN = theXHistN;
377  }
378  else
379  {
380  if(singlex)
381  {
382  if(hitbigx) theXHistN = 1 * 100000 + cotalphaHistBin * 100 + cotbetaHistBin ;
383  else theXHistN = 1 * 100000 + 1 * 1000 + cotalphaHistBin * 100 + cotbetaHistBin ;
384  }
385  else
386  {
387  if(hasBigPixelInX) theXHistN = 1 * 1000000 + 1 * 100000 + cotalphaHistBin * 1000 + cotbetaHistBin * 10 + (nqbin+1);
388  else theXHistN = 1 * 1000000 + 1 * 100000 + 1 * 10000 + cotalphaHistBin * 1000 + cotbetaHistBin * 10 + (nqbin+1);
389  }
390  if(singley)
391  {
392  if(hitbigy) theYHistN = 1 * 100000 + cotalphaHistBin * 100 + cotbetaHistBin ;
393  else theYHistN = 1 * 100000 + 1 * 1000 + cotalphaHistBin * 100 + cotbetaHistBin ;
394  }
395  else
396  {
397  if(hasBigPixelInY) theYHistN = 1 * 1000000 + 1 * 100000 + cotalphaHistBin * 1000 + cotbetaHistBin * 10 + (nqbin+1);
398  else theYHistN = 1 * 1000000 + 1 * 100000 + 1 * 10000 + cotalphaHistBin * 1000 + cotbetaHistBin * 10 + (nqbin+1);
399  }
400  }
401  }
402  else
403  {
404  if(edge)
405  {
406  theXHistN = cotalphaHistBin * 1000 + cotbetaHistBin * 10 + (nqbin+1);
407  theYHistN = theXHistN;
408  }
409  else
410  {
411  if( singlex )
412  if( hitbigx )
413  theXHistN = 100000 + cotalphaHistBin * 100 + cotbetaHistBin;
414  else
415  theXHistN = 100000 + 1000 + cotalphaHistBin * 100 + cotbetaHistBin;
416  else
417  theXHistN = 100000 + 10000 + cotalphaHistBin * 1000 + cotbetaHistBin * 10 + (nqbin+1);
418  if( singley )
419  if( hitbigy )
420  theYHistN = 100000 + cotalphaHistBin * 100 + cotbetaHistBin;
421  else
422  theYHistN = 100000 + 1000 + cotalphaHistBin * 100 + cotbetaHistBin;
423  else
424  theYHistN = 100000 + 10000 + cotalphaHistBin * 1000 + cotbetaHistBin * 10 + (nqbin+1);
425  }
426  }
427  unsigned int counter = 0;
428  do {
429  //
430  // Smear the hit Position
431  thePositionX = theXHistos[theXHistN]->generate();
432  thePositionY = theYHistos[theYHistN]->generate();
433  if( isForward ) thePositionY *= sign;
434  thePositionZ = 0.0; // set at the centre of the active area
435  //protect from empty resolution histograms
436  //if( thePositionX > 0.0799 ) thePositionX = 0;
437  //if( thePositionY > 0.0799 ) thePositionY = 0;
438  thePosition =
439  Local3DPoint(simHit.localPosition().x() + thePositionX ,
440  simHit.localPosition().y() + thePositionY ,
441  simHit.localPosition().z() + thePositionZ );
442 #ifdef FAMOS_DEBUG
443  std::cout << " Detector bounds: "
444  << "\t\tx = " << boundX
445  << "\ty = " << boundY
446  << std::endl;
447  std::cout << " Generated local position "
448  << "\tx = " << thePosition.x()
449  << "\ty = " << thePosition.y()
450  << std::endl;
451 #endif
452  counter++;
453  if(counter > 20) {
454  thePosition = Local3DPoint(simHit.localPosition().x(),
455  simHit.localPosition().y(),
456  simHit.localPosition().z());
457  break;
458  }
459  } while(fabs(thePosition.x()) > boundX || fabs(thePosition.y()) > boundY);
460 
461 }
int i
Definition: DBlmapReader.cc:9
#define BXSIZE
T y() const
Definition: PV2DBase.h:46
#define BYSIZE
virtual int ncolumns() const =0
LocalVector momentumAtEntry() const
The momentum of the track that produced the hit, at entry point.
Definition: PSimHit.h:47
T y() const
Definition: PV3DBase.h:63
virtual int nrows() const =0
std::map< unsigned, const SimpleHistogramGenerator * > theXHistos
bool interpolate(int id, float cotalpha, float cotbeta, float locBz)
virtual PixelGeomDetType & specificType() const
Local3DPoint localPosition() const
Definition: PSimHit.h:44
void xtemp(int fxbin, int lxbin, float xtemplate[41][BXSIZE])
#define BHX
T z() const
Definition: PV3DBase.h:64
void temperrors(int id, float cotalpha, float cotbeta, int qBin, float &sigmay, float &sigmax, float &sy1, float &sy2, float &sx1, float &sx2)
void ytemp(int fybin, int lybin, float ytemplate[41][BYSIZE])
float s50()
1/2 of the pixel threshold signal in electrons
LocalVector localDirection() const
Obsolete. Same as momentumAtEntry().unit(), for backward compatibility.
Definition: PSimHit.h:52
Point3DBase< float, LocalTag > Local3DPoint
Definition: LocalPoint.h:9
Vector3DBase unit() const
Definition: Vector3DBase.h:57
void qbin_dist(int id, float cotalpha, float cotbeta, float qbin_frac[4], float &ny1_frac, float &ny2_frac, float &nx1_frac, float &nx2_frac)
std::map< unsigned, const SimpleHistogramGenerator * > theYHistos
double flatShoot(double xmin=0.0, double xmax=1.0) const
Definition: RandomEngine.h:30
#define BHY
tuple cout
Definition: gather_cfg.py:121
T x() const
Definition: PV2DBase.h:45
T x() const
Definition: PV3DBase.h:62

Member Data Documentation

bool SiPixelGaussianSmearingRecHitConverterAlgorithm::isForward
private
edm::ParameterSet SiPixelGaussianSmearingRecHitConverterAlgorithm::pset_
private
const RandomEngine* SiPixelGaussianSmearingRecHitConverterAlgorithm::random
private
double SiPixelGaussianSmearingRecHitConverterAlgorithm::rescotAlpha_binMin
private
unsigned int SiPixelGaussianSmearingRecHitConverterAlgorithm::rescotAlpha_binN
private
double SiPixelGaussianSmearingRecHitConverterAlgorithm::rescotAlpha_binWidth
private
double SiPixelGaussianSmearingRecHitConverterAlgorithm::rescotBeta_binMin
private
unsigned int SiPixelGaussianSmearingRecHitConverterAlgorithm::rescotBeta_binN
private
double SiPixelGaussianSmearingRecHitConverterAlgorithm::rescotBeta_binWidth
private
int SiPixelGaussianSmearingRecHitConverterAlgorithm::resqbin_binMin
private
unsigned int SiPixelGaussianSmearingRecHitConverterAlgorithm::resqbin_binN
private
int SiPixelGaussianSmearingRecHitConverterAlgorithm::resqbin_binWidth
private
int SiPixelGaussianSmearingRecHitConverterAlgorithm::tempId
private
SiPixelTemplate SiPixelGaussianSmearingRecHitConverterAlgorithm::templ
private
unsigned int SiPixelGaussianSmearingRecHitConverterAlgorithm::theClslenx
private
unsigned int SiPixelGaussianSmearingRecHitConverterAlgorithm::theClsleny
private
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
GeomDetType::SubDetector SiPixelGaussianSmearingRecHitConverterAlgorithm::thePixelPart
private
TFile* SiPixelGaussianSmearingRecHitConverterAlgorithm::thePixelResolutionFile1
private
TFile* SiPixelGaussianSmearingRecHitConverterAlgorithm::thePixelResolutionFile2
private
std::string SiPixelGaussianSmearingRecHitConverterAlgorithm::thePixelResolutionFileName1
private
std::string SiPixelGaussianSmearingRecHitConverterAlgorithm::thePixelResolutionFileName2
private
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
std::map<unsigned,const SimpleHistogramGenerator*> SiPixelGaussianSmearingRecHitConverterAlgorithm::theYHistos
private
bool SiPixelGaussianSmearingRecHitConverterAlgorithm::useCMSSWPixelParameterization
private