#include <PixelTemplateSmearerBase.h>

Inheritance diagram for PixelTemplateSmearerBase:

Classes
struct	MergeGroup

Public Member Functions
void	beginRun (edm::Run const &run, const edm::EventSetup &eventSetup, const SiPixelTemplateDBObject *pixelTemplateDBObjectPtr, std::vector< SiPixelTemplateStore > &tempStoreRef) override

bool	hitsMerge (const PSimHit &simHit1, const PSimHit &simHit2) const

	PixelTemplateSmearerBase (const std::string &name, const edm::ParameterSet &config, edm::ConsumesCollector &consumesCollector)

TrackingRecHitProductPtr	process (TrackingRecHitProductPtr product) const override

TrackingRecHitProductPtr	processMergeGroups (std::vector< MergeGroup * > &mergeGroups, TrackingRecHitProductPtr product, const PixelGeomDetUnit detUnit, const double boundX, const double boundY, RandomEngineAndDistribution const random) const

TrackingRecHitProductPtr	processUnmergedHits (std::vector< TrackingRecHitProduct::SimHitIdPair > &unmergedHits, TrackingRecHitProductPtr product, const PixelGeomDetUnit detUnit, const double boundX, const double boundY, RandomEngineAndDistribution const random) const

FastSingleTrackerRecHit	smearHit (const PSimHit &simHit, const PixelGeomDetUnit detUnit, const double boundX, const double boundY, RandomEngineAndDistribution const ) const

FastSingleTrackerRecHit	smearMergeGroup (MergeGroup mg, const PixelGeomDetUnit detUnit, const double boundX, const double boundY, const RandomEngineAndDistribution *random) const

	~PixelTemplateSmearerBase () override

Public Member Functions inherited from TrackingRecHitAlgorithm
virtual void	beginEvent (edm::Event &event, const edm::EventSetup &eventSetup)

virtual void	beginStream (const edm::StreamID &id)

virtual void	endEvent (edm::Event &event, const edm::EventSetup &eventSetup)

virtual void	endStream ()

const TrackerGeometry &	getMisalignedGeometry () const

const std::string &	getName () const

const RandomEngineAndDistribution &	getRandomEngine () const

const std::string &	getSelectionString () const

const TrackerGeometry &	getTrackerGeometry () const

const TrackerTopology &	getTrackerTopology () const

	TrackingRecHitAlgorithm (const std::string &name, const edm::ParameterSet &config, edm::ConsumesCollector &consumesCollector)

virtual	~TrackingRecHitAlgorithm ()

Protected Attributes
bool	isBarrel

bool	mergeHitsOn = false

const SiPixelTemplateDBObject *	pixelTemplateDBObject_ = 0

int	templateId = -1

std::string	theBigPixelResolutionFileName

std::shared_ptr< PixelResolutionHistograms >	theBigPixelResolutions

std::string	theEdgePixelResolutionFileName

std::shared_ptr< PixelResolutionHistograms >	theEdgePixelResolutions

std::unique_ptr< TFile >	theMergedPixelResolutionXFile

std::string	theMergedPixelResolutionXFileName

std::unique_ptr< TFile >	theMergedPixelResolutionYFile

std::string	theMergedPixelResolutionYFileName

std::unique_ptr< TFile >	theMergingProbabilityFile

std::string	theMergingProbabilityFileName

std::vector< SiPixelTemplateStore >	thePixelTemp_

std::vector< SiPixelTemplateStore > &	thePixelTempRef = thePixelTemp_

std::string	theRegularPixelResolutionFileName

std::shared_ptr< PixelResolutionHistograms >	theRegularPixelResolutions

Detailed Description

Definition at line 40 of file PixelTemplateSmearerBase.h.

Constructor & Destructor Documentation

PixelTemplateSmearerBase::PixelTemplateSmearerBase	(	const std::string &	name,
		const edm::ParameterSet &	config,
		edm::ConsumesCollector &	consumesCollector
	)

explicit

Definition at line 47 of file PixelTemplateSmearerBase.cc.

