Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:26:17 

0001 #include "RecoLocalTracker/SiPhase2VectorHitBuilder/interface/VectorHitBuilderAlgorithm.h"
0002 #include "RecoLocalTracker/ClusterParameterEstimator/interface/ClusterParameterEstimator.h"
0003 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0004 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0005 #include "DataFormats/TrackerRecHit2D/interface/VectorHit2D.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 
0008 void VectorHitBuilderAlgorithm::run(edm::Handle<edmNew::DetSetVector<Phase2TrackerCluster1D>> clusters,
0009                                     VectorHitCollection& vhAcc,
0010                                     VectorHitCollection& vhRej,
0011                                     edmNew::DetSetVector<Phase2TrackerCluster1D>& clustersAcc,
0012                                     edmNew::DetSetVector<Phase2TrackerCluster1D>& clustersRej) const {
0013   LogDebug("VectorHitBuilderAlgorithm") << "Run VectorHitBuilderAlgorithm ... \n";
0014   const auto* clustersPhase2Collection = clusters.product();
0015 
0016   //loop over the DetSetVector
0017   LogDebug("VectorHitBuilderAlgorithm") << "with #clusters : " << clustersPhase2Collection->size() << std::endl;
0018   for (auto dSViter : *clustersPhase2Collection) {
0019     unsigned int rawDetId1(dSViter.detId());
0020     DetId detId1(rawDetId1);
0021     DetId lowerDetId, upperDetId;
0022     if (tkTopo_->isLower(detId1)) {
0023       lowerDetId = detId1;
0024       upperDetId = tkTopo_->partnerDetId(detId1);
0025     } else
0026       continue;
0027 
0028     DetId detIdStack = tkTopo_->stack(detId1);
0029 
0030     //debug
0031     LogDebug("VectorHitBuilderAlgorithm") << "  DetId stack : " << detIdStack.rawId() << std::endl;
0032     LogDebug("VectorHitBuilderAlgorithm") << "  DetId lower set of clusters  : " << lowerDetId.rawId();
0033     LogDebug("VectorHitBuilderAlgorithm") << "  DetId upper set of clusters  : " << upperDetId.rawId() << std::endl;
0034 
0035     const GeomDet* gd;
0036     const StackGeomDet* stackDet;
0037     const auto& it_detLower = dSViter;
0038     const auto& it_detUpper = clustersPhase2Collection->find(upperDetId);
0039 
0040     if (it_detUpper != clustersPhase2Collection->end()) {
0041       gd = tkGeom_->idToDet(detIdStack);
0042       stackDet = dynamic_cast<const StackGeomDet*>(gd);
0043       buildVectorHits(vhAcc, vhRej, detIdStack, stackDet, clusters, it_detLower, *it_detUpper);
0044     }
0045   }
0046   LogDebug("VectorHitBuilderAlgorithm") << "End run VectorHitBuilderAlgorithm ... \n";
0047 }
0048 
0049 bool VectorHitBuilderAlgorithm::checkClustersCompatibilityBeforeBuilding(
0050     edm::Handle<edmNew::DetSetVector<Phase2TrackerCluster1D>> clusters,
0051     const Detset& theLowerDetSet,
0052     const Detset& theUpperDetSet) const {
0053   if (theLowerDetSet.size() == 1 && theUpperDetSet.size() == 1)
0054     return true;
0055 
0056   //order lower clusters in u
0057   std::vector<Phase2TrackerCluster1D> lowerClusters;
0058   lowerClusters.reserve(theLowerDetSet.size());
0059   if (theLowerDetSet.size() > 1)
0060     LogDebug("VectorHitBuilderAlgorithm") << " more than 1 lower cluster! " << std::endl;
0061   if (theUpperDetSet.size() > 1)
0062     LogDebug("VectorHitBuilderAlgorithm") << " more than 1 upper cluster! " << std::endl;
0063   for (const_iterator cil = theLowerDetSet.begin(); cil != theLowerDetSet.end(); ++cil) {
0064     Phase2TrackerCluster1DRef clusterLower = edmNew::makeRefTo(clusters, cil);
0065     lowerClusters.push_back(*clusterLower);
0066   }
0067   return true;
0068 }
0069 
0070 bool VectorHitBuilderAlgorithm::checkClustersCompatibility(Local3DPoint& poslower,
0071                                                            Local3DPoint& posupper,
0072                                                            LocalError& errlower,
0073                                                            LocalError& errupper) const {
0074   return true;
0075 }
0076 
0077 //----------------------------------------------------------------------------
0078 //ERICA::in the DT code the global position is used to compute the alpha angle and put a cut on that.
0079 void VectorHitBuilderAlgorithm::buildVectorHits(VectorHitCollection& vhAcc,
0080                                                 VectorHitCollection& vhRej,
0081                                                 DetId detIdStack,
0082                                                 const StackGeomDet* stack,
0083                                                 edm::Handle<edmNew::DetSetVector<Phase2TrackerCluster1D>> clusters,
0084                                                 const Detset& theLowerDetSet,
0085                                                 const Detset& theUpperDetSet,
0086                                                 const std::vector<bool>& phase2OTClustersToSkip) const {
0087   if (checkClustersCompatibilityBeforeBuilding(clusters, theLowerDetSet, theUpperDetSet)) {
0088     LogDebug("VectorHitBuilderAlgorithm") << "  compatible -> continue ... " << std::endl;
0089   } else {
0090     LogTrace("VectorHitBuilderAlgorithm") << "  not compatible, going to the next cluster";
0091   }
0092 
0093   edmNew::DetSetVector<VectorHit>::FastFiller vh_colAcc(vhAcc, detIdStack);
0094   edmNew::DetSetVector<VectorHit>::FastFiller vh_colRej(vhRej, detIdStack);
0095 
0096   unsigned int layerStack = tkTopo_->layer(stack->geographicalId());
0097   if (stack->subDetector() == GeomDetEnumerators::SubDetector::P2OTB)
0098     LogDebug("VectorHitBuilderAlgorithm") << " \t is barrel.    " << std::endl;
0099   if (stack->subDetector() == GeomDetEnumerators::SubDetector::P2OTEC)
0100     LogDebug("VectorHitBuilderAlgorithm") << " \t is endcap.    " << std::endl;
0101   LogDebug("VectorHitBuilderAlgorithm") << " \t layer is : " << layerStack << std::endl;
0102 
0103   float cut = 0.0;
0104   if (stack->subDetector() == GeomDetEnumerators::SubDetector::P2OTB)
0105     cut = barrelCut_.at(layerStack);
0106   if (stack->subDetector() == GeomDetEnumerators::SubDetector::P2OTEC)
0107     cut = endcapCut_.at(layerStack);
0108   LogDebug("VectorHitBuilderAlgorithm") << " \t the cut is:" << cut << std::endl;
0109 
0110   //only cache local parameters for upper cluster as we loop over lower clusters only once anyway
0111   std::vector<std::pair<LocalPoint, LocalError>> localParamsUpper;
0112   std::vector<const PixelGeomDetUnit*> localGDUUpper;
0113 
0114   const PixelGeomDetUnit* gduUpp = dynamic_cast<const PixelGeomDetUnit*>(stack->upperDet());
0115   for (auto const& clusterUpper : theUpperDetSet) {
0116     localGDUUpper.push_back(gduUpp);
0117     localParamsUpper.push_back(cpe_->localParameters(clusterUpper, *gduUpp));
0118   }
0119   int upperIterator = 0;
0120   const PixelGeomDetUnit* gduLow = dynamic_cast<const PixelGeomDetUnit*>(stack->lowerDet());
0121   for (const_iterator cil = theLowerDetSet.begin(); cil != theLowerDetSet.end(); ++cil) {
0122     LogDebug("VectorHitBuilderAlgorithm") << " lower clusters " << std::endl;
0123     Phase2TrackerCluster1DRef cluL = edmNew::makeRefTo(clusters, cil);
0124 #ifdef EDM_ML_DEBUG
0125     printCluster(stack->lowerDet(), &*cluL);
0126 #endif
0127     auto&& lparamsLow = cpe_->localParameters(*cluL, *gduLow);
0128     upperIterator = 0;
0129     for (const_iterator ciu = theUpperDetSet.begin(); ciu != theUpperDetSet.end(); ++ciu) {
0130       LogDebug("VectorHitBuilderAlgorithm") << "\t upper clusters " << std::endl;
0131       Phase2TrackerCluster1DRef cluU = edmNew::makeRefTo(clusters, ciu);
0132 #ifdef EDM_ML_DEBUG
0133       printCluster(stack->upperDet(), &*cluU);
0134 #endif
0135       //applying the parallax correction
0136       double pC = computeParallaxCorrection(
0137           gduLow, lparamsLow.first, localGDUUpper[upperIterator], localParamsUpper[upperIterator].first);
0138       LogDebug("VectorHitBuilderAlgorithm") << " \t parallax correction:" << pC << std::endl;
0139       double lpos_upp_corr = 0.0;
0140       double lpos_low_corr = 0.0;
0141       auto const localUpperX = localParamsUpper[upperIterator].first.x();
0142       if (localUpperX > lparamsLow.first.x()) {
0143         if (localUpperX > 0) {
0144           lpos_low_corr = lparamsLow.first.x();
0145           lpos_upp_corr = localParamsUpper[upperIterator].first.x() - std::abs(pC);
0146         } else if (localUpperX < 0) {
0147           lpos_low_corr = lparamsLow.first.x() + std::abs(pC);
0148           lpos_upp_corr = localUpperX;
0149         }
0150       } else if (localUpperX < lparamsLow.first.x()) {
0151         if (localUpperX > 0) {
0152           lpos_low_corr = lparamsLow.first.x() - std::abs(pC);
0153           lpos_upp_corr = localUpperX;
0154         } else if (localUpperX < 0) {
0155           lpos_low_corr = lparamsLow.first.x();
0156           lpos_upp_corr = localUpperX + std::abs(pC);
0157         }
0158       } else {
0159         if (localUpperX > 0) {
0160           lpos_low_corr = lparamsLow.first.x();
0161           lpos_upp_corr = localUpperX - std::abs(pC);
0162         } else if (localUpperX < 0) {
0163           lpos_low_corr = lparamsLow.first.x();
0164           lpos_upp_corr = localUpperX + std::abs(pC);
0165         }
0166       }
0167 
0168       LogDebug("VectorHitBuilderAlgorithm") << " \t local pos upper corrected (x):" << lpos_upp_corr << std::endl;
0169       LogDebug("VectorHitBuilderAlgorithm") << " \t local pos lower corrected (x):" << lpos_low_corr << std::endl;
0170 
0171       double width = lpos_low_corr - lpos_upp_corr;
0172       LogDebug("VectorHitBuilderAlgorithm") << " \t width: " << width << std::endl;
0173 
0174       //old cut: indipendent from layer
0175       //building my tolerance : 10*sigma
0176       //double delta = 10.0 * sqrt(lparamsLow.second.xx() + localParamsUpper[upperIterator].second.xx());
0177       //LogDebug("VectorHitBuilderAlgorithm") << " \t delta: " << delta << std::endl;
0178       //if( (lpos_upp_corr < lpos_low_corr + delta) &&
0179       //    (lpos_upp_corr > lpos_low_corr - delta) ){
0180       //new cut: dependent on layers
0181       if (std::abs(width) < cut) {
0182         LogDebug("VectorHitBuilderAlgorithm") << " accepting VH! " << std::endl;
0183         VectorHit vh = buildVectorHit(stack, cluL, cluU);
0184         //protection: the VH can also be empty!!
0185         if (vh.isValid()) {
0186           vh_colAcc.push_back(vh);
0187         }
0188 
0189       } else {
0190         LogDebug("VectorHitBuilderAlgorithm") << " rejecting VH: " << std::endl;
0191         //storing vh rejected for combinatiorial studies
0192         VectorHit vh = buildVectorHit(stack, cluL, cluU);
0193         if (vh.isValid()) {
0194           vh_colRej.push_back(vh);
0195         }
0196       }
0197       upperIterator += 1;
0198     }
0199   }
0200 }
0201 
0202 VectorHit VectorHitBuilderAlgorithm::buildVectorHit(const StackGeomDet* stack,
0203                                                     Phase2TrackerCluster1DRef lower,
0204                                                     Phase2TrackerCluster1DRef upper) const {
0205   LogTrace("VectorHitBuilderAlgorithm") << "Build VH with: ";
0206 #ifdef EDM_ML_DEBUG
0207   printCluster(stack->upperDet(), &*upper);
0208 #endif
0209   const PixelGeomDetUnit* geomDetLower = static_cast<const PixelGeomDetUnit*>(stack->lowerDet());
0210   const PixelGeomDetUnit* geomDetUpper = static_cast<const PixelGeomDetUnit*>(stack->upperDet());
0211 
0212   auto&& lparamsLower = cpe_->localParameters(*lower, *geomDetLower);  // x, y, z, e2_xx, e2_xy, e2_yy
0213   Global3DPoint gparamsLower = geomDetLower->surface().toGlobal(lparamsLower.first);
0214   LogTrace("VectorHitBuilderAlgorithm") << "\t lower global pos: " << gparamsLower;
0215 
0216   auto&& lparamsUpper = cpe_->localParameters(*upper, *geomDetUpper);
0217   Global3DPoint gparamsUpper = geomDetUpper->surface().toGlobal(lparamsUpper.first);
0218   LogTrace("VectorHitBuilderAlgorithm") << "\t upper global pos: " << gparamsUpper;
0219 
0220   //local parameters of upper cluster in lower system of reference
0221   Local3DPoint lparamsUpperInLower = geomDetLower->surface().toLocal(gparamsUpper);
0222 
0223   LogTrace("VectorHitBuilderAlgorithm") << "\t lower global pos: " << gparamsLower;
0224   LogTrace("VectorHitBuilderAlgorithm") << "\t upper global pos: " << gparamsUpper;
0225 
0226   LogTrace("VectorHitBuilderAlgorithm") << "A:\t lower local pos: " << lparamsLower.first
0227                                         << " with error: " << lparamsLower.second << std::endl;
0228   LogTrace("VectorHitBuilderAlgorithm") << "A:\t upper local pos in the lower sof " << lparamsUpperInLower
0229                                         << " with error: " << lparamsUpper.second << std::endl;
0230 
0231   bool ok =
0232       checkClustersCompatibility(lparamsLower.first, lparamsUpper.first, lparamsLower.second, lparamsUpper.second);
0233 
0234   if (ok) {
0235     AlgebraicSymMatrix22 covMat2Dzx;
0236     double chi22Dzx = 0.0;
0237     Local3DPoint pos2Dzx;
0238     Local3DVector dir2Dzx;
0239     fit2Dzx(lparamsLower.first,
0240             lparamsUpperInLower,
0241             lparamsLower.second,
0242             lparamsUpper.second,
0243             pos2Dzx,
0244             dir2Dzx,
0245             covMat2Dzx,
0246             chi22Dzx);
0247     LogTrace("VectorHitBuilderAlgorithm") << "\t  pos2Dzx: " << pos2Dzx;
0248     LogTrace("VectorHitBuilderAlgorithm") << "\t  dir2Dzx: " << dir2Dzx;
0249     LogTrace("VectorHitBuilderAlgorithm") << "\t  cov2Dzx: " << covMat2Dzx;
0250     VectorHit2D vh2Dzx = VectorHit2D(pos2Dzx, dir2Dzx, covMat2Dzx, chi22Dzx);
0251 
0252     AlgebraicSymMatrix22 covMat2Dzy;
0253     double chi22Dzy = 0.0;
0254     Local3DPoint pos2Dzy;
0255     Local3DVector dir2Dzy;
0256     fit2Dzy(lparamsLower.first,
0257             lparamsUpperInLower,
0258             lparamsLower.second,
0259             lparamsUpper.second,
0260             pos2Dzy,
0261             dir2Dzy,
0262             covMat2Dzy,
0263             chi22Dzy);
0264     LogTrace("VectorHitBuilderAlgorithm") << "\t  pos2Dzy: " << pos2Dzy;
0265     LogTrace("VectorHitBuilderAlgorithm") << "\t  dir2Dzy: " << dir2Dzy;
0266     LogTrace("VectorHitBuilderAlgorithm") << "\t  cov2Dzy: " << covMat2Dzy;
0267     VectorHit2D vh2Dzy = VectorHit2D(pos2Dzy, dir2Dzy, covMat2Dzy, chi22Dzy);
0268 
0269     OmniClusterRef lowerOmni(lower);
0270     OmniClusterRef upperOmni(upper);
0271 
0272     Global3DPoint gPositionLower = VectorHit::phase2clusterGlobalPos(geomDetLower, lower);
0273     Global3DPoint gPositionUpper = VectorHit::phase2clusterGlobalPos(geomDetUpper, upper);
0274     GlobalError gErrorLower = VectorHit::phase2clusterGlobalPosErr(geomDetLower);
0275     GlobalError gErrorUpper = VectorHit::phase2clusterGlobalPosErr(geomDetUpper);
0276 
0277     if (gPositionLower.perp() > gPositionUpper.perp()) {
0278       std::swap(gPositionLower, gPositionUpper);
0279       std::swap(gErrorLower, gErrorUpper);
0280     }
0281 
0282     const CurvatureAndPhi curvatureAndPhi = curvatureANDphi(gPositionLower, gPositionUpper, gErrorLower, gErrorUpper);
0283     VectorHit vh = VectorHit(*stack,
0284                              vh2Dzx,
0285                              vh2Dzy,
0286                              lowerOmni,
0287                              upperOmni,
0288                              curvatureAndPhi.curvature,
0289                              curvatureAndPhi.curvatureError,
0290                              curvatureAndPhi.phi);
0291     return vh;
0292   }
0293 
0294   return VectorHit();
0295 }
0296 
0297 void VectorHitBuilderAlgorithm::fit2Dzx(const Local3DPoint lpCI,
0298                                         const Local3DPoint lpCO,
0299                                         const LocalError leCI,
0300                                         const LocalError leCO,
0301                                         Local3DPoint& pos,
0302                                         Local3DVector& dir,
0303                                         AlgebraicSymMatrix22& covMatrix,
0304                                         double& chi2) const {
0305   float x[2] = {lpCI.z(), lpCO.z()};
0306   float y[2] = {lpCI.x(), lpCO.x()};
0307   float sqCI = sqrt(leCI.xx());
0308   float sqCO = sqrt(leCO.xx());
0309   float sigy2[2] = {sqCI * sqCI, sqCO * sqCO};
0310 
0311   fit(x, y, sigy2, pos, dir, covMatrix, chi2);
0312 
0313   return;
0314 }
0315 
0316 void VectorHitBuilderAlgorithm::fit2Dzy(const Local3DPoint lpCI,
0317                                         const Local3DPoint lpCO,
0318                                         const LocalError leCI,
0319                                         const LocalError leCO,
0320                                         Local3DPoint& pos,
0321                                         Local3DVector& dir,
0322                                         AlgebraicSymMatrix22& covMatrix,
0323                                         double& chi2) const {
0324   float x[2] = {lpCI.z(), lpCO.z()};
0325   float y[2] = {lpCI.y(), lpCO.y()};
0326   float sqCI = sqrt(leCI.yy());
0327   float sqCO = sqrt(leCO.yy());
0328   float sigy2[2] = {sqCI * sqCI, sqCO * sqCO};
0329 
0330   fit(x, y, sigy2, pos, dir, covMatrix, chi2);
0331 
0332   return;
0333 }
0334 
0335 void VectorHitBuilderAlgorithm::fit(float x[2],
0336                                     float y[2],
0337                                     float sigy2[2],
0338                                     Local3DPoint& pos,
0339                                     Local3DVector& dir,
0340                                     AlgebraicSymMatrix22& covMatrix,
0341                                     double& chi2) const {
0342   float slope = 0.;
0343   float intercept = 0.;
0344   float covss = 0.;
0345   float covii = 0.;
0346   float covsi = 0.;
0347 
0348   linearFit(x, y, 2, sigy2, slope, intercept, covss, covii, covsi);
0349 
0350   covMatrix[0][0] = covss;  // this is var(dy/dz)
0351   covMatrix[1][1] = covii;  // this is var(y)
0352   covMatrix[1][0] = covsi;  // this is cov(dy/dz,y)
0353 
0354   for (unsigned int j = 0; j < 2; j++) {
0355     const double ypred = intercept + slope * x[j];
0356     const double dy = (y[j] - ypred) / sqrt(sigy2[j]);
0357     chi2 += dy * dy;
0358   }
0359 
0360   pos = Local3DPoint(intercept, 0., 0.);
0361   //difference in z is the difference of the lowermost and the uppermost cluster z pos
0362   float slopeZ = x[1] - x[0];
0363   dir = LocalVector(slope, 0., slopeZ);
0364 }
0365 
0366 VectorHitBuilderAlgorithm::CurvatureAndPhi VectorHitBuilderAlgorithm::curvatureANDphi(Global3DPoint gPositionLower,
0367                                                                                       Global3DPoint gPositionUpper,
0368                                                                                       GlobalError gErrorLower,
0369                                                                                       GlobalError gErrorUpper) const {
0370   VectorHitBuilderAlgorithm::CurvatureAndPhi result;
0371 
0372   float curvature = -999.;
0373   float errorCurvature = -999.;
0374   float phi = -999.;
0375 
0376   float h1 = gPositionLower.x() * gPositionUpper.y() - gPositionUpper.x() * gPositionLower.y();
0377 
0378   //determine sign of curvature
0379   AlgebraicVector2 n1;
0380   n1[0] = -gPositionLower.y();
0381   n1[1] = gPositionLower.x();
0382   AlgebraicVector2 n2;
0383   n2[0] = gPositionUpper.x() - gPositionLower.x();
0384   n2[1] = gPositionUpper.y() - gPositionLower.y();
0385 
0386   double n3 = n1[0] * n2[0] + n1[1] * n2[1];
0387   double signCurv = -copysign(1.0, n3);
0388   double phi1 = atan2(gPositionUpper.y() - gPositionLower.y(), gPositionUpper.x() - gPositionLower.x());
0389 
0390   double x2Low = pow(gPositionLower.x(), 2);
0391   double y2Low = pow(gPositionLower.y(), 2);
0392   double x2Up = pow(gPositionUpper.x(), 2);
0393   double y2Up = pow(gPositionUpper.y(), 2);
0394 
0395   if (h1 != 0) {
0396     double h2 = 2 * h1;
0397     double h2Inf = 1. / (2 * h1);
0398     double r12 = gPositionLower.perp2();
0399     double r22 = gPositionUpper.perp2();
0400     double h3 = pow(n2[0], 2) + pow(n2[1], 2);
0401     double h4 = -x2Low * gPositionUpper.x() + gPositionLower.x() * x2Up + gPositionLower.x() * y2Up -
0402                 gPositionUpper.x() * y2Low;
0403     double h5 =
0404         x2Low * gPositionUpper.y() - x2Up * gPositionLower.y() + y2Low * gPositionUpper.y() - gPositionLower.y() * y2Up;
0405 
0406     //radius of circle
0407     double invRho2 = (4. * h1 * h1) / (r12 * r22 * h3);
0408     curvature = sqrt(invRho2);
0409 
0410     //center of circle
0411     double xcentre = h5 / h2;
0412     double ycentre = h4 / h2;
0413 
0414     //to compute phi at the cluster points
0415     double xtg = gPositionLower.y() - ycentre;
0416     double ytg = -(gPositionLower.x() - xcentre);
0417 
0418     //to compute phi at the origin
0419     phi = atan2(ytg, xtg);
0420 
0421     AlgebraicROOTObject<4, 4>::Matrix jacobian;
0422 
0423     double denom1 = 1. / sqrt(r12 * r22 * h3);
0424     double denom2 = 1. / (pow(r12 * r22 * h3, 1.5));
0425     jacobian[0][0] = 1.0;  // dx1/dx1 dx1/dy1 dx2/dx1 dy2/dx1
0426     jacobian[1][1] = 1.0;  //dy1/dx1 dy1/dy1 dy2/dx1 dy2/dx1
0427     jacobian[2][0] =
0428         -2. * ((h1 * (gPositionLower.x() * r22 * h3 + (gPositionLower.x() - gPositionUpper.x()) * r12 * r22)) * denom2 -
0429                (gPositionUpper.y()) * denom1);  // dkappa/dx1
0430     jacobian[2][1] =
0431         -2. * ((gPositionUpper.x()) * denom1 +
0432                (h1 * (gPositionLower.y() * r22 * h3 + r12 * r22 * (gPositionLower.y() - gPositionUpper.y()))) *
0433                    denom2);  // dkappa/dy1
0434     jacobian[2][2] =
0435         -2. * ((gPositionLower.y()) * denom1 +
0436                (h1 * (gPositionUpper.x() * r12 * h3 - (gPositionLower.x() - gPositionUpper.x()) * r12 * r22)) *
0437                    denom2);  // dkappa/dx2
0438     jacobian[2][3] =
0439         -2. * ((h1 * (gPositionUpper.y() * r12 * h3 - r12 * r22 * (gPositionLower.y() - gPositionUpper.y()))) * denom2 -
0440                (gPositionLower.x()) * denom1);  // dkappa/dy2
0441     AlgebraicVector2 mVector;
0442     //to compute phi at the cluster points
0443     mVector[0] = (gPositionLower.y() - ycentre) * invRho2;   // dphi/dxcentre
0444     mVector[1] = -(gPositionLower.x() - xcentre) * invRho2;  // dphi/dycentre
0445     //to compute phi at the origin
0446 
0447     double h22Inv = 1. / pow(h2, 2);
0448 
0449     AlgebraicROOTObject<2, 4>::Matrix kMatrix;
0450     kMatrix[0][0] =
0451         2. * ((gPositionLower.x() * gPositionUpper.y()) * h2Inf - (gPositionUpper.y() * h5) * h22Inv);  // dxm/dx1
0452     kMatrix[0][1] = (2. * gPositionUpper.x() * h5) * h22Inv -
0453                     (x2Up + y2Up - 2. * gPositionLower.y() * gPositionUpper.y()) * h2Inf;  // dxm/dy1
0454     kMatrix[0][2] =
0455         2. * ((gPositionLower.y() * h5) * h22Inv - (gPositionUpper.x() * gPositionLower.y()) * h2Inf);  // dxm/dx2
0456     kMatrix[0][3] = (x2Low + y2Low - 2. * gPositionUpper.y() * gPositionLower.y()) * h2Inf -
0457                     (2. * gPositionLower.x() * h5) * h22Inv;  // dxm/dy2
0458     kMatrix[1][0] = (x2Up - 2. * gPositionLower.x() * gPositionUpper.x() + y2Up) * h2Inf -
0459                     (2. * gPositionUpper.y() * h4) * h22Inv;  // dym/dx1
0460     kMatrix[1][1] =
0461         2. * ((gPositionUpper.x() * h4) * h22Inv - (gPositionUpper.x() * gPositionLower.y()) * h2Inf);  // dym/dy1
0462     kMatrix[1][2] = (2. * gPositionLower.y() * h4) * h22Inv -
0463                     (x2Low - 2. * gPositionUpper.x() * gPositionLower.x() + y2Low) * h2Inf;  // dym/dx2
0464     kMatrix[1][3] =
0465         2. * (gPositionLower.x() * gPositionUpper.y()) * h2Inf - (gPositionLower.x() * h4) * h22Inv;  // dym/dy2
0466 
0467     AlgebraicVector4 nMatrix = mVector * kMatrix;
0468     jacobian[3][0] = nMatrix[0];  // dphi/(dx1,dy1,dx2,dy2)
0469     jacobian[3][1] = nMatrix[1];  // dphi/(dx1,dy1,dx2,dy2)
0470     jacobian[3][2] = nMatrix[2];  // dphi/(dx1,dy1,dx2,dy2)
0471     jacobian[3][3] = nMatrix[3];  // dphi/(dx1,dy1,dx2,dy2)
0472 
0473     //assign correct sign to the curvature errors
0474     if ((signCurv < 0 && curvature > 0) || (signCurv > 0 && curvature < 0)) {
0475       curvature = -curvature;
0476       for (int i = 0; i < 4; i++) {
0477         jacobian[2][i] = -jacobian[2][i];
0478       }
0479     }
0480 
0481     // bring phi in the same quadrant as phi1
0482     if (deltaPhi(phi, phi1) > M_PI / 2.) {
0483       phi = phi + M_PI;
0484       if (phi > M_PI)
0485         phi = phi - 2. * M_PI;
0486     }
0487 
0488     //computing the curvature error
0489     AlgebraicVector4 curvatureJacobian;
0490     for (int i = 0; i < 4; i++) {
0491       curvatureJacobian[i] = jacobian[2][i];
0492     }
0493 
0494     AlgebraicROOTObject<4, 4>::Matrix gErrors;
0495 
0496     gErrors[0][0] = gErrorLower.cxx();
0497     gErrors[0][1] = gErrorLower.cyx();
0498     gErrors[1][0] = gErrorLower.cyx();
0499     gErrors[1][1] = gErrorLower.cyy();
0500     gErrors[2][2] = gErrorUpper.cxx();
0501     gErrors[2][3] = gErrorUpper.cyx();
0502     gErrors[3][2] = gErrorUpper.cyx();
0503     gErrors[3][3] = gErrorUpper.cyy();
0504 
0505     AlgebraicVector4 temp = curvatureJacobian;
0506     temp = temp * gErrors;
0507     errorCurvature = temp[0] * curvatureJacobian[0] + temp[1] * curvatureJacobian[1] + temp[2] * curvatureJacobian[2] +
0508                      temp[3] * curvatureJacobian[3];
0509   }
0510 
0511   result.curvature = curvature;
0512   result.curvatureError = errorCurvature;
0513   result.phi = phi;
0514   return result;
0515 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!