Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-02-01 23:47:19

0001 #include <cuda_runtime.h>
0002 
0003 #include "CondFormats/SiPixelTransient/interface/SiPixelTemplate.h"
0004 #include "DataFormats/DetId/interface/DetId.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0007 #include "Geometry/TrackerGeometryBuilder/interface/RectangularPixelTopology.h"
0008 #include "Geometry/CommonTopologies/interface/SimplePixelTopology.h"
0009 #include "HeterogeneousCore/CUDAUtilities/interface/cudaCheck.h"
0010 #include "MagneticField/Engine/interface/MagneticField.h"
0011 #include "RecoLocalTracker/SiPixelRecHits/interface/PixelCPEFast.h"
0012 
0013 // Services
0014 // this is needed to get errors from templates
0015 
0016 namespace {
0017   constexpr float micronsToCm = 1.0e-4;
0018 }
0019 
0020 //-----------------------------------------------------------------------------
0021 //!  The constructor.
0022 //-----------------------------------------------------------------------------
0023 PixelCPEFast::PixelCPEFast(edm::ParameterSet const& conf,
0024                            const MagneticField* mag,
0025                            const TrackerGeometry& geom,
0026                            const TrackerTopology& ttopo,
0027                            const SiPixelLorentzAngle* lorentzAngle,
0028                            const SiPixelGenErrorDBObject* genErrorDBObject,
0029                            const SiPixelLorentzAngle* lorentzAngleWidth)
0030     : PixelCPEGenericBase(conf, mag, geom, ttopo, lorentzAngle, genErrorDBObject, lorentzAngleWidth) {
0031   // Use errors from templates or from GenError
0032   if (useErrorsFromTemplates_) {
0033     if (!SiPixelGenError::pushfile(*genErrorDBObject_, thePixelGenError_))
0034       throw cms::Exception("InvalidCalibrationLoaded")
0035           << "ERROR: GenErrors not filled correctly. Check the sqlite file. Using SiPixelTemplateDBObject version "
0036           << (*genErrorDBObject_).version();
0037   }
0038 
0039   isPhase2_ = conf.getParameter<bool>("isPhase2");
0040 
0041   fillParamsForGpu();
0042 
0043   cpuData_ = {
0044       &commonParamsGPU_,
0045       detParamsGPU_.data(),
0046       &layerGeometry_,
0047       &averageGeometry_,
0048   };
0049 }
0050 
0051 const pixelCPEforGPU::ParamsOnGPU* PixelCPEFast::getGPUProductAsync(cudaStream_t cudaStream) const {
0052   const auto& data = gpuData_.dataForCurrentDeviceAsync(cudaStream, [this](GPUData& data, cudaStream_t stream) {
0053     // and now copy to device...
0054 
0055     cudaCheck(cudaMalloc((void**)&data.paramsOnGPU_h.m_commonParams, sizeof(pixelCPEforGPU::CommonParams)));
0056     cudaCheck(cudaMalloc((void**)&data.paramsOnGPU_h.m_detParams,
0057                          this->detParamsGPU_.size() * sizeof(pixelCPEforGPU::DetParams)));
0058     cudaCheck(cudaMalloc((void**)&data.paramsOnGPU_h.m_averageGeometry, sizeof(pixelCPEforGPU::AverageGeometry)));
0059     cudaCheck(cudaMalloc((void**)&data.paramsOnGPU_h.m_layerGeometry, sizeof(pixelCPEforGPU::LayerGeometry)));
0060     cudaCheck(cudaMalloc((void**)&data.paramsOnGPU_d, sizeof(pixelCPEforGPU::ParamsOnGPU)));
0061     cudaCheck(cudaMemcpyAsync(
0062         data.paramsOnGPU_d, &data.paramsOnGPU_h, sizeof(pixelCPEforGPU::ParamsOnGPU), cudaMemcpyDefault, stream));
0063     cudaCheck(cudaMemcpyAsync((void*)data.paramsOnGPU_h.m_commonParams,
0064                               &this->commonParamsGPU_,
0065                               sizeof(pixelCPEforGPU::CommonParams),
0066                               cudaMemcpyDefault,
0067                               stream));
0068     cudaCheck(cudaMemcpyAsync((void*)data.paramsOnGPU_h.m_averageGeometry,
0069                               &this->averageGeometry_,
0070                               sizeof(pixelCPEforGPU::AverageGeometry),
0071                               cudaMemcpyDefault,
0072                               stream));
0073     cudaCheck(cudaMemcpyAsync((void*)data.paramsOnGPU_h.m_layerGeometry,
0074                               &this->layerGeometry_,
0075                               sizeof(pixelCPEforGPU::LayerGeometry),
0076                               cudaMemcpyDefault,
0077                               stream));
0078     cudaCheck(cudaMemcpyAsync((void*)data.paramsOnGPU_h.m_detParams,
0079                               this->detParamsGPU_.data(),
0080                               this->detParamsGPU_.size() * sizeof(pixelCPEforGPU::DetParams),
0081                               cudaMemcpyDefault,
0082                               stream));
0083   });
0084   return data.paramsOnGPU_d;
0085 }
0086 
0087 void PixelCPEFast::fillParamsForGpu() {
0088   //
0089   // this code executes only once per job, computation inefficiency is not an issue
0090   // many code blocks are repeated: better keep the computation local and self oconsistent as blocks may in future move around, be deleted ...
0091   // It is valid only for Phase1 and the version of GenError in DB used in late 2018 and in 2021
0092 
0093   commonParamsGPU_.theThicknessB = m_DetParams.front().theThickness;
0094   commonParamsGPU_.theThicknessE = m_DetParams.back().theThickness;
0095   commonParamsGPU_.thePitchX = m_DetParams[0].thePitchX;
0096   commonParamsGPU_.thePitchY = m_DetParams[0].thePitchY;
0097 
0098   commonParamsGPU_.numberOfLaddersInBarrel =
0099       isPhase2_ ? phase2PixelTopology::numberOfLaddersInBarrel : phase1PixelTopology::numberOfLaddersInBarrel;
0100   commonParamsGPU_.isPhase2 = isPhase2_;
0101 
0102   LogDebug("PixelCPEFast") << "pitch & thickness " << commonParamsGPU_.thePitchX << ' ' << commonParamsGPU_.thePitchY
0103                            << "  " << commonParamsGPU_.theThicknessB << ' ' << commonParamsGPU_.theThicknessE;
0104 
0105   // zero average geometry
0106   memset(&averageGeometry_, 0, sizeof(pixelCPEforGPU::AverageGeometry));
0107 
0108   uint32_t oldLayer = 0;
0109   uint32_t oldLadder = 0;
0110   float rl = 0;
0111   float zl = 0;
0112   float miz = 500, mxz = 0;
0113   float pl = 0;
0114   int nl = 0;
0115   detParamsGPU_.resize(m_DetParams.size());
0116 
0117   for (auto i = 0U; i < m_DetParams.size(); ++i) {
0118     auto& p = m_DetParams[i];
0119     auto& g = detParamsGPU_[i];
0120 
0121     if (!isPhase2_) {
0122       g.nRowsRoc = phase1PixelTopology::numRowsInRoc;
0123       g.nColsRoc = phase1PixelTopology::numColsInRoc;
0124       g.nRows = phase1PixelTopology::numRowsInModule;
0125       g.nCols = phase1PixelTopology::numColsInModule;
0126 
0127       g.numPixsInModule = g.nRows * g.nCols;
0128 
0129     } else {
0130       g.nRowsRoc = p.theDet->specificTopology().rowsperroc();
0131       g.nColsRoc = p.theDet->specificTopology().colsperroc();
0132       g.nRows = p.theDet->specificTopology().rocsX() * g.nRowsRoc;
0133       g.nCols = p.theDet->specificTopology().rocsY() * g.nColsRoc;
0134 
0135       g.numPixsInModule = g.nRows * g.nCols;
0136     }
0137 
0138     assert(p.theDet->index() == int(i));
0139     assert(commonParamsGPU_.thePitchY == p.thePitchY);
0140     assert(commonParamsGPU_.thePitchX == p.thePitchX);
0141 
0142     g.isBarrel = GeomDetEnumerators::isBarrel(p.thePart);
0143     g.isPosZ = p.theDet->surface().position().z() > 0;
0144     g.layer = ttopo_.layer(p.theDet->geographicalId());
0145     g.index = i;  // better be!
0146     g.rawId = p.theDet->geographicalId();
0147     auto thickness = g.isBarrel ? commonParamsGPU_.theThicknessB : commonParamsGPU_.theThicknessE;
0148     assert(thickness == p.theThickness);
0149 
0150     auto ladder = ttopo_.pxbLadder(p.theDet->geographicalId());
0151     if (oldLayer != g.layer) {
0152       oldLayer = g.layer;
0153       LogDebug("PixelCPEFast") << "new layer at " << i << (g.isBarrel ? " B  " : (g.isPosZ ? " E+ " : " E- "))
0154                                << g.layer << " starting at " << g.rawId << '\n'
0155                                << "old layer had " << nl << " ladders";
0156       nl = 0;
0157     }
0158     if (oldLadder != ladder) {
0159       oldLadder = ladder;
0160       LogDebug("PixelCPEFast") << "new ladder at " << i << (g.isBarrel ? " B  " : (g.isPosZ ? " E+ " : " E- "))
0161                                << ladder << " starting at " << g.rawId << '\n'
0162                                << "old ladder ave z,r,p mz " << zl / 8.f << " " << rl / 8.f << " " << pl / 8.f << ' '
0163                                << miz << ' ' << mxz;
0164       rl = 0;
0165       zl = 0;
0166       pl = 0;
0167       miz = isPhase2_ ? 500 : 90;
0168       mxz = 0;
0169       nl++;
0170     }
0171 
0172     g.shiftX = 0.5f * p.lorentzShiftInCmX;
0173     g.shiftY = 0.5f * p.lorentzShiftInCmY;
0174     g.chargeWidthX = p.lorentzShiftInCmX * p.widthLAFractionX;
0175     g.chargeWidthY = p.lorentzShiftInCmY * p.widthLAFractionY;
0176 
0177     g.x0 = p.theOrigin.x();
0178     g.y0 = p.theOrigin.y();
0179     g.z0 = p.theOrigin.z();
0180 
0181     auto vv = p.theDet->surface().position();
0182     auto rr = pixelCPEforGPU::Rotation(p.theDet->surface().rotation());
0183     g.frame = pixelCPEforGPU::Frame(vv.x(), vv.y(), vv.z(), rr);
0184 
0185     zl += vv.z();
0186     miz = std::min(miz, std::abs(vv.z()));
0187     mxz = std::max(mxz, std::abs(vv.z()));
0188     rl += vv.perp();
0189     pl += vv.phi();  // (not obvious)
0190 
0191     // errors .....
0192     ClusterParamGeneric cp;
0193 
0194     cp.with_track_angle = false;
0195 
0196     auto lape = p.theDet->localAlignmentError();
0197     if (lape.invalid())
0198       lape = LocalError();  // zero....
0199 
0200     g.apeXX = lape.xx();
0201     g.apeYY = lape.yy();
0202 
0203     auto toMicron = [&](float x) { return std::min(511, int(x * 1.e4f + 0.5f)); };
0204 
0205     // average angle
0206     auto gvx = p.theOrigin.x() + 40.f * commonParamsGPU_.thePitchX;
0207     auto gvy = p.theOrigin.y();
0208     auto gvz = 1.f / p.theOrigin.z();
0209     //--- Note that the normalization is not required as only the ratio used
0210 
0211     {
0212       // calculate angles (fed into errorFromTemplates)
0213       cp.cotalpha = gvx * gvz;
0214       cp.cotbeta = gvy * gvz;
0215 
0216       if (!isPhase2_)
0217         errorFromTemplates(p, cp, 20000.);
0218       else
0219         cp.qBin_ = 0.f;
0220     }
0221 
0222 #ifdef EDM_ML_DEBUG
0223     auto m = 10000.f;
0224     for (float qclus = 15000; qclus < 35000; qclus += 15000) {
0225       errorFromTemplates(p, cp, qclus);
0226       LogDebug("PixelCPEFast") << i << ' ' << qclus << ' ' << cp.pixmx << ' ' << m * cp.sigmax << ' ' << m * cp.sx1
0227                                << ' ' << m * cp.sx2 << ' ' << m * cp.sigmay << ' ' << m * cp.sy1 << ' ' << m * cp.sy2;
0228     }
0229     LogDebug("PixelCPEFast") << i << ' ' << m * std::sqrt(lape.xx()) << ' ' << m * std::sqrt(lape.yy());
0230 #endif  // EDM_ML_DEBUG
0231 
0232     g.pixmx = std::max(0, cp.pixmx);
0233     g.sx2 = toMicron(cp.sx2);
0234     g.sy1 = std::max(21, toMicron(cp.sy1));  // for some angles sy1 is very small
0235     g.sy2 = std::max(55, toMicron(cp.sy2));  // sometimes sy2 is smaller than others (due to angle?)
0236 
0237     // sample xerr as function of position
0238     auto const xoff = float(phase1PixelTopology::xOffset) * commonParamsGPU_.thePitchX;
0239 
0240     for (int ix = 0; ix < CPEFastParametrisation::kNumErrorBins; ++ix) {
0241       auto x = xoff * (1.f - (0.5f + float(ix)) / 8.f);
0242       auto gvx = p.theOrigin.x() - x;
0243       auto gvy = p.theOrigin.y();
0244       auto gvz = 1.f / p.theOrigin.z();
0245       cp.cotbeta = gvy * gvz;
0246       cp.cotalpha = gvx * gvz;
0247       errorFromTemplates(p, cp, 20000.f);
0248       g.sigmax[ix] = toMicron(cp.sigmax);
0249       g.sigmax1[ix] = toMicron(cp.sx1);
0250       LogDebug("PixelCPEFast") << "sigmax vs x " << i << ' ' << x << ' ' << cp.cotalpha << ' ' << int(g.sigmax[ix])
0251                                << ' ' << int(g.sigmax1[ix]) << ' ' << 10000.f * cp.sigmay << std::endl;
0252     }
0253 #ifdef EDM_ML_DEBUG
0254     // sample yerr as function of position
0255     auto const yoff = float(phase1PixelTopology::yOffset) * commonParamsGPU_.thePitchY;
0256     for (int ix = 0; ix < CPEFastParametrisation::kNumErrorBins; ++ix) {
0257       auto y = yoff * (1.f - (0.5f + float(ix)) / 8.f);
0258       auto gvx = p.theOrigin.x() + 40.f * commonParamsGPU_.thePitchY;
0259       auto gvy = p.theOrigin.y() - y;
0260       auto gvz = 1.f / p.theOrigin.z();
0261       cp.cotbeta = gvy * gvz;
0262       cp.cotalpha = gvx * gvz;
0263       errorFromTemplates(p, cp, 20000.f);
0264       LogDebug("PixelCPEFast") << "sigmay vs y " << i << ' ' << y << ' ' << cp.cotbeta << ' ' << 10000.f * cp.sigmay
0265                                << std::endl;
0266     }
0267 #endif  // EDM_ML_DEBUG
0268 
0269     // calculate angles (repeated)
0270     cp.cotalpha = gvx * gvz;
0271     cp.cotbeta = gvy * gvz;
0272     auto aveCB = cp.cotbeta;
0273 
0274     // sample x by charge
0275     int qbin = CPEFastParametrisation::kGenErrorQBins;  // low charge
0276     int k = 0;
0277     for (int qclus = 1000; qclus < 200000; qclus += 1000) {
0278       errorFromTemplates(p, cp, qclus);
0279       if (cp.qBin_ == qbin)
0280         continue;
0281       qbin = cp.qBin_;
0282       g.xfact[k] = cp.sigmax;
0283       g.yfact[k] = cp.sigmay;
0284       g.minCh[k++] = qclus;
0285 #ifdef EDM_ML_DEBUG
0286       LogDebug("PixelCPEFast") << i << ' ' << g.rawId << ' ' << cp.cotalpha << ' ' << qclus << ' ' << cp.qBin_ << ' '
0287                                << cp.pixmx << ' ' << m * cp.sigmax << ' ' << m * cp.sx1 << ' ' << m * cp.sx2 << ' '
0288                                << m * cp.sigmay << ' ' << m * cp.sy1 << ' ' << m * cp.sy2 << std::endl;
0289 #endif  // EDM_ML_DEBUG
0290     }
0291 
0292     assert(k <= CPEFastParametrisation::kGenErrorQBins);
0293 
0294     // fill the rest  (sometimes bin 4 is missing)
0295     for (int kk = k; kk < CPEFastParametrisation::kGenErrorQBins; ++kk) {
0296       g.xfact[kk] = g.xfact[k - 1];
0297       g.yfact[kk] = g.yfact[k - 1];
0298       g.minCh[kk] = g.minCh[k - 1];
0299     }
0300     auto detx = 1.f / g.xfact[0];
0301     auto dety = 1.f / g.yfact[0];
0302     for (int kk = 0; kk < CPEFastParametrisation::kGenErrorQBins; ++kk) {
0303       g.xfact[kk] *= detx;
0304       g.yfact[kk] *= dety;
0305     }
0306     // sample y in "angle"  (estimated from cluster size)
0307     float ys = 8.f - 4.f;  // apperent bias of half pixel (see plot)
0308     // plot: https://indico.cern.ch/event/934821/contributions/3974619/attachments/2091853/3515041/DigilessReco.pdf page 25
0309     // sample yerr as function of "size"
0310     for (int iy = 0; iy < CPEFastParametrisation::kNumErrorBins; ++iy) {
0311       ys += 1.f;  // first bin 0 is for size 9  (and size is in fixed point 2^3)
0312       if (CPEFastParametrisation::kNumErrorBins - 1 == iy)
0313         ys += 8.f;  // last bin for "overflow"
0314       // cp.cotalpha = ys*(commonParamsGPU_.thePitchX/(8.f*thickness));  //  use this to print sampling in "x"  (and comment the line below)
0315       cp.cotbeta = std::copysign(ys * (commonParamsGPU_.thePitchY / (8.f * thickness)), aveCB);
0316       errorFromTemplates(p, cp, 20000.f);
0317       g.sigmay[iy] = toMicron(cp.sigmay);
0318       LogDebug("PixelCPEFast") << "sigmax/sigmay " << i << ' ' << (ys + 4.f) / 8.f << ' ' << cp.cotalpha << '/'
0319                                << cp.cotbeta << ' ' << 10000.f * cp.sigmax << '/' << int(g.sigmay[iy]) << std::endl;
0320     }
0321   }  // loop over det
0322 
0323   const int numberOfModulesInLadder =
0324       isPhase2_ ? int(phase2PixelTopology::numberOfModulesInLadder) : int(phase1PixelTopology::numberOfModulesInLadder);
0325   const int numberOfModulesInBarrel =
0326       isPhase2_ ? int(phase2PixelTopology::numberOfModulesInBarrel) : int(phase1PixelTopology::numberOfModulesInBarrel);
0327   const int numberOfLaddersInBarrel = commonParamsGPU_.numberOfLaddersInBarrel;
0328 
0329   const int firstEndcapPos = 4, firstEndcapNeg = isPhase2_ ? 16 : 7;
0330   const float ladderFactor = 1.f / float(numberOfModulesInLadder);
0331 
0332   // compute ladder baricenter (only in global z) for the barrel
0333   //
0334   auto& aveGeom = averageGeometry_;
0335   int il = 0;
0336   for (int im = 0, nm = numberOfModulesInBarrel; im < nm; ++im) {
0337     auto const& g = detParamsGPU_[im];
0338     il = im / numberOfModulesInLadder;
0339     assert(il < int(numberOfLaddersInBarrel));
0340     auto z = g.frame.z();
0341     aveGeom.ladderZ[il] += ladderFactor * z;
0342     aveGeom.ladderMinZ[il] = std::min(aveGeom.ladderMinZ[il], z);
0343     aveGeom.ladderMaxZ[il] = std::max(aveGeom.ladderMaxZ[il], z);
0344     aveGeom.ladderX[il] += ladderFactor * g.frame.x();
0345     aveGeom.ladderY[il] += ladderFactor * g.frame.y();
0346     aveGeom.ladderR[il] += ladderFactor * sqrt(g.frame.x() * g.frame.x() + g.frame.y() * g.frame.y());
0347   }
0348   assert(il + 1 == int(numberOfLaddersInBarrel));
0349   // add half_module and tollerance
0350   const float module_length = isPhase2_ ? 4.345f : 6.7f;
0351   constexpr float module_tolerance = 0.2f;
0352   for (int il = 0, nl = numberOfLaddersInBarrel; il < nl; ++il) {
0353     aveGeom.ladderMinZ[il] -= (0.5f * module_length - module_tolerance);
0354     aveGeom.ladderMaxZ[il] += (0.5f * module_length - module_tolerance);
0355   }
0356 
0357   // compute "max z" for first layer in endcap (should we restrict to the outermost ring?)
0358   if (!isPhase2_) {
0359     for (auto im = phase1PixelTopology::layerStart[firstEndcapPos];
0360          im < phase1PixelTopology::layerStart[firstEndcapPos + 1];
0361          ++im) {
0362       auto const& g = detParamsGPU_[im];
0363       aveGeom.endCapZ[0] = std::max(aveGeom.endCapZ[0], g.frame.z());
0364     }
0365     for (auto im = phase1PixelTopology::layerStart[firstEndcapNeg];
0366          im < phase1PixelTopology::layerStart[firstEndcapNeg + 1];
0367          ++im) {
0368       auto const& g = detParamsGPU_[im];
0369       aveGeom.endCapZ[1] = std::min(aveGeom.endCapZ[1], g.frame.z());
0370     }
0371     // correct for outer ring being closer
0372     aveGeom.endCapZ[0] -= 1.5f;
0373     aveGeom.endCapZ[1] += 1.5f;
0374   } else {
0375     for (auto im = phase2PixelTopology::layerStart[firstEndcapPos];
0376          im < phase2PixelTopology::layerStart[firstEndcapPos + 1];
0377          ++im) {
0378       auto const& g = detParamsGPU_[im];
0379       aveGeom.endCapZ[0] = std::max(aveGeom.endCapZ[0], g.frame.z());
0380     }
0381     for (auto im = phase2PixelTopology::layerStart[firstEndcapNeg];
0382          im < phase2PixelTopology::layerStart[firstEndcapNeg + 1];
0383          ++im) {
0384       auto const& g = detParamsGPU_[im];
0385       aveGeom.endCapZ[1] = std::min(aveGeom.endCapZ[1], g.frame.z());
0386     }
0387   }
0388 #ifdef EDM_ML_DEBUG
0389   for (int jl = 0, nl = numberOfLaddersInBarrel; jl < nl; ++jl) {
0390     LogDebug("PixelCPEFast") << jl << ':' << aveGeom.ladderR[jl] << '/'
0391                              << std::sqrt(aveGeom.ladderX[jl] * aveGeom.ladderX[jl] +
0392                                           aveGeom.ladderY[jl] * aveGeom.ladderY[jl])
0393                              << ',' << aveGeom.ladderZ[jl] << ',' << aveGeom.ladderMinZ[jl] << ','
0394                              << aveGeom.ladderMaxZ[jl] << '\n';
0395   }
0396   LogDebug("PixelCPEFast") << aveGeom.endCapZ[0] << ' ' << aveGeom.endCapZ[1];
0397 #endif  // EDM_ML_DEBUG
0398 
0399   // fill Layer and ladders geometry
0400   memset(&layerGeometry_, 0, sizeof(pixelCPEforGPU::LayerGeometry));
0401   if (!isPhase2_) {
0402     memcpy(layerGeometry_.layerStart, phase1PixelTopology::layerStart, sizeof(phase1PixelTopology::layerStart));
0403     memcpy(layerGeometry_.layer, phase1PixelTopology::layer.data(), phase1PixelTopology::layer.size());
0404     layerGeometry_.maxModuleStride = phase1PixelTopology::maxModuleStride;
0405   } else {
0406     memcpy(layerGeometry_.layerStart, phase2PixelTopology::layerStart, sizeof(phase2PixelTopology::layerStart));
0407     memcpy(layerGeometry_.layer, phase2PixelTopology::layer.data(), phase2PixelTopology::layer.size());
0408     layerGeometry_.maxModuleStride = phase2PixelTopology::maxModuleStride;
0409   }
0410 }
0411 
0412 PixelCPEFast::GPUData::~GPUData() {
0413   if (paramsOnGPU_d != nullptr) {
0414     cudaFree((void*)paramsOnGPU_h.m_commonParams);
0415     cudaFree((void*)paramsOnGPU_h.m_detParams);
0416     cudaFree((void*)paramsOnGPU_h.m_averageGeometry);
0417     cudaFree((void*)paramsOnGPU_h.m_layerGeometry);
0418     cudaFree(paramsOnGPU_d);
0419   }
0420 }
0421 
0422 void PixelCPEFast::errorFromTemplates(DetParam const& theDetParam,
0423                                       ClusterParamGeneric& theClusterParam,
0424                                       float qclus) const {
0425   float locBz = theDetParam.bz;
0426   float locBx = theDetParam.bx;
0427   LogDebug("PixelCPEFast") << "PixelCPEFast::localPosition(...) : locBz = " << locBz;
0428 
0429   theClusterParam.pixmx = std::numeric_limits<int>::max();  // max pixel charge for truncation of 2-D cluster
0430 
0431   theClusterParam.sigmay = -999.9;  // CPE Generic y-error for multi-pixel cluster
0432   theClusterParam.sigmax = -999.9;  // CPE Generic x-error for multi-pixel cluster
0433   theClusterParam.sy1 = -999.9;     // CPE Generic y-error for single single-pixel
0434   theClusterParam.sy2 = -999.9;     // CPE Generic y-error for single double-pixel cluster
0435   theClusterParam.sx1 = -999.9;     // CPE Generic x-error for single single-pixel cluster
0436   theClusterParam.sx2 = -999.9;     // CPE Generic x-error for single double-pixel cluster
0437 
0438   float dummy;
0439 
0440   SiPixelGenError gtempl(thePixelGenError_);
0441   int gtemplID = theDetParam.detTemplateId;
0442 
0443   theClusterParam.qBin_ = gtempl.qbin(gtemplID,
0444                                       theClusterParam.cotalpha,
0445                                       theClusterParam.cotbeta,
0446                                       locBz,
0447                                       locBx,
0448                                       qclus,
0449                                       false,
0450                                       theClusterParam.pixmx,
0451                                       theClusterParam.sigmay,
0452                                       dummy,
0453                                       theClusterParam.sigmax,
0454                                       dummy,
0455                                       theClusterParam.sy1,
0456                                       dummy,
0457                                       theClusterParam.sy2,
0458                                       dummy,
0459                                       theClusterParam.sx1,
0460                                       dummy,
0461                                       theClusterParam.sx2,
0462                                       dummy);
0463 
0464   theClusterParam.sigmax = theClusterParam.sigmax * micronsToCm;
0465   theClusterParam.sx1 = theClusterParam.sx1 * micronsToCm;
0466   theClusterParam.sx2 = theClusterParam.sx2 * micronsToCm;
0467 
0468   theClusterParam.sigmay = theClusterParam.sigmay * micronsToCm;
0469   theClusterParam.sy1 = theClusterParam.sy1 * micronsToCm;
0470   theClusterParam.sy2 = theClusterParam.sy2 * micronsToCm;
0471 }
0472 
0473 //-----------------------------------------------------------------------------
0474 //! Hit position in the local frame (in cm).  Unlike other CPE's, this
0475 //! one converts everything from the measurement frame (in channel numbers)
0476 //! into the local frame (in centimeters).
0477 //-----------------------------------------------------------------------------
0478 LocalPoint PixelCPEFast::localPosition(DetParam const& theDetParam, ClusterParam& theClusterParamBase) const {
0479   ClusterParamGeneric& theClusterParam = static_cast<ClusterParamGeneric&>(theClusterParamBase);
0480 
0481   assert(!theClusterParam.with_track_angle);
0482 
0483   if (useErrorsFromTemplates_) {
0484     errorFromTemplates(theDetParam, theClusterParam, theClusterParam.theCluster->charge());
0485   } else {
0486     theClusterParam.qBin_ = 0;
0487   }
0488 
0489   int q_f_X;  //!< Q of the first  pixel  in X
0490   int q_l_X;  //!< Q of the last   pixel  in X
0491   int q_f_Y;  //!< Q of the first  pixel  in Y
0492   int q_l_Y;  //!< Q of the last   pixel  in Y
0493   collect_edge_charges(theClusterParam, q_f_X, q_l_X, q_f_Y, q_l_Y, useErrorsFromTemplates_ && truncatePixelCharge_);
0494 
0495   // do GPU like ...
0496   pixelCPEforGPU::ClusParams cp;
0497 
0498   cp.minRow[0] = theClusterParam.theCluster->minPixelRow();
0499   cp.maxRow[0] = theClusterParam.theCluster->maxPixelRow();
0500   cp.minCol[0] = theClusterParam.theCluster->minPixelCol();
0501   cp.maxCol[0] = theClusterParam.theCluster->maxPixelCol();
0502 
0503   cp.q_f_X[0] = q_f_X;
0504   cp.q_l_X[0] = q_l_X;
0505   cp.q_f_Y[0] = q_f_Y;
0506   cp.q_l_Y[0] = q_l_Y;
0507 
0508   cp.charge[0] = theClusterParam.theCluster->charge();
0509 
0510   auto ind = theDetParam.theDet->index();
0511   pixelCPEforGPU::position(commonParamsGPU_, detParamsGPU_[ind], cp, 0);
0512   auto xPos = cp.xpos[0];
0513   auto yPos = cp.ypos[0];
0514 
0515   // set the error  (mind ape....)
0516   pixelCPEforGPU::errorFromDB(commonParamsGPU_, detParamsGPU_[ind], cp, 0);
0517   theClusterParam.sigmax = cp.xerr[0];
0518   theClusterParam.sigmay = cp.yerr[0];
0519 
0520   LogDebug("PixelCPEFast") << " in PixelCPEFast:localPosition - pos = " << xPos << " " << yPos << " size "
0521                            << cp.maxRow[0] - cp.minRow[0] << ' ' << cp.maxCol[0] - cp.minCol[0];
0522 
0523   //--- Now put the two together
0524   LocalPoint pos_in_local(xPos, yPos);
0525   return pos_in_local;
0526 }
0527 
0528 //==============  INFLATED ERROR AND ERRORS FROM DB BELOW  ================
0529 
0530 //-------------------------------------------------------------------------
0531 //  Hit error in the local frame
0532 //-------------------------------------------------------------------------
0533 LocalError PixelCPEFast::localError(DetParam const& theDetParam, ClusterParam& theClusterParamBase) const {
0534   ClusterParamGeneric& theClusterParam = static_cast<ClusterParamGeneric&>(theClusterParamBase);
0535 
0536   auto xerr = theClusterParam.sigmax;
0537   auto yerr = theClusterParam.sigmay;
0538 
0539   LogDebug("PixelCPEFast") << " errors  " << xerr << " " << yerr;
0540 
0541   auto xerr_sq = xerr * xerr;
0542   auto yerr_sq = yerr * yerr;
0543 
0544   return LocalError(xerr_sq, 0, yerr_sq);
0545 }
0546 
0547 void PixelCPEFast::fillPSetDescription(edm::ParameterSetDescription& desc) {
0548   // call PixelCPEGenericBase fillPSetDescription to add common rechit errors
0549   PixelCPEGenericBase::fillPSetDescription(desc);
0550   desc.add<bool>("isPhase2", false);
0551 }