References Exception, edm::ParameterSet::exists(), edm::FileInPath::fullPath(), edm::ParameterSet::getParameter(), isBarrel, mergeHitsOn, SiPixelTemplate::pushfile(), mps_update::status, templateId, theBigPixelResolutionFileName, theBigPixelResolutions, theEdgePixelResolutionFileName, theEdgePixelResolutions, theMergedPixelResolutionXFile, theMergedPixelResolutionXFileName, theMergedPixelResolutionYFile, theMergedPixelResolutionYFileName, theMergingProbabilityFile, theMergingProbabilityFileName, thePixelTemp_, theRegularPixelResolutionFileName, theRegularPixelResolutions, and funct::true.

  :
   TrackingRecHitAlgorithm(name,config,consumesCollector)
 {
     //--- Basic stuff
     mergeHitsOn = config.getParameter<bool>("MergeHitsOn");
     isBarrel    = config.getParameter<bool> ( "isBarrel" );
     int detType = (isBarrel) ? 1 : 0;    // 1 for barrel, 0 for forward  (or could we just promote bool into int...?)
 
     //--- Resolution file names.
     theBigPixelResolutionFileName     = config.getParameter<string>( "BigPixelResolutionFile" );
     theBigPixelResolutionFileName     = edm::FileInPath( theBigPixelResolutionFileName ).fullPath();
 
     theEdgePixelResolutionFileName    = config.getParameter<string>( "EdgePixelResolutionFile" );
     theEdgePixelResolutionFileName    = edm::FileInPath( theEdgePixelResolutionFileName ).fullPath();
 
     theRegularPixelResolutionFileName = config.getParameter<string>( "RegularPixelResolutionFile" );
     theRegularPixelResolutionFileName = edm::FileInPath( theRegularPixelResolutionFileName ).fullPath();
 
 
     //--- Create the resolution histogram objects, which will load the histograms
     //    and initialize random number generators.
     //
     int status = 0;
     theRegularPixelResolutions = 
       std::make_shared<PixelResolutionHistograms>( theRegularPixelResolutionFileName, "" );
     if (( status = theRegularPixelResolutions->status()) != 0 ) {
       throw cms::Exception("PixelTemplateSmearerBase:")
     << " constructing PixelResolutionHistograms file " << theRegularPixelResolutionFileName
     << " failed with status = " << status << std::endl;
     }
     
     theBigPixelResolutions =
       std::make_shared<PixelResolutionHistograms>( theBigPixelResolutionFileName, "", detType, (!isBarrel), false, true );  // can miss qBin
     if (( status = theBigPixelResolutions->status()) != 0 ) {
       throw cms::Exception("PixelTemplateSmearerBase:")
     << " constructing PixelResolutionHistograms file " << theBigPixelResolutionFileName
     << " failed with status = " << status << std::endl;
     }
     
     theEdgePixelResolutions =
       std::make_shared<PixelResolutionHistograms>( theEdgePixelResolutionFileName, "", detType, false, true, true );  // can miss both single & qBin
     if (( status = theEdgePixelResolutions->status()) != 0 ) {
       throw cms::Exception("PixelTemplateSmearerBase:")
     << " constructing PixelResolutionHistograms file " << theEdgePixelResolutionFileName
     << " failed with status = " << status << std::endl;
     }
 
 
     
     //--- Merging info.
     theMergingProbabilityFileName     = config.getParameter<string>( "MergingProbabilityFile" );
     theMergingProbabilityFileName     = edm::FileInPath( theMergingProbabilityFileName ).fullPath();
     theMergingProbabilityFile         = std::make_unique<TFile>( theMergingProbabilityFileName.c_str(), "READ");
 
     theMergedPixelResolutionXFileName = config.getParameter<string>( "MergedPixelResolutionXFile" );
     theMergedPixelResolutionXFileName = edm::FileInPath( theMergedPixelResolutionXFileName ).fullPath();
     theMergedPixelResolutionXFile     = std::make_unique<TFile>( theMergedPixelResolutionXFileName.c_str(), "READ");
 
     theMergedPixelResolutionYFileName = config.getParameter<string>( "MergedPixelResolutionYFile" );
     theMergedPixelResolutionYFileName = edm::FileInPath( theMergedPixelResolutionYFileName ).fullPath();
     theMergedPixelResolutionYFile     = std::make_unique<TFile>( theMergedPixelResolutionYFileName.c_str(), "READ");
 
 
     // const SiPixelTemplateDBObject & dbobject;
     // const SiPixelTemplateDBObject dbobject;        // dummy, just to make it compile &&&
 
     //--- Load the templates.
     if ( config.exists("templateId") ) {
       //--- Load template with ID=templateId from a local ascii file.
       templateId  = config.getParameter<int> ( "templateId" );
       if ( templateId > 0 ) {
     if ( !SiPixelTemplate::pushfile(templateId, thePixelTemp_) ) {
       throw cms::Exception("PixelTemplateSmearerBase:")
         <<"SiPixel Template " << templateId << " Not Loaded Correctly!" << std::endl;
     }
       }
     }
     
     //--- Else... The templates will be loaded from the DB... 
     //    (They are needed for data and full sim MC, so in a production FastSim
     //    run, everything should already be in the DB.)
     //
     //    But note that we can do it only at the beginning of the
     //    event.  So nothing happens now.
 }

