Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-09-07 04:36:19

0001 
0002 /** PixelTemplateSmearerBase.cc
0003  * ---------------------------------------------------------------------
0004  * Base class for FastSim plugins to simulate all simHits on one DetUnit.
0005  * 
0006  * Petar Maksimovic (JHU), based the code by 
0007  * Guofan Hu (JHU) from SiPixelGaussianSmearingRecHitConverterAlgorithm.cc
0008  * Alice Sady (JHU): new pixel resolutions (2015) and hit merging code.
0009  * ---------------------------------------------------------------------
0010  */
0011 
0012 // SiPixel Gaussian Smearing
0013 #include "FastSimulation/TrackingRecHitProducer/interface/PixelTemplateSmearerBase.h"
0014 #include "FastSimulation/TrackingRecHitProducer/interface/TrackingRecHitAlgorithmFactory.h"
0015 #include "FastSimulation/TrackingRecHitProducer/interface/TrackingRecHitProduct.h"
0016 #include "FastSimulation/TrackingRecHitProducer/interface/PixelResolutionHistograms.h"
0017 
0018 // Pixel related stuff
0019 #include "CondFormats/SiPixelObjects/interface/SiPixelTemplateDBObject.h"
0020 
0021 // Geometry
0022 /// #include "Geometry/CommonDetUnit/interface/GeomDetUnit.h"  // Keep... needed if we backport to CMSSW_9
0023 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0024 #include "Geometry/TrackerGeometryBuilder/interface/RectangularPixelTopology.h"
0025 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0026 #include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementPoint.h"
0027 
0028 // Famos
0029 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0030 #include "FastSimulation/Utilities/interface/SimpleHistogramGenerator.h"
0031 
0032 // Framework (includes ESHandle<>)
0033 #include "FWCore/Framework/interface/EventSetup.h"
0034 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0035 #include "CalibTracker/Records/interface/SiPixelTemplateDBObjectESProducerRcd.h"
0036 
0037 // ROOT
0038 #include <TFile.h>
0039 #include <TH1F.h>
0040 #include <TH2F.h>
0041 
0042 using namespace std;
0043 
0044 const double microntocm = 0.0001;
0045 
0046 PixelTemplateSmearerBase::PixelTemplateSmearerBase(const std::string& name,
0047                                                    const edm::ParameterSet& config,
0048                                                    edm::ConsumesCollector& consumesCollector)
0049     : TrackingRecHitAlgorithm(name, config, consumesCollector) {
0050   //--- Basic stuff
0051   mergeHitsOn = config.getParameter<bool>("MergeHitsOn");
0052   isBarrel = config.getParameter<bool>("isBarrel");
0053   int detType = (isBarrel) ? 1 : 0;  // 1 for barrel, 0 for forward  (or could we just promote bool into int...?)
0054 
0055   //--- Resolution file names.
0056   theBigPixelResolutionFileName = config.getParameter<string>("BigPixelResolutionFile");
0057   theBigPixelResolutionFileName = edm::FileInPath(theBigPixelResolutionFileName).fullPath();
0058 
0059   theEdgePixelResolutionFileName = config.getParameter<string>("EdgePixelResolutionFile");
0060   theEdgePixelResolutionFileName = edm::FileInPath(theEdgePixelResolutionFileName).fullPath();
0061 
0062   theRegularPixelResolutionFileName = config.getParameter<string>("RegularPixelResolutionFile");
0063   theRegularPixelResolutionFileName = edm::FileInPath(theRegularPixelResolutionFileName).fullPath();
0064 
0065   //--- Create the resolution histogram objects, which will load the histograms
0066   //    and initialize random number generators.
0067   //
0068   int status = 0;
0069   theRegularPixelResolutions = std::make_shared<PixelResolutionHistograms>(theRegularPixelResolutionFileName, "");
0070   if ((status = theRegularPixelResolutions->status()) != 0) {
0071     throw cms::Exception("PixelTemplateSmearerBase:")
0072         << " constructing PixelResolutionHistograms file " << theRegularPixelResolutionFileName
0073         << " failed with status = " << status << std::endl;
0074   }
0075 
0076   theBigPixelResolutions = std::make_shared<PixelResolutionHistograms>(
0077       theBigPixelResolutionFileName, "", detType, (!isBarrel), false, true);  // can miss qBin
0078   if ((status = theBigPixelResolutions->status()) != 0) {
0079     throw cms::Exception("PixelTemplateSmearerBase:")
0080         << " constructing PixelResolutionHistograms file " << theBigPixelResolutionFileName
0081         << " failed with status = " << status << std::endl;
0082   }
0083 
0084   theEdgePixelResolutions = std::make_shared<PixelResolutionHistograms>(
0085       theEdgePixelResolutionFileName, "", detType, false, true, true);  // can miss both single & qBin
0086   if ((status = theEdgePixelResolutions->status()) != 0) {
0087     throw cms::Exception("PixelTemplateSmearerBase:")
0088         << " constructing PixelResolutionHistograms file " << theEdgePixelResolutionFileName
0089         << " failed with status = " << status << std::endl;
0090   }
0091 
0092   //--- Merging info.
0093   theMergingProbabilityFileName = config.getParameter<string>("MergingProbabilityFile");
0094   theMergingProbabilityFileName = edm::FileInPath(theMergingProbabilityFileName).fullPath();
0095   theMergingProbabilityFile = std::make_unique<TFile>(theMergingProbabilityFileName.c_str(), "READ");
0096 
0097   theMergedPixelResolutionXFileName = config.getParameter<string>("MergedPixelResolutionXFile");
0098   theMergedPixelResolutionXFileName = edm::FileInPath(theMergedPixelResolutionXFileName).fullPath();
0099   theMergedPixelResolutionXFile = std::make_unique<TFile>(theMergedPixelResolutionXFileName.c_str(), "READ");
0100 
0101   theMergedPixelResolutionYFileName = config.getParameter<string>("MergedPixelResolutionYFile");
0102   theMergedPixelResolutionYFileName = edm::FileInPath(theMergedPixelResolutionYFileName).fullPath();
0103   theMergedPixelResolutionYFile = std::make_unique<TFile>(theMergedPixelResolutionYFileName.c_str(), "READ");
0104 
0105   // const SiPixelTemplateDBObject & dbobject;
0106   // const SiPixelTemplateDBObject dbobject;        // dummy, just to make it compile &&&
0107 
0108   //--- Load the templates.
0109   if (config.exists("templateId")) {
0110     //--- Load template with ID=templateId from a local ascii file.
0111     templateId = config.getParameter<int>("templateId");
0112     if (templateId > 0) {
0113       if (!SiPixelTemplate::pushfile(templateId, thePixelTemp_)) {
0114         throw cms::Exception("PixelTemplateSmearerBase:")
0115             << "SiPixel Template " << templateId << " Not Loaded Correctly!" << std::endl;
0116       }
0117     }
0118   }
0119 
0120   //--- Else... The templates will be loaded from the DB...
0121   //    (They are needed for data and full sim MC, so in a production FastSim
0122   //    run, everything should already be in the DB.)
0123   //
0124   //    But note that we can do it only at the beginning of the
0125   //    event.  So nothing happens now.
0126 }
0127 
0128 PixelTemplateSmearerBase::~PixelTemplateSmearerBase() {}
0129 
0130 //-------------------------------------------------------------------------------
0131 //   beginRun(); the templates are loaded in TrackingRecHitProducer, and unpacked
0132 //   into the template store.  We get their references here, and use them.  However,
0133 //   if we are loading a dedicated template ID from an ascii file just for this
0134 //   rechit smearing algorithm, then we use our own template store.
0135 //-------------------------------------------------------------------------------
0136 void PixelTemplateSmearerBase::beginRun(edm::Run const& run,
0137                                         const edm::EventSetup& eventSetup,
0138                                         const SiPixelTemplateDBObject* pixelTemplateDBObjectPtr,
0139                                         const std::vector<SiPixelTemplateStore>& tempStoreRef) {
0140   //--- Check if we need to use the template from the DB (namely if
0141   //    id == -1).  Otherwise the template has already been loaded from
0142   //    the ascii file in constructor, and thePixelTempRef wakes up
0143   //    pointing to thePixelTemp_, so then we use our own store.
0144   //
0145   if (templateId == -1) {
0146     thePixelTempRef = &tempStoreRef;                    // we use the store from TrackingRecHitProducer
0147     pixelTemplateDBObject_ = pixelTemplateDBObjectPtr;  // needed for template<-->DetId map.
0148   }
0149 
0150   //--- Commented code below (the DB interface) should say here, in case we need it:
0151   // edm::ESHandle<SiPixelTemplateDBObject> templateDBobject;
0152   // eventSetup.get<SiPixelTemplateDBObjectESProducerRcd>().get(templateDBobject);
0153   // pixelTemplateDBObject_ = templateDBobject.product();
0154 
0155   // //--- Now that we have the DB object, load the correct templates from the DB.
0156   // //    (They are needed for data and full sim MC, so in a production FastSim
0157   // //    run, everything should already be in the DB.)
0158   // if ( !SiPixelTemplate::pushfile( *pixelTemplateDBObject_ , thePixelTemp_) ) {
0159   //   throw cms::Exception("PixelTemplateSmearerPlugin:")
0160   //    <<"SiPixel Template " << templateId << " Not Loaded Correctly!"<<endl;
0161   // }
0162 }
0163 
0164 //-------------------------------------------------------------------------------
0165 //   Simulate one DetUnit:
0166 //   1. figure out where the hits are
0167 //   2. figure out which hits merge; merge them into "merge groups"
0168 //   3. smear all individual (unmerged hits)
0169 //   4. smear all merge groups.
0170 //-------------------------------------------------------------------------------
0171 TrackingRecHitProductPtr PixelTemplateSmearerBase::process(TrackingRecHitProductPtr product) const {
0172   std::vector<std::pair<unsigned int, const PSimHit*>>& simHitIdPairs = product->getSimHitIdPairs();
0173   std::vector<const PSimHit*> simHits(simHitIdPairs.size());
0174   for (unsigned int ihit = 0; ihit < simHitIdPairs.size(); ++ihit) {
0175     simHits[ihit] = simHitIdPairs[ihit].second;
0176   }
0177 
0178   RandomEngineAndDistribution const& randomEngine = getRandomEngine();
0179 
0180   const GeomDet* geomDet = getTrackerGeometry().idToDetUnit(product->getDetId());
0181   const PixelGeomDetUnit* pixelGeomDet = dynamic_cast<const PixelGeomDetUnit*>(geomDet);
0182   if (pixelGeomDet == nullptr) {
0183     throw cms::Exception("FastSimulation/TrackingRecHitProducer")
0184         << "The GeomDetUnit is not a PixelGeomDetUnit.  This should never happen!";
0185   }
0186   const BoundPlane& theDetPlane = pixelGeomDet->surface();
0187   const Bounds& theBounds = theDetPlane.bounds();
0188   const double boundX = theBounds.width() / 2.;
0189   const double boundY = theBounds.length() / 2.;
0190 
0191   std::vector<TrackingRecHitProduct::SimHitIdPair> listOfUnmergedHits;
0192   std::vector<MergeGroup*> listOfMergeGroups;
0193   int nHits = simHits.size();
0194 
0195   // fixed size array, 0 if hit is unmerged
0196   MergeGroup* mergeGroupByHit[nHits];
0197 
0198   if (nHits == 0) {
0199     return product;
0200   } else if (nHits == 1) {
0201     listOfUnmergedHits.push_back(simHitIdPairs[0]);
0202   } else {
0203     if (mergeHitsOn) {
0204       for (int i = 0; i < nHits; ++i) {
0205         //initialize this cell to a NULL pointer here
0206         mergeGroupByHit[i] = nullptr;
0207       }
0208       for (int i = 0; i < nHits - 1; ++i) {
0209         for (int j = i + 1; j < nHits; ++j) {
0210           //--- Calculate the distance between hits i and j:
0211           bool merged = hitsMerge(*simHitIdPairs[i].second, *simHitIdPairs[j].second);
0212 
0213           if (merged) {
0214             // First, check if the other guy (j) is in some merge group already
0215             if (mergeGroupByHit[j] != nullptr) {
0216               if (mergeGroupByHit[i] == nullptr) {
0217                 mergeGroupByHit[i] = mergeGroupByHit[j];
0218                 mergeGroupByHit[i]->group.push_back(simHitIdPairs[i]);
0219                 mergeGroupByHit[i]->smearIt = true;
0220               } else {
0221                 if (mergeGroupByHit[i] != mergeGroupByHit[j]) {
0222                   for (auto hit_it = mergeGroupByHit[j]->group.begin(); hit_it != mergeGroupByHit[j]->group.end();
0223                        ++hit_it) {
0224                     mergeGroupByHit[i]->group.push_back(*hit_it);
0225                     mergeGroupByHit[i]->smearIt = true;
0226                   }
0227 
0228                   // Step 2: iterate over all hits, replace mgbh[j] by mgbh[i] (so that nobody points to i)
0229                   MergeGroup* mgbhj = mergeGroupByHit[j];
0230                   for (int k = 0; k < nHits; ++k) {
0231                     if (mgbhj == mergeGroupByHit[k]) {
0232                       // Hit k also uses the same merge group, tell them to switch to mgbh[i]
0233                       mergeGroupByHit[k] = mergeGroupByHit[i];
0234                     }
0235                   }
0236                   mgbhj->smearIt = false;
0237                   mergeGroupByHit[i]->smearIt = true;
0238 
0239                   //  Step 3 would have been to delete mgbh[j]... however, we'll do that at the end anyway.
0240                   //  The key was to prevent mgbh[j] from being accessed further, and we have done that,
0241                   //  since now no mergeGroupByHit[] points to mgbhj any more.  Note that the above loop
0242                   //  also set mergeGroupByHit[i] = mergeGroupByHit[j], too.
0243                 }
0244               }
0245             } else {
0246               // j is not merged.  Check if i is merged with another hit yet.
0247               //
0248               if (mergeGroupByHit[i] == nullptr) {
0249                 // This is the first time we realized i is merged with any
0250                 // other hit.  Create a new merge group for i and j
0251                 mergeGroupByHit[i] = new MergeGroup();
0252                 listOfMergeGroups.push_back(mergeGroupByHit[i]);  // keep track of it
0253                 //
0254                 // Add hit i as the first to its own merge group
0255                 // (simHits[i] is a const pointer to PSimHit).
0256                 mergeGroupByHit[i]->group.push_back(simHitIdPairs[i]);
0257                 mergeGroupByHit[i]->smearIt = true;
0258               }
0259               //--- Add hit j as well
0260               mergeGroupByHit[i]->group.push_back(simHitIdPairs[j]);
0261               mergeGroupByHit[i]->smearIt = true;
0262 
0263               mergeGroupByHit[j] = mergeGroupByHit[i];
0264 
0265             }  // --- end of else if ( j has merge group )
0266 
0267           }  //--- end of if (merged)
0268 
0269         }  //--- end of loop over j
0270 
0271         //--- At this point, there are two possibilities.  Either hit i
0272         //    was already chosen to be merged with some hit prior to it,
0273         //    or the loop over j found another merged hit.  In either
0274         //    case, if mergeGroupByHit[i] is empty, then the hit is
0275         //    unmerged.
0276         //
0277         if (mergeGroupByHit[i] == nullptr) {
0278           //--- Keep track of it.
0279           listOfUnmergedHits.push_back(simHitIdPairs[i]);
0280         }
0281       }  //--- end of loop over i
0282     }  // --- end of if (mergeHitsOn)
0283     else {
0284       // Now we've turned off hit merging, so all hits should be pushed
0285       // back to listOfUnmergedHits
0286       for (int i = 0; i < nHits; ++i) {
0287         listOfUnmergedHits.push_back(simHitIdPairs[i]);
0288       }
0289     }
0290   }  // --- end of if (nHits == 1) else {...}
0291 
0292   //--- We now have two lists: a list of hits that are unmerged, and
0293   //    the list of merge groups.  Process each separately.
0294   //
0295   product = processUnmergedHits(listOfUnmergedHits, product, pixelGeomDet, boundX, boundY, &randomEngine);
0296 
0297   product = processMergeGroups(listOfMergeGroups, product, pixelGeomDet, boundX, boundY, &randomEngine);
0298 
0299   //--- We're done with this det unit, and ought to clean up used
0300   //    memory.  We don't own the PSimHits, and the vector of
0301   //    listOfUnmergedHits simply goes out of scope.  However, we
0302   //    created the MergeGroups and thus we need to get rid of them.
0303   //
0304   for (auto mg_it = listOfMergeGroups.begin(); mg_it != listOfMergeGroups.end(); ++mg_it) {
0305     delete *mg_it;  // each MergeGroup is deleted; its ptrs to PSimHits we do not own...
0306   }
0307 
0308   return product;
0309 }
0310 
0311 //------------------------------------------------------------------------------
0312 //   Smear one hit.  The main action is in here.
0313 //------------------------------------------------------------------------------
0314 FastSingleTrackerRecHit PixelTemplateSmearerBase::smearHit(const PSimHit& simHit,
0315                                                            const PixelGeomDetUnit* detUnit,
0316                                                            const double boundX,
0317                                                            const double boundY,
0318                                                            RandomEngineAndDistribution const* random) const {
0319   //--- At the beginning the position is the Local Point in the local pixel module reference frame
0320   //    same code as in PixelCPEBase
0321   //
0322   LocalVector localDir = simHit.momentumAtEntry();  // don't need .unit(), we will take the ratio
0323   float locx = localDir.x();
0324   float locy = localDir.y();
0325   float locz = localDir.z();
0326 
0327   //--- cotangent of local angles \alpha and \beta.
0328   //    alpha: angle with respect to local x axis in local (x,z) plane
0329   //    beta: angle with respect to local y axis in local (y,z) plane
0330   //
0331   float cotalpha = locx / locz;
0332   float cotbeta = locy / locz;
0333 
0334   //--- Save the original signs of cot\alpha and cot\beta
0335   int signOfCotalpha = (cotalpha < 0) ? -1 : 1;  // sign(cotalpha);
0336   int signOfCotbeta = (cotbeta < 0) ? -1 : 1;    // sign(cotbeta);
0337   //
0338   //--- Use absolute values to find the templates from the list
0339   cotalpha *= signOfCotalpha;  // = abs(cotalpha)
0340   cotbeta *= signOfCotbeta;    // = abs(cotbeta)
0341 
0342   LogDebug("SmearHit") << "LocalVector=" << locx << "," << locy << "," << locz << "   momentum=" << localDir.mag()
0343                        << "   cotalpha=" << cotalpha << ",  cotbeta=" << cotbeta;
0344 
0345   const PixelTopology* theSpecificTopology = &(detUnit->specificType().specificTopology());
0346   const RectangularPixelTopology* rectPixelTopology = static_cast<const RectangularPixelTopology*>(theSpecificTopology);
0347 
0348   const int nrows = theSpecificTopology->nrows();
0349   const int ncolumns = theSpecificTopology->ncolumns();
0350 
0351   const Local3DPoint lp = simHit.localPosition();
0352   //Transform local position to measurement position
0353   const MeasurementPoint mp = rectPixelTopology->measurementPosition(lp);
0354   float mpy = mp.y();
0355   float mpx = mp.x();
0356   //Get the center of the struck pixel in measurement position
0357   float pixelCenterY = 0.5 + (int)mpy;
0358   float pixelCenterX = 0.5 + (int)mpx;
0359 
0360   const MeasurementPoint mpCenter(pixelCenterX, pixelCenterY);
0361   //Transform the center of the struck pixel back into local position
0362   const Local3DPoint lpCenter = rectPixelTopology->localPosition(mpCenter);
0363 
0364   //Get the relative position of struck point to the center of the struck pixel
0365   float xtrk = lp.x() - lpCenter.x();
0366   float ytrk = lp.y() - lpCenter.y();
0367   //Pixel Y, X pitch
0368   const float ysize = {0.015}, xsize = {0.01};
0369   //Variables for SiPixelTemplate input, see SiPixelTemplate reco
0370   float yhit = 20. + 8. * (ytrk / ysize);
0371   float xhit = 20. + 8. * (xtrk / xsize);
0372   int ybin = (int)yhit;
0373   int xbin = (int)xhit;
0374   float yfrac = yhit - (float)ybin;
0375   float xfrac = xhit - (float)xbin;
0376   //Protect againt ybin, xbin being outside of range [0-39]  // &&& Why limit of 39?
0377   if (ybin < 0)
0378     ybin = 0;
0379   if (ybin > 39)
0380     ybin = 39;
0381   if (xbin < 0)
0382     xbin = 0;
0383   if (xbin > 39)
0384     xbin = 39;
0385 
0386   int ID = templateId;
0387   if (templateId == -1) {
0388     // We have loaded the whole template set from the DB,
0389     // so ask the DB object to find us the right one.
0390     ID = pixelTemplateDBObject_->getTemplateID(detUnit->geographicalId());  // need uint32_t detid
0391     //                  theDetParam.theDet->geographicalId());
0392   }
0393 
0394   //--- Make the template object
0395   SiPixelTemplate templ(*thePixelTempRef);
0396 
0397   //--- Produce the template that corresponds to our local angles.
0398   templ.interpolate(ID, cotalpha, cotbeta);
0399 
0400   //Variables for SiPixelTemplate output
0401   //qBin -- normalized pixel charge deposition
0402   float qbin_frac[4];
0403   //Single pixel cluster projection possibility
0404   float ny1_frac, ny2_frac, nx1_frac, nx2_frac;
0405   bool singlex = false, singley = false;
0406   templ.qbin_dist(ID, cotalpha, cotbeta, qbin_frac, ny1_frac, ny2_frac, nx1_frac, nx2_frac);
0407   int nqbin;
0408 
0409   double xsizeProbability = random->flatShoot();
0410   double ysizeProbability = random->flatShoot();
0411   bool hitbigx = rectPixelTopology->isItBigPixelInX((int)mpx);  // pixel we hit in x
0412   bool hitbigy = rectPixelTopology->isItBigPixelInY((int)mpy);  // pixel we hit in y
0413 
0414   if (hitbigx)
0415     if (xsizeProbability < nx2_frac)
0416       singlex = true;
0417     else
0418       singlex = false;
0419   else if (xsizeProbability < nx1_frac)
0420     singlex = true;
0421   else
0422     singlex = false;
0423 
0424   if (hitbigy)
0425     if (ysizeProbability < ny2_frac)
0426       singley = true;
0427     else
0428       singley = false;
0429   else if (ysizeProbability < ny1_frac)
0430     singley = true;
0431   else
0432     singley = false;
0433 
0434   // random multiplicity for alpha and beta
0435   double qbinProbability = random->flatShoot();
0436   for (int i = 0; i < 4; ++i) {
0437     nqbin = i;
0438     if (qbinProbability < qbin_frac[i]) {
0439       break;
0440     }
0441   }
0442 
0443   //Store interpolated pixel cluster profile
0444   //BYSIZE, BXSIZE, const definition from SiPixelTemplate
0445   float ytempl[41][BYSIZE] = {{0}}, xtempl[41][BXSIZE] = {{0}};
0446   templ.ytemp(0, 40, ytempl);
0447   templ.xtemp(0, 40, xtempl);
0448 
0449   std::vector<double> ytemp(BYSIZE);
0450   for (int i = 0; i < BYSIZE; ++i) {
0451     ytemp[i] = (1. - yfrac) * ytempl[ybin][i] + yfrac * ytempl[ybin + 1][i];
0452   }
0453 
0454   std::vector<double> xtemp(BXSIZE);
0455   for (int i = 0; i < BXSIZE; ++i) {
0456     xtemp[i] = (1. - xfrac) * xtempl[xbin][i] + xfrac * xtempl[xbin + 1][i];
0457   }
0458 
0459   //Pixel readout threshold
0460   const float qThreshold = templ.s50() * 2.0;
0461 
0462   //Cut away pixels below readout threshold
0463   //For cluster lengths calculation
0464   int offsetX1 = 0, offsetX2 = 0, offsetY1 = 0, offsetY2 = 0;
0465   int firstY, lastY, firstX, lastX;
0466   for (firstY = 0; firstY < BYSIZE; ++firstY) {
0467     bool yCluster = ytemp[firstY] > qThreshold;
0468     if (yCluster) {
0469       offsetY1 = BHY - firstY;
0470       break;
0471     }
0472   }
0473   for (lastY = firstY; lastY < BYSIZE; ++lastY) {
0474     bool yCluster = ytemp[lastY] > qThreshold;
0475     if (!yCluster) {
0476       lastY = lastY - 1;
0477       offsetY2 = lastY - BHY;
0478       break;
0479     }
0480   }
0481 
0482   for (firstX = 0; firstX < BXSIZE; ++firstX) {
0483     bool xCluster = xtemp[firstX] > qThreshold;
0484     if (xCluster) {
0485       offsetX1 = BHX - firstX;
0486       break;
0487     }
0488   }
0489   for (lastX = firstX; lastX < BXSIZE; ++lastX) {
0490     bool xCluster = xtemp[lastX] > qThreshold;
0491     if (!xCluster) {
0492       lastX = lastX - 1;
0493       offsetX2 = lastX - BHX;
0494       break;
0495     }
0496   }
0497 
0498   //--- Prepare to return results
0499   Local3DPoint thePosition;
0500   double theShiftInX;
0501   double theShiftInY;
0502   double theShiftInZ;
0503   LocalError theError;
0504   double theErrorX;
0505   double theErrorY;
0506 
0507   //------------------------------
0508   //  Check if the cluster is near an edge.  If it protrudes
0509   //  outside the edge of the sensor, the truncate it and it will
0510   //  get significantly messed up.
0511   //------------------------------
0512   bool edge, edgex, edgey;
0513   //  bool bigx, bigy;
0514 
0515   int firstPixelInX = (int)mpx - offsetX1;
0516   int firstPixelInY = (int)mpy - offsetY1;
0517   int lastPixelInX = (int)mpx + offsetX2;
0518   int lastPixelInY = (int)mpy + offsetY2;
0519   firstPixelInX = (firstPixelInX >= 0) ? firstPixelInX : 0;
0520   firstPixelInY = (firstPixelInY >= 0) ? firstPixelInY : 0;
0521   lastPixelInX = (lastPixelInX < nrows) ? lastPixelInX : nrows - 1;
0522   lastPixelInY = (lastPixelInY < ncolumns) ? lastPixelInY : ncolumns - 1;
0523 
0524   edgex = rectPixelTopology->isItEdgePixelInX(firstPixelInX) || rectPixelTopology->isItEdgePixelInX(lastPixelInX);
0525   edgey = rectPixelTopology->isItEdgePixelInY(firstPixelInY) || rectPixelTopology->isItEdgePixelInY(lastPixelInY);
0526   edge = edgex || edgey;
0527 
0528   //  bigx = rectPixelTopology->isItBigPixelInX( firstPixelInX ) || rectPixelTopology->isItBigPixelInX( lastPixelInX );
0529   //  bigy = rectPixelTopology->isItBigPixelInY( firstPixelInY ) || rectPixelTopology->isItBigPixelInY( lastPixelInY );
0530   bool hasBigPixelInX = rectPixelTopology->containsBigPixelInX(firstPixelInX, lastPixelInX);
0531   bool hasBigPixelInY = rectPixelTopology->containsBigPixelInY(firstPixelInY, lastPixelInY);
0532 
0533   //Variables for SiPixelTemplate pixel hit error output
0534   float sigmay, sigmax, sy1, sy2, sx1, sx2;
0535   templ.temperrors(ID,
0536                    cotalpha,
0537                    cotbeta,
0538                    nqbin,  // inputs
0539                    sigmay,
0540                    sigmax,
0541                    sy1,
0542                    sy2,
0543                    sx1,
0544                    sx2  // outputs
0545   );
0546 
0547   if (edge) {
0548     if (edgex && !edgey) {
0549       theErrorX = 23.0 * microntocm;
0550       theErrorY = 39.0 * microntocm;
0551     } else if (!edgex && edgey) {
0552       theErrorX = 24.0 * microntocm;
0553       theErrorY = 96.0 * microntocm;
0554     } else {
0555       theErrorX = 31.0 * microntocm;
0556       theErrorY = 90.0 * microntocm;
0557     }
0558   } else {
0559     if (singlex) {
0560       if (hitbigx) {
0561         theErrorX = sx2 * microntocm;
0562       } else {
0563         theErrorX = sx1 * microntocm;
0564       }
0565     } else {
0566       theErrorX = sigmax * microntocm;
0567     }
0568     if (singley) {
0569       if (hitbigy) {
0570         theErrorY = sy2 * microntocm;
0571       } else {
0572         theErrorY = sy1 * microntocm;
0573       }
0574     } else {
0575       theErrorY = sigmay * microntocm;
0576     }
0577   }
0578 
0579   //add misalignment error
0580   const TrackerGeomDet* misalignmentDetUnit = getMisalignedGeometry().idToDet(detUnit->geographicalId());
0581   const LocalError& alignmentError = misalignmentDetUnit->localAlignmentError();
0582   if (alignmentError.valid()) {
0583     theError = LocalError(
0584         theErrorX * theErrorX + alignmentError.xx(), alignmentError.xy(), theErrorY * theErrorY + alignmentError.yy());
0585   } else {
0586     theError = LocalError(theErrorX * theErrorX, 0.0, theErrorY * theErrorY);
0587   }
0588 
0589   // Local Error is 2D: (xx,xy,yy), square of sigma in first an third position
0590   // as for resolution matrix
0591 
0592   //--- Next, we need to generate the smeared position.  First we need to figure
0593   //    out which kind of histograms we are supposed to use for this particular hit.
0594   //    These are pointers to the set of histograms used to generate the rec hit
0595   //    positions.  (We need to handle X and Y separately.)
0596   shared_ptr<PixelResolutionHistograms> resHistsX = nullptr;
0597   shared_ptr<PixelResolutionHistograms> resHistsY = nullptr;
0598 
0599   if (edge) {
0600     resHistsX = resHistsY = theEdgePixelResolutions;
0601     singlex = singley = false;  // no single resolutions for Edge
0602   } else {
0603     //--- Decide resolution histogram set for X
0604     if ((singlex && hitbigx) || (isBarrel && hasBigPixelInX)) {
0605       resHistsX = theBigPixelResolutions;
0606     } else {
0607       resHistsX = theRegularPixelResolutions;
0608     }
0609     //--- Decide resolution histogram set for Y
0610     if ((singley && hitbigy) || (isBarrel && hasBigPixelInY)) {
0611       resHistsY = theBigPixelResolutions;
0612     } else {
0613       resHistsY = theRegularPixelResolutions;
0614     }
0615   }
0616 
0617   //--- Get generators, separately for X and for Y.
0618   const SimpleHistogramGenerator* xgen = resHistsX->getGeneratorX(cotalpha, cotbeta, nqbin, singlex);
0619   const SimpleHistogramGenerator* ygen = resHistsY->getGeneratorY(cotalpha, cotbeta, nqbin, singley);
0620 
0621   //--- Check if we found a histogram.  If nullptr, then throw up.
0622   if (!xgen || !ygen) {
0623     throw cms::Exception("FastSimulation/TrackingRecHitProducer")
0624         << "Histogram (cot\alpha=" << cotalpha << ", cot\beta=" << cotbeta << ", nQbin=" << nqbin
0625         << ") was not found for PixelTemplateSmearer. Check if the smearing resolution histogram exists.";
0626   }
0627 
0628   //--- Smear the hit Position.  We do it in the do-while loop in order to
0629   //--- allow multiple tries, in case we generate a rec hit which is outside
0630   //--- of the boundaries of the sensor.
0631   unsigned int retry = 0;
0632 
0633   do {
0634     // Generate the position (x,y of the rec hit).
0635     theShiftInX = xgen->generate(random);
0636     theShiftInY = ygen->generate(random);
0637 
0638     // Now multiply by the sign of the cotangent of appropriate angle
0639     theShiftInX *= signOfCotalpha;
0640     theShiftInY *= signOfCotbeta;
0641 
0642     theShiftInZ = 0.0;  // set to the mid-plane of the sensor.
0643 
0644     thePosition = Local3DPoint(simHit.localPosition().x() + theShiftInX,
0645                                simHit.localPosition().y() + theShiftInY,
0646                                simHit.localPosition().z() + theShiftInZ);
0647     retry++;
0648     if (retry > 10) {
0649       // If we tried to generate thePosition, and it's out of the bounds
0650       // for 10 times, then take and return the simHit's location.
0651       thePosition = Local3DPoint(simHit.localPosition().x(), simHit.localPosition().y(), simHit.localPosition().z());
0652       break;
0653     }
0654   } while (fabs(thePosition.x()) > boundX || fabs(thePosition.y()) > boundY);
0655 
0656   FastSingleTrackerRecHit recHit(thePosition, theError, *detUnit, fastTrackerRecHitType::siPixel);
0657   return recHit;
0658 }
0659 
0660 //------------------------------------------------------------------------------
0661 //   Smear all umerged hits on this DetUnit
0662 //------------------------------------------------------------------------------
0663 TrackingRecHitProductPtr PixelTemplateSmearerBase::processUnmergedHits(
0664     std::vector<TrackingRecHitProduct::SimHitIdPair>& unmergedHits,
0665     TrackingRecHitProductPtr product,
0666     const PixelGeomDetUnit* detUnit,
0667     const double boundX,
0668     const double boundY,
0669     RandomEngineAndDistribution const* random) const {
0670   for (auto simHitIdPair : unmergedHits) {
0671     FastSingleTrackerRecHit recHit = smearHit(*simHitIdPair.second, detUnit, boundX, boundY, random);
0672     product->addRecHit(recHit, {simHitIdPair});
0673   }
0674   return product;
0675 }
0676 
0677 //------------------------------------------------------------------------------
0678 //   Smear all MERGED hits on this DetUnit
0679 //------------------------------------------------------------------------------
0680 TrackingRecHitProductPtr PixelTemplateSmearerBase::processMergeGroups(std::vector<MergeGroup*>& mergeGroups,
0681                                                                       TrackingRecHitProductPtr product,
0682                                                                       const PixelGeomDetUnit* detUnit,
0683                                                                       const double boundX,
0684                                                                       const double boundY,
0685                                                                       RandomEngineAndDistribution const* random) const {
0686   for (auto mg_it = mergeGroups.begin(); mg_it != mergeGroups.end(); ++mg_it) {
0687     if ((*mg_it)->smearIt) {
0688       FastSingleTrackerRecHit recHit = smearMergeGroup(*mg_it, detUnit, boundX, boundY, random);
0689       product->addRecHit(recHit, (*mg_it)->group);
0690     }
0691   }
0692   return product;
0693 }
0694 
0695 //------------------------------------------------------------------------------
0696 //   Smear all hits MERGED together.  This is called a MergeGroup.
0697 //------------------------------------------------------------------------------
0698 FastSingleTrackerRecHit PixelTemplateSmearerBase::smearMergeGroup(MergeGroup* mg,
0699                                                                   const PixelGeomDetUnit* detUnit,
0700                                                                   const double boundX,
0701                                                                   const double boundY,
0702                                                                   RandomEngineAndDistribution const* random) const {
0703   float loccx = 0;
0704   float loccy = 0;
0705   float loccz = 0;
0706   float nHit = 0;
0707   float locpx = 0;
0708   float locpy = 0;
0709   float locpz = 0;
0710 
0711   for (auto hit_it = mg->group.begin(); hit_it != mg->group.end(); ++hit_it) {
0712     const PSimHit simHit = *hit_it->second;
0713     //getting local momentum and adding all of the hits' momentums up
0714     LocalVector localDir = simHit.momentumAtEntry().unit();
0715     loccx += localDir.x();
0716     loccy += localDir.y();
0717     loccz += localDir.z();
0718     //getting local position and adding all of the hits' positions up
0719     const Local3DPoint lpos = simHit.localPosition();
0720     locpx += lpos.x();
0721     locpy += lpos.y();
0722     locpz += lpos.z();
0723     //counting how many sim hits are in the merge group
0724     nHit += 1;
0725   }
0726   //averaging the momentums by diving momentums added up/number of hits
0727   float locx = loccx / nHit;
0728   float locy = loccy / nHit;
0729   float locz = loccz / nHit;
0730 
0731   //--- cotangent of local angles \alpha and \beta.
0732   //    alpha: angle with respect to local x axis in local (x,z) plane
0733   //    beta: angle with respect to local y axis in local (y,z) plane
0734   //
0735   float cotalpha = locx / locz;
0736   float cotbeta = locy / locz;
0737 
0738   //--- Save the original signs of cot\alpha and cot\beta
0739   int signOfCotalpha = (cotalpha < 0) ? -1 : 1;  // sign(cotalpha);
0740   int signOfCotbeta = (cotbeta < 0) ? -1 : 1;    // sign(cotbeta);
0741   //
0742   //--- Use absolute values to find the templates from the list
0743   cotalpha *= signOfCotalpha;  // = abs(cotalpha)
0744   cotbeta *= signOfCotbeta;    // = abs(cotbeta)
0745 
0746   float lpx = locpx / nHit;
0747   float lpy = locpy / nHit;
0748   float lpz = locpz / nHit;
0749 
0750   //Get the relative position of struck point to the center of the struck pixel
0751   float xtrk = lpx;
0752   float ytrk = lpy;
0753   //Pixel Y, X pitch
0754   const float ysize = {0.015}, xsize = {0.01};
0755   //Variables for SiPixelTemplate input, see SiPixelTemplate reco
0756   float yhit = 20. + 8. * (ytrk / ysize);
0757   float xhit = 20. + 8. * (xtrk / xsize);
0758   int ybin = (int)yhit;
0759   int xbin = (int)xhit;
0760   float yfrac = yhit - (float)ybin;
0761   float xfrac = xhit - (float)xbin;
0762   // Protect againt ybin, xbin being outside of range [0-39]
0763   if (ybin < 0)
0764     ybin = 0;
0765   if (ybin > 39)
0766     ybin = 39;
0767   if (xbin < 0)
0768     xbin = 0;
0769   if (xbin > 39)
0770     xbin = 39;
0771 
0772   int ID = templateId;
0773   if (templateId == -1) {
0774     // We have loaded the whole template set from the DB,
0775     // so ask the DB object to find us the right one.
0776     ID = pixelTemplateDBObject_->getTemplateID(detUnit->geographicalId());  // need uint32_t detid
0777     //                  theDetParam.theDet->geographicalId());
0778   }
0779 
0780   //--- Make the template object
0781   SiPixelTemplate templ(*thePixelTempRef);
0782 
0783   //--- Produce the template that corresponds to our local angles.
0784   templ.interpolate(ID, cotalpha, cotbeta);
0785 
0786   // Variables for SiPixelTemplate output
0787   // qBin -- normalized pixel charge deposition
0788   float qbin_frac[4];
0789   // Single pixel cluster projection possibility
0790   float ny1_frac, ny2_frac, nx1_frac, nx2_frac;
0791   bool singlex = false, singley = false;
0792   templ.qbin_dist(ID, cotalpha, cotbeta, qbin_frac, ny1_frac, ny2_frac, nx1_frac, nx2_frac);
0793   int nqbin;
0794 
0795   //  double xsizeProbability = random->flatShoot();
0796   //double ysizeProbability = random->flatShoot();
0797   bool hitbigx = false;
0798   bool hitbigy = false;
0799 
0800   // random multiplicity for alpha and beta
0801 
0802   double qbinProbability = random->flatShoot();
0803   for (int i = 0; i < 4; ++i) {
0804     nqbin = i;
0805     if (qbinProbability < qbin_frac[i])
0806       break;
0807   }
0808 
0809   //Store interpolated pixel cluster profile
0810   //BYSIZE, BXSIZE, const definition from SiPixelTemplate
0811   float ytempl[41][BYSIZE] = {{0}}, xtempl[41][BXSIZE] = {{0}};
0812   templ.ytemp(0, 40, ytempl);
0813   templ.xtemp(0, 40, xtempl);
0814 
0815   std::vector<double> ytemp(BYSIZE);
0816   for (int i = 0; i < BYSIZE; ++i) {
0817     ytemp[i] = (1. - yfrac) * ytempl[ybin][i] + yfrac * ytempl[ybin + 1][i];
0818   }
0819 
0820   std::vector<double> xtemp(BXSIZE);
0821   for (int i = 0; i < BXSIZE; ++i) {
0822     xtemp[i] = (1. - xfrac) * xtempl[xbin][i] + xfrac * xtempl[xbin + 1][i];
0823   }
0824 
0825   //--- Prepare to return results
0826   Local3DPoint thePosition;
0827   double theShiftInX;
0828   double theShiftInY;
0829   double theShiftInZ;
0830   LocalError theError;
0831   double theErrorX;
0832   double theErrorY;
0833 
0834   //------------------------------
0835   //  Check if the cluster is near an edge.  If it protrudes
0836   //  outside the edge of the sensor, the truncate it and it will
0837   //  get significantly messed up.
0838   //------------------------------
0839   bool edge = false;
0840   bool edgex = false;
0841   bool edgey = false;
0842 
0843   //Variables for SiPixelTemplate pixel hit error output
0844   float sigmay, sigmax, sy1, sy2, sx1, sx2;
0845   templ.temperrors(ID,
0846                    cotalpha,
0847                    cotbeta,
0848                    nqbin,  // inputs
0849                    sigmay,
0850                    sigmax,
0851                    sy1,
0852                    sy2,
0853                    sx1,
0854                    sx2);  // outputs
0855 
0856   // define private mebers --> Errors
0857   if (edge) {
0858     if (edgex && !edgey) {
0859       theErrorX = 23.0 * microntocm;
0860       theErrorY = 39.0 * microntocm;
0861     } else if (!edgex && edgey) {
0862       theErrorX = 24.0 * microntocm;
0863       theErrorY = 96.0 * microntocm;
0864     } else {
0865       theErrorX = 31.0 * microntocm;
0866       theErrorY = 90.0 * microntocm;
0867     }
0868 
0869   } else {
0870     if (singlex) {
0871       if (hitbigx) {
0872         theErrorX = sx2 * microntocm;
0873       } else {
0874         theErrorX = sx1 * microntocm;
0875       }
0876     } else {
0877       theErrorX = sigmax * microntocm;
0878     }
0879 
0880     if (singley) {
0881       if (hitbigy) {
0882         theErrorY = sy2 * microntocm;
0883       } else {
0884         theErrorY = sy1 * microntocm;
0885       }
0886     } else {
0887       theErrorY = sigmay * microntocm;
0888     }
0889   }
0890 
0891   theError = LocalError(theErrorX * theErrorX, 0., theErrorY * theErrorY);
0892 
0893   unsigned int retry = 0;
0894   do {
0895     const SimpleHistogramGenerator* xgen =
0896         new SimpleHistogramGenerator((TH1F*)theMergedPixelResolutionXFile->Get("th1x"));
0897     const SimpleHistogramGenerator* ygen =
0898         new SimpleHistogramGenerator((TH1F*)theMergedPixelResolutionYFile->Get("th1y"));
0899 
0900     // Generate the position (x,y of the rec hit).
0901     theShiftInX = xgen->generate(random);
0902     theShiftInY = ygen->generate(random);
0903 
0904     // Now multiply by the sign of the cotangent of appropriate angle
0905     theShiftInX *= signOfCotalpha;
0906     theShiftInY *= signOfCotbeta;
0907 
0908     theShiftInZ = 0.0;  // set at the centre of the active area
0909 
0910     thePosition = Local3DPoint(lpx + theShiftInX, lpy + theShiftInY, lpz + theShiftInZ);
0911 
0912     retry++;
0913     if (retry > 10) {
0914       // If we tried to generate thePosition, and it's out of the bounds
0915       // for 10 times, then take and return the simHit's location.
0916       thePosition = Local3DPoint(lpx, lpy, lpz);
0917       break;
0918     }
0919   } while (fabs(thePosition.x()) > boundX || fabs(thePosition.y()) > boundY);
0920 
0921   FastSingleTrackerRecHit recHit(thePosition, theError, *detUnit, fastTrackerRecHitType::siPixel);
0922   return recHit;
0923 }
0924 
0925 bool PixelTemplateSmearerBase::hitsMerge(const PSimHit& simHit1, const PSimHit& simHit2) const {
0926   LocalVector localDir = simHit1.momentumAtEntry().unit();
0927   float locy1 = localDir.y();
0928   float locz1 = localDir.z();
0929   float cotbeta = locy1 / locz1;
0930   float loceta = fabs(-log((double)(-cotbeta + sqrt((double)(1. + cotbeta * cotbeta)))));
0931 
0932   const Local3DPoint lp1 = simHit1.localPosition();
0933   const Local3DPoint lp2 = simHit2.localPosition();
0934   float lpy1 = lp1.y();
0935   float lpx1 = lp1.x();
0936   float lpy2 = lp2.y();
0937   float lpx2 = lp2.x();
0938   float locdis = 10000. * sqrt(pow(lpx1 - lpx2, 2) + pow(lpy1 - lpy2, 2));
0939   TH2F* probhisto = (TH2F*)theMergingProbabilityFile->Get("h2bc");
0940   float prob =
0941       probhisto->GetBinContent(probhisto->GetXaxis()->FindFixBin(locdis), probhisto->GetYaxis()->FindFixBin(loceta));
0942   return prob > 0;
0943 }