#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_ = nullptr

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 46 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, and theRegularPixelResolutions.

     : 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 128 of file PixelTemplateSmearerBase.cc.

128 {}

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 136 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 925 of file PixelTemplateSmearerBase.cc.

References PSimHit::localPosition(), log, PSimHit::momentumAtEntry(), funct::pow(), 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 171 of file PixelTemplateSmearerBase.cc.

References Exception, TrackingRecHitAlgorithm::getRandomEngine(), TrackingRecHitAlgorithm::getTrackerGeometry(), PixelTemplateSmearerBase::MergeGroup::group, watchdog::group, hitsMerge(), mps_fire::i, TrackerGeometry::idToDetUnit(), dqmiolumiharvest::j, isotrackApplyRegressor::k, Bounds::length(), isotrackNtupler::merged, mergeHitsOn, nHits, 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 680 of file PixelTemplateSmearerBase.cc.

References smearMergeGroup().

Referenced by process().

                                                                                                                        {
   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 663 of file PixelTemplateSmearerBase.cc.

References smearHit().

Referenced by process().

                                                      {
   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 314 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 698 of file PixelTemplateSmearerBase.cc.

References BXSIZE, BYSIZE, RandomEngineAndDistribution::flatShoot(), SimpleHistogramGenerator::generate(), GeomDet::geographicalId(), SiPixelTemplateDBObject::getTemplateID(), PixelTemplateSmearerBase::MergeGroup::group, mps_fire::i, SiPixelTemplate::interpolate(), PSimHit::localPosition(), microntocm, PSimHit::momentumAtEntry(), pixelTemplateDBObject_, SiPixelTemplate::qbin_dist(), 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;
 }