PixelTemplateSmearerBase::~PixelTemplateSmearerBase ( )

override

Definition at line 138 of file PixelTemplateSmearerBase.cc.

References thePixelTemp_, and x.

 {
   //--- Delete the templates. This is safe even if thePixelTemp_ vector is empty.
   for (auto x : thePixelTemp_) x.destroy();
 }

Member Function Documentation

void PixelTemplateSmearerBase::beginRun	(	edm::Run const &	run,
		const edm::EventSetup &	eventSetup,
		const SiPixelTemplateDBObject *	pixelTemplateDBObjectPtr,
		std::vector< SiPixelTemplateStore > &	tempStoreRef
	)

overridevirtual

Reimplemented from TrackingRecHitAlgorithm.

Definition at line 151 of file PixelTemplateSmearerBase.cc.

References pixelTemplateDBObject_, templateId, and thePixelTempRef.

 {
   //--- Check if we need to use the template from the DB (namely if
   //    id == -1).  Otherwise the template has already been loaded from
   //    the ascii file in constructor, and thePixelTempRef wakes up
   //    pointing to thePixelTemp_, so then we use our own store.
   //
   if ( templateId == -1 ) {
     thePixelTempRef = tempStoreRef;        // we use the store from TrackingRecHitProducer
     pixelTemplateDBObject_ = pixelTemplateDBObjectPtr;   // needed for template<-->DetId map.
   }
   
 
   //--- Commented code below (the DB interface) should say here, in case we need it:
   // edm::ESHandle<SiPixelTemplateDBObject> templateDBobject;
   // eventSetup.get<SiPixelTemplateDBObjectESProducerRcd>().get(templateDBobject);
   // pixelTemplateDBObject_ = templateDBobject.product();
   
   // //--- Now that we have the DB object, load the correct templates from the DB.  
   // //    (They are needed for data and full sim MC, so in a production FastSim
   // //    run, everything should already be in the DB.)
   // if ( !SiPixelTemplate::pushfile( *pixelTemplateDBObject_ , thePixelTemp_) ) {
   //   throw cms::Exception("PixelTemplateSmearerPlugin:")
   //    <<"SiPixel Template " << templateId << " Not Loaded Correctly!"<<endl;
   // }
 }

bool PixelTemplateSmearerBase::hitsMerge	(	const PSimHit &	simHit1,
		const PSimHit &	simHit2
	)		const

Definition at line 1081 of file PixelTemplateSmearerBase.cc.

References PSimHit::localPosition(), cmsBatch::log, PSimHit::momentumAtEntry(), funct::pow(), TtFullHadEvtBuilder_cfi::prob, mathSSE::sqrt(), theMergingProbabilityFile, Vector3DBase< T, FrameTag >::unit(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by process().

 {
     LocalVector localDir = simHit1.momentumAtEntry().unit();
     float locy1 = localDir.y();
     float locz1 = localDir.z();
     float cotbeta = locy1/locz1;
     float loceta = fabs(-log((double)(-cotbeta+sqrt((double)(1.+cotbeta*cotbeta)))));
 
     const Local3DPoint lp1 = simHit1.localPosition();
     const Local3DPoint lp2 = simHit2.localPosition();
     float lpy1 = lp1.y();
     float lpx1 = lp1.x();
     float lpy2 = lp2.y();
     float lpx2 = lp2.x();
     float locdis = 10000.* sqrt(pow(lpx1 - lpx2, 2) + pow(lpy1 - lpy2, 2));
     TH2F * probhisto = (TH2F*)theMergingProbabilityFile->Get("h2bc");
     float prob = probhisto->GetBinContent(probhisto->GetXaxis()->FindFixBin(locdis),probhisto->GetYaxis()->FindFixBin(loceta));
     return prob > 0;
 }

TrackingRecHitProductPtr PixelTemplateSmearerBase::process ( TrackingRecHitProductPtr product ) const

overridevirtual

Reimplemented from TrackingRecHitAlgorithm.

Definition at line 189 of file PixelTemplateSmearerBase.cc.

References Exception, TrackingRecHitAlgorithm::getRandomEngine(), TrackingRecHitAlgorithm::getTrackerGeometry(), edmIntegrityCheck::group, PixelTemplateSmearerBase::MergeGroup::group, hitsMerge(), mps_fire::i, TrackerGeometry::idToDetUnit(), gen::k, Bounds::length(), mergeHitsOn, processMergeGroups(), processUnmergedHits(), randomEngine, edm::second(), trackerHits::simHits, PixelTemplateSmearerBase::MergeGroup::smearIt, GeomDet::surface(), and Bounds::width().

 {
     std::vector<std::pair<unsigned int,const PSimHit*>> & simHitIdPairs = product->getSimHitIdPairs();
     std::vector<const PSimHit*> simHits(simHitIdPairs.size());
     for (unsigned int ihit = 0; ihit<simHitIdPairs.size();++ihit)
     {
         simHits[ihit]=simHitIdPairs[ihit].second;
     }
 
     RandomEngineAndDistribution const & randomEngine = getRandomEngine();
 
 
     const GeomDet* geomDet = getTrackerGeometry().idToDetUnit(product->getDetId());
     const PixelGeomDetUnit * pixelGeomDet = dynamic_cast< const PixelGeomDetUnit* >( geomDet );
     if (pixelGeomDet == nullptr)
     {
         throw cms::Exception("FastSimulation/TrackingRecHitProducer") << "The GeomDetUnit is not a PixelGeomDetUnit.  This should never happen!";
     }
     const BoundPlane& theDetPlane = pixelGeomDet->surface();
     const Bounds& theBounds = theDetPlane.bounds();
     const double boundX = theBounds.width()/2.;
     const double boundY = theBounds.length()/2.;
 
 
     std::vector<TrackingRecHitProduct::SimHitIdPair> listOfUnmergedHits; 
     std::vector< MergeGroup* > listOfMergeGroups;
     int nHits = simHits.size();          
     
     // fixed size array, 0 if hit is unmerged             
     MergeGroup* mergeGroupByHit[nHits];             
 
 
     if (nHits == 0)
     {
         return product;
     }
     else if (nHits == 1)
     {
         listOfUnmergedHits.push_back(simHitIdPairs[0]);
     }
     else 
     {
         if ( mergeHitsOn )
         {
    
             for (int i = 0; i < nHits; ++i )
             {
                 //initialize this cell to a NULL pointer here
                 mergeGroupByHit[i] = nullptr;
             }
             for ( int i = 0; i < nHits-1; ++i )
             {
                 for ( int j = i+1 ; j < nHits; ++j ) {
 
                     //--- Calculate the distance between hits i and j:
                     bool merged = hitsMerge( *simHitIdPairs[i].second, *simHitIdPairs[j].second );
 
 
                     if ( merged )
                     {
                         // First, check if the other guy (j) is in some merge group already
                         if ( mergeGroupByHit[j] != nullptr ) 
                         {
                             if (mergeGroupByHit[i] == nullptr ) 
                             {
                                 mergeGroupByHit[i] = mergeGroupByHit[j];
                                 mergeGroupByHit[i]->group.push_back(simHitIdPairs[i]);
                                 mergeGroupByHit[i]->smearIt = true;
                             }
                             else
                             {
                                 if (mergeGroupByHit[i] != mergeGroupByHit[j])
                                 {
                                     for (auto hit_it = mergeGroupByHit[j]->group.begin(); hit_it != mergeGroupByHit[j]->group.end(); ++hit_it)
                                     {
                                         mergeGroupByHit[i]->group.push_back( *hit_it );
                                         mergeGroupByHit[i]->smearIt = true;
                                     }
 
                                     // Step 2: iterate over all hits, replace mgbh[j] by mgbh[i] (so that nobody points to i)
                                     MergeGroup * mgbhj = mergeGroupByHit[j];                    
                                     for ( int k = 0; k < nHits; ++k )
                                     {
                                             if ( mgbhj == mergeGroupByHit[k])
                                             {
                                                 // Hit k also uses the same merge group, tell them to switch to mgbh[i]
                                                 mergeGroupByHit[k] = mergeGroupByHit[i];
                         }
                                     }
                                     mgbhj->smearIt = false;
                                     mergeGroupByHit[i]->smearIt = true;
 
                                     //  Step 3 would have been to delete mgbh[j]... however, we'll do that at the end anyway.
                                     //  The key was to prevent mgbh[j] from being accessed further, and we have done that,
                                     //  since now no mergeGroupByHit[] points to mgbhj any more.  Note that the above loop
                                     //  also set mergeGroupByHit[i] = mergeGroupByHit[j], too. 
                                 }
                             }
                         }
                         else 
                         { 
                             // j is not merged.  Check if i is merged with another hit yet.
                             //
                             if ( mergeGroupByHit[i] == nullptr )
                             {
                                 // This is the first time we realized i is merged with any
                                 // other hit.  Create a new merge group for i and j
                                 mergeGroupByHit[i] = new MergeGroup();
                                 listOfMergeGroups.push_back( mergeGroupByHit[i] );   // keep track of it
                                 //
                                 // Add hit i as the first to its own merge group
                                 // (simHits[i] is a const pointer to PSimHit).
                                 mergeGroupByHit[i]->group.push_back( simHitIdPairs[i] );
                                 mergeGroupByHit[i]->smearIt = true;
                             }
                             //--- Add hit j as well
                             mergeGroupByHit[i]->group.push_back( simHitIdPairs[j] );
                             mergeGroupByHit[i]->smearIt = true;
                             
                             mergeGroupByHit[j] = mergeGroupByHit[i];
 
                         } // --- end of else if ( j has merge group )
 
                     } //--- end of if (merged)
 
                 } //--- end of loop over j
 
                 //--- At this point, there are two possibilities.  Either hit i
                 //    was already chosen to be merged with some hit prior to it,
                 //    or the loop over j found another merged hit.  In either
                 //    case, if mergeGroupByHit[i] is empty, then the hit is
                 //    unmerged.
                 //
                 if ( mergeGroupByHit[i] == nullptr )
                 {
                     //--- Keep track of it.
                     listOfUnmergedHits.push_back( simHitIdPairs[i] );
                 }
             } //--- end of loop over i
         } // --- end of if (mergeHitsOn)
         else
         {
             // Now we've turned off hit merging, so all hits should be pushed
             // back to listOfUnmergedHits
             for (int i = 0; i < nHits; ++i )
             {
                 listOfUnmergedHits.push_back( simHitIdPairs[i] );
             }
         }
     } // --- end of if (nHits == 1) else {...}
 
 
     //--- We now have two lists: a list of hits that are unmerged, and
     //    the list of merge groups.  Process each separately.
     //
     product = processUnmergedHits( listOfUnmergedHits, product,
     pixelGeomDet, boundX, boundY,
     &randomEngine );
 
 
     product = processMergeGroups(  listOfMergeGroups,  product,
     pixelGeomDet, boundX, boundY,
     &randomEngine );
 
     //--- We're done with this det unit, and ought to clean up used
     //    memory.  We don't own the PSimHits, and the vector of
     //    listOfUnmergedHits simply goes out of scope.  However, we
     //    created the MergeGroups and thus we need to get rid of them.
     //
     for (auto mg_it = listOfMergeGroups.begin(); mg_it != listOfMergeGroups.end(); ++mg_it )
     {
         delete *mg_it;    // each MergeGroup is deleted; its ptrs to PSimHits we do not own...
     }
 
     return product;
 }

TrackingRecHitProductPtr PixelTemplateSmearerBase::processMergeGroups	(	std::vector< MergeGroup * > &	mergeGroups,
		TrackingRecHitProductPtr	product,
		const PixelGeomDetUnit *	detUnit,
		const double	boundX,
		const double	boundY,
		RandomEngineAndDistribution const *	random
	)		const

Definition at line 801 of file PixelTemplateSmearerBase.cc.

References rpcPointValidation_cfi::recHit, and smearMergeGroup().

Referenced by process(), and processUnmergedHits().

 {
     for ( auto mg_it = mergeGroups.begin(); mg_it != mergeGroups.end(); ++mg_it) 
     {
         if ((*mg_it)->smearIt)
         {
             FastSingleTrackerRecHit recHit = smearMergeGroup( *mg_it, detUnit, boundX, boundY, random );
             product->addRecHit( recHit, (*mg_it)->group);
         }
     }
     return product;
 }

TrackingRecHitProductPtr PixelTemplateSmearerBase::processUnmergedHits	(	std::vector< TrackingRecHitProduct::SimHitIdPair > &	unmergedHits,
		TrackingRecHitProductPtr	product,
		const PixelGeomDetUnit *	detUnit,
		const double	boundX,
		const double	boundY,
		RandomEngineAndDistribution const *	random
	)		const

Definition at line 779 of file PixelTemplateSmearerBase.cc.

References processMergeGroups(), rpcPointValidation_cfi::recHit, and smearHit().

Referenced by process(), and smearHit().

 {
     for (auto simHitIdPair: unmergedHits)
     {
         FastSingleTrackerRecHit recHit = smearHit( *simHitIdPair.second, detUnit, boundX, boundY, random );
         product->addRecHit( recHit ,{simHitIdPair});
     }
     return product;
 }

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

Definition at line 370 of file PixelTemplateSmearerBase.cc.

Referenced by processUnmergedHits().

 {
 
     //--- 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();  // don't need .unit(), we will take the ratio
     float locx = localDir.x();
     float locy = localDir.y();
     float locz = localDir.z();
 
     //--- cotangent of local angles \alpha and \beta.
     //    alpha: angle with respect to local x axis in local (x,z) plane
     //    beta: angle with respect to local y axis in local (y,z) plane
     //
     float cotalpha = locx/locz;
     float cotbeta = locy/locz;
 
     //--- Save the original signs of cot\alpha and cot\beta
     int signOfCotalpha = (cotalpha < 0) ? -1 : 1;   // sign(cotalpha);
     int signOfCotbeta  = (cotbeta  < 0) ? -1 : 1;   // sign(cotbeta);
     //
     //--- Use absolute values to find the templates from the list
     cotalpha *= signOfCotalpha;  // = abs(cotalpha)
     cotbeta  *= signOfCotbeta;   // = abs(cotbeta)
 
     LogDebug ("SmearHit") << "LocalVector=" << locx << "," << locy << "," << locz
               << "   momentum=" << localDir.mag() 
               << "   cotalpha=" << cotalpha << ",  cotbeta=" << cotbeta;
 
     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;
 
     const MeasurementPoint mpCenter(pixelCenterX, pixelCenterY);
     //Transform the center of the struck pixel back into local position
     const Local3DPoint lpCenter = rectPixelTopology->localPosition( mpCenter );
 
     //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]  // &&& Why limit of 39?
     if( ybin < 0 )    ybin = 0;
     if( ybin > 39 )   ybin = 39;
     if( xbin < 0 )    xbin = 0;
     if( xbin > 39 )   xbin = 39; 
 
 
     int ID = templateId;
     if ( templateId == -1 ) {
       // We have loaded the whole template set from the DB,
       // so ask the DB object to find us the right one.
       ID = pixelTemplateDBObject_->getTemplateID( detUnit->geographicalId() );  // need uint32_t detid
       //                    theDetParam.theDet->geographicalId());
     }
 
     //--- Make the template object
     SiPixelTemplate templ( thePixelTempRef );
 
     //--- Produce the template that corresponds to our local angles.
     templ.interpolate( ID, cotalpha, cotbeta );
 
     //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;
     templ.qbin_dist( ID, 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 );  // pixel we hit in x
     bool hitbigy = rectPixelTopology->isItBigPixelInY( (int)mpy );  // pixel we hit in y
 
     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;
         }
     }
 
 
     //--- Prepare to return results
     Local3DPoint thePosition;
     double       theShiftInX; 
     double       theShiftInY; 
     double       theShiftInZ; 
     LocalError   theError;     
     double       theErrorX;    
     double       theErrorY;    
 
     //------------------------------
     //  Check if the cluster is near an edge.  If it protrudes
     //  outside the edge of the sensor, the truncate it and it will
     //  get significantly messed up.
     //------------------------------
     bool edge, edgex, edgey;
     //  bool bigx, bigy;
 
     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(
         ID, cotalpha, cotbeta, nqbin,          // inputs
         sigmay, sigmax, sy1, sy2, sx1, sx2     // outputs
     );
 
     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;
         }
     }
 
 
     //add misalignment error
     const TrackerGeomDet* misalignmentDetUnit = getMisalignedGeometry().idToDet(detUnit->geographicalId());
     const LocalError& alignmentError = misalignmentDetUnit->localAlignmentError();
     if (alignmentError.valid())
     {
         theError = LocalError( 
             theErrorX*theErrorX+alignmentError.xx(), 
             alignmentError.xy(), 
             theErrorY*theErrorY+alignmentError.yy()
         );
     }
     else
     {
         theError = LocalError( 
             theErrorX*theErrorX, 
             0.0,
             theErrorY*theErrorY
         );
     }
 
     // Local Error is 2D: (xx,xy,yy), square of sigma in first an third position 
     // as for resolution matrix
 
     //--- Next, we need to generate the smeared position.  First we need to figure
     //    out which kind of histograms we are supposed to use for this particular hit.
     //    These are pointers to the set of histograms used to generate the rec hit 
     //    positions.  (We need to handle X and Y separately.)
     shared_ptr<PixelResolutionHistograms> resHistsX = nullptr;
     shared_ptr<PixelResolutionHistograms> resHistsY = nullptr;
 
     
     if (edge) {
       resHistsX = resHistsY = theEdgePixelResolutions;
       singlex = singley = false;          // no single resolutions for Edge
     }
     else {
       //--- Decide resolution histogram set for X
       if ( (singlex && hitbigx) || (isBarrel && hasBigPixelInX) ) {
     resHistsX = theBigPixelResolutions;
       }
       else {
     resHistsX = theRegularPixelResolutions;
       }
       //--- Decide resolution histogram set for Y
       if ( (singley && hitbigy) || (isBarrel && hasBigPixelInY) ) {
     resHistsY = theBigPixelResolutions;
       }
       else {
     resHistsY = theRegularPixelResolutions;
       }
     }
 
     //--- Get generators, separately for X and for Y.
     const SimpleHistogramGenerator * xgen
       = resHistsX->getGeneratorX( cotalpha, cotbeta, nqbin, singlex );
     const SimpleHistogramGenerator * ygen
       = resHistsY->getGeneratorY( cotalpha, cotbeta, nqbin, singley );
 
     //--- Check if we found a histogram.  If nullptr, then throw up.
     if ( !xgen || !ygen ) {
       throw cms::Exception("FastSimulation/TrackingRecHitProducer")
     << "Histogram (cot\alpha=" << cotalpha
     <<", cot\beta="<< cotbeta << ", nQbin=" << nqbin
     << ") was not found for PixelTemplateSmearer. Check if the smearing resolution histogram exists.";
     }
 
     
     //--- Smear the hit Position.  We do it in the do-while loop in order to
     //--- allow multiple tries, in case we generate a rec hit which is outside
     //--- of the boundaries of the sensor.
     unsigned int retry = 0;
     
     do 
     {
         // Generate the position (x,y of the rec hit).
         theShiftInX = xgen->generate(random);
         theShiftInY = ygen->generate(random);
 
     // Now multiply by the sign of the cotangent of appropriate angle
     theShiftInX *= signOfCotalpha;
     theShiftInY *= signOfCotbeta;
 
     theShiftInZ = 0.0; // set to the mid-plane of the sensor.
 
         thePosition = Local3DPoint(
             simHit.localPosition().x() + theShiftInX, 
             simHit.localPosition().y() + theShiftInY, 
             simHit.localPosition().z() + theShiftInZ
         );
         retry++;
         if (retry > 10) 
         {
             // If we tried to generate thePosition, and it's out of the bounds
             // for 10 times, then take and return the simHit's location.
             thePosition = Local3DPoint(
                 simHit.localPosition().x(), 
                 simHit.localPosition().y(), 
                 simHit.localPosition().z()
             );
             break;
         }
     } while(fabs(thePosition.x()) > boundX  || fabs(thePosition.y()) > boundY);
     
     
     FastSingleTrackerRecHit recHit(
         thePosition,
         theError,
         *detUnit,
         fastTrackerRecHitType::siPixel
     );
     return recHit;
 }

FastSingleTrackerRecHit PixelTemplateSmearerBase::smearMergeGroup	(	MergeGroup *	mg,
		const PixelGeomDetUnit *	detUnit,
		const double	boundX,
		const double	boundY,
		const RandomEngineAndDistribution *	random
	)		const

Definition at line 826 of file PixelTemplateSmearerBase.cc.

References BXSIZE, BYSIZE, RandomEngineAndDistribution::flatShoot(), objects.autophobj::float, SimpleHistogramGenerator::generate(), GeomDet::geographicalId(), SiPixelTemplateDBObject::getTemplateID(), PixelTemplateSmearerBase::MergeGroup::group, mps_fire::i, createfilelist::int, SiPixelTemplate::interpolate(), PSimHit::localPosition(), microntocm, PSimHit::momentumAtEntry(), pixelTemplateDBObject_, SiPixelTemplate::qbin_dist(), rpcPointValidation_cfi::recHit, rpcPointValidation_cfi::simHit, fastTrackerRecHitType::siPixel, SiPixelTemplate::temperrors(), templateId, theMergedPixelResolutionXFile, theMergedPixelResolutionYFile, thePixelTempRef, Vector3DBase< T, FrameTag >::unit(), PV3DBase< T, PVType, FrameType >::x(), xsize, SiPixelTemplate::xtemp(), PV3DBase< T, PVType, FrameType >::y(), ysize, SiPixelTemplate::ytemp(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by processMergeGroups().

 {
 
     float loccx = 0;
     float loccy = 0;
     float loccz = 0;
     float nHit = 0;
     float locpx = 0;
     float locpy = 0;
     float locpz = 0;
 
     for (auto hit_it = mg->group.begin(); hit_it != mg->group.end(); ++hit_it) 
     {
         const PSimHit simHit = *hit_it->second; 
         //getting local momentum and adding all of the hits' momentums up
         LocalVector localDir = simHit.momentumAtEntry().unit();
         loccx += localDir.x();
         loccy += localDir.y();
         loccz += localDir.z();
         //getting local position and adding all of the hits' positions up
         const Local3DPoint lpos = simHit.localPosition();
         locpx += lpos.x();
         locpy += lpos.y();
         locpz += lpos.z();
         //counting how many sim hits are in the merge group
         nHit += 1;
     }
     //averaging the momentums by diving momentums added up/number of hits
     float locx = loccx/nHit;
     float locy = loccy/nHit;
     float locz = loccz/nHit;
 
     //--- cotangent of local angles \alpha and \beta.
     //    alpha: angle with respect to local x axis in local (x,z) plane
     //    beta: angle with respect to local y axis in local (y,z) plane
     //
     float cotalpha = locx/locz;
     float cotbeta = locy/locz;
 
     //--- Save the original signs of cot\alpha and cot\beta
     int signOfCotalpha = (cotalpha < 0) ? -1 : 1;   // sign(cotalpha);
     int signOfCotbeta  = (cotbeta  < 0) ? -1 : 1;   // sign(cotbeta);
     //
     //--- Use absolute values to find the templates from the list
     cotalpha *= signOfCotalpha;  // = abs(cotalpha)
     cotbeta  *= signOfCotbeta;   // = abs(cotbeta)
 
     float lpx = locpx/nHit;
     float lpy = locpy/nHit;
     float lpz = locpz/nHit;
 
     //Get the relative position of struck point to the center of the struck pixel
     float xtrk = lpx;
     float ytrk = lpy;
     //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; 
 
     int ID = templateId;
     if ( templateId == -1 ) {
       // We have loaded the whole template set from the DB,
       // so ask the DB object to find us the right one.
       ID = pixelTemplateDBObject_->getTemplateID( detUnit->geographicalId() );  // need uint32_t detid
       //                    theDetParam.theDet->geographicalId());
     }
 
     //--- Make the template object
     SiPixelTemplate templ( thePixelTempRef );
 
     //--- Produce the template that corresponds to our local angles.
     templ.interpolate( ID, cotalpha, cotbeta );
 
     // 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;
     templ.qbin_dist( ID, cotalpha, cotbeta, qbin_frac, ny1_frac, ny2_frac, nx1_frac, nx2_frac );
     int  nqbin;
 
     //  double xsizeProbability = random->flatShoot();
     //double ysizeProbability = random->flatShoot();
     bool hitbigx = false;
     bool hitbigy = 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];
     }
 
     //--- Prepare to return results
     Local3DPoint thePosition;  
     double       theShiftInX; 
     double       theShiftInY; 
     double       theShiftInZ; 
     LocalError   theError;     
     double       theErrorX;    
     double       theErrorY;    
 
 
     //------------------------------
     //  Check if the cluster is near an edge.  If it protrudes
     //  outside the edge of the sensor, the truncate it and it will
     //  get significantly messed up.
     //------------------------------
     bool edge = false;
     bool edgex = false;
     bool edgey = false;
 
 
     //Variables for SiPixelTemplate pixel hit error output
     float sigmay, sigmax, sy1, sy2, sx1, sx2;  
     templ.temperrors( ID, cotalpha, cotbeta, nqbin,       // inputs
            sigmay, sigmax, sy1, sy2, sx1, sx2 );      // outputs
 
     // 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;
         }
     }
     
     theError = LocalError( theErrorX*theErrorX, 0., theErrorY*theErrorY);
 
     
     unsigned int retry = 0;
     do
     {
         const SimpleHistogramGenerator* xgen = new SimpleHistogramGenerator( (TH1F*) theMergedPixelResolutionXFile-> Get("th1x")); 
         const SimpleHistogramGenerator* ygen = new SimpleHistogramGenerator( (TH1F*) theMergedPixelResolutionYFile-> Get("th1y")); 
 
         // Generate the position (x,y of the rec hit).
         theShiftInX = xgen->generate(random);
         theShiftInY = ygen->generate(random);
 
     // Now multiply by the sign of the cotangent of appropriate angle
     theShiftInX *= signOfCotalpha;
     theShiftInY *= signOfCotbeta;
 
         theShiftInZ = 0.0; // set at the centre of the active area
 
         thePosition =
       Local3DPoint( lpx + theShiftInX , 
             lpy + theShiftInY , 
             lpz + theShiftInZ );
         
         retry++;
         if (retry > 10)
         {
             // If we tried to generate thePosition, and it's out of the bounds
             // for 10 times, then take and return the simHit's location.
             thePosition = Local3DPoint(lpx,lpy,lpz);
             break;
         }
     } while(fabs(thePosition.x()) > boundX  || fabs(thePosition.y()) > boundY);
 
     FastSingleTrackerRecHit recHit(
         thePosition,
         theError,
         *detUnit,
         fastTrackerRecHitType::siPixel
     );
     return recHit;
 }