diff --git a/RecoPixelVertexing/PixelTriplets/plugins/CACell.h b/RecoPixelVertexing/PixelTriplets/plugins/CACell.h index b561c94f063b7..7a5667eded54c 100644 --- a/RecoPixelVertexing/PixelTriplets/plugins/CACell.h +++ b/RecoPixelVertexing/PixelTriplets/plugins/CACell.h @@ -10,286 +10,291 @@ #include #include -class CACell { -public: - using Hit = RecHitsSortedInPhi::Hit; - using CAntuplet = std::vector; - - CACell(const HitDoublets* doublets, int doubletId, const unsigned int cellId, const int innerHitId, const int outerHitId) : - theCAState(0), theInnerHitId(innerHitId), theOuterHitId(outerHitId), theCellId(cellId), hasSameStateNeighbors(0), theDoublets(doublets), theDoubletId(doubletId), - theInnerR(doublets->r(doubletId, HitDoublets::inner)), theOuterR(doublets->r(doubletId, HitDoublets::outer)), - theInnerZ(doublets->z(doubletId, HitDoublets::inner)), theOuterZ(doublets->z(doubletId, HitDoublets::outer)) { - } - - unsigned int getCellId() const { - return theCellId; - } - - Hit const & getInnerHit() const { - return theDoublets->hit(theDoubletId, HitDoublets::inner); - } - - Hit const & getOuterHit() const { - return theDoublets->hit(theDoubletId, HitDoublets::outer); - } - - float getInnerX() const { - return theDoublets->x(theDoubletId, HitDoublets::inner); - } - - float getOuterX() const { - return theDoublets->x(theDoubletId, HitDoublets::outer); - } - - float getInnerY() const { - return theDoublets->y(theDoubletId, HitDoublets::inner); - } - - float getOuterY() const { - return theDoublets->y(theDoubletId, HitDoublets::outer); - } - - float getInnerZ() const { - return theInnerZ; - } - - float getOuterZ() const { - return theOuterZ; - } - - float getInnerR() const { - return theInnerR; - } - - float getOuterR() const { - return theOuterR; - } - - float getInnerPhi() const { - return theDoublets->phi(theDoubletId, HitDoublets::inner); - } - - float getOuterPhi() const { - return theDoublets->phi(theDoubletId, HitDoublets::outer); - } - - unsigned int getInnerHitId() const { - return theInnerHitId; - } - - unsigned int getOuterHitId() const { - return theOuterHitId; - } - - void evolve() { - - hasSameStateNeighbors = 0; - unsigned int numberOfNeighbors = theOuterNeighbors.size(); - - for (unsigned int i = 0; i < numberOfNeighbors; ++i) { - - if (theOuterNeighbors[i]->getCAState() == theCAState) { - - hasSameStateNeighbors = 1; - - break; - } - } - - } - - - void checkAlignmentAndTag(CACell* innerCell, const float ptmin, const float region_origin_x, - const float region_origin_y, const float region_origin_radius, const float thetaCut, - const float phiCut, const float hardPtCut) { - - if (areAlignedRZ(innerCell, ptmin, thetaCut) && - haveSimilarCurvature(innerCell,ptmin, region_origin_x, region_origin_y, - region_origin_radius, phiCut, hardPtCut)) { - tagAsInnerNeighbor(innerCell); - innerCell->tagAsOuterNeighbor(this); - } - } - - void checkAlignmentAndPushTriplet(CACell* innerCell, std::vector& foundTriplets, - const float ptmin, const float region_origin_x, const float region_origin_y, - const float region_origin_radius, const float thetaCut, const float phiCut, - const float hardPtCut) { - - if (areAlignedRZ(innerCell, ptmin, thetaCut) && - haveSimilarCurvature(innerCell,ptmin, region_origin_x, region_origin_y, - region_origin_radius, phiCut, hardPtCut)) { - foundTriplets.emplace_back(CACell::CAntuplet{innerCell,this}); - - } - } - - - - bool areAlignedRZ(const CACell* otherCell, const float ptmin, const float thetaCut) const - { - - float r1 = otherCell->getInnerR(); - float z1 = otherCell->getInnerZ(); - float radius_diff = std::abs(r1 - theOuterR); - float distance_13_squared = radius_diff*radius_diff + (z1 - theOuterZ)*(z1 - theOuterZ); - - float pMin = ptmin*std::sqrt(distance_13_squared); //this needs to be divided by radius_diff later - - float tan_12_13_half_mul_distance_13_squared = fabs(z1 * (theInnerR - theOuterR) + theInnerZ * (theOuterR - r1) + theOuterZ * (r1 - theInnerR)) ; - return tan_12_13_half_mul_distance_13_squared * pMin <= thetaCut * distance_13_squared * radius_diff; - } - - void tagAsOuterNeighbor(CACell* otherCell) - { - theOuterNeighbors.push_back(otherCell); - } - - void tagAsInnerNeighbor(CACell* otherCell) - { - theInnerNeighbors.push_back(otherCell); - } - - bool haveSimilarCurvature(const CACell* otherCell, const float ptmin, - const float region_origin_x, const float region_origin_y, const float region_origin_radius, const float phiCut, const float hardPtCut) const - { - - - auto x1 = otherCell->getInnerX(); - auto y1 = otherCell->getInnerY(); - - auto x2 = getInnerX(); - auto y2 = getInnerY(); - - auto x3 = getOuterX(); - auto y3 = getOuterY(); - - float distance_13_squared = (x1 - x3)*(x1 - x3) + (y1 - y3)*(y1 - y3); - float tan_12_13_half_mul_distance_13_squared = std::abs(y1 * (x2 - x3) + y2 * (x3 - x1) + y3 * (x1 - x2)) ; - if(tan_12_13_half_mul_distance_13_squared * ptmin <= 1.0e-4f*distance_13_squared) - { - - float distance_3_beamspot_squared = (x3-region_origin_x) * (x3-region_origin_x) + (y3-region_origin_y) * (y3-region_origin_y); - - float dot_bs3_13 = ((x1 - x3)*( region_origin_x - x3) + (y1 - y3) * (region_origin_y-y3)); - float proj_bs3_on_13_squared = dot_bs3_13*dot_bs3_13/distance_13_squared; +class CACellStatus { - float distance_13_beamspot_squared = distance_3_beamspot_squared - proj_bs3_on_13_squared; - - if(distance_13_beamspot_squared > (region_origin_radius+phiCut)*(region_origin_radius+phiCut) ) - { - return false; - } - return true; - } - - //87 cm/GeV = 1/(3.8T * 0.3) - - //take less than radius given by the hardPtCut and reject everything below - float minRadius = hardPtCut*87.f; - - auto det = (x1 - x2) * (y2 - y3) - (x2 - x3) * (y1 - y2); - - - auto offset = x2 * x2 + y2*y2; - - auto bc = (x1 * x1 + y1 * y1 - offset)*0.5f; - - auto cd = (offset - x3 * x3 - y3 * y3)*0.5f; - - - - auto idet = 1.f / det; - - auto x_center = (bc * (y2 - y3) - cd * (y1 - y2)) * idet; - auto y_center = (cd * (x1 - x2) - bc * (x2 - x3)) * idet; - - auto radius = std::sqrt((x2 - x_center)*(x2 - x_center) + (y2 - y_center)*(y2 - y_center)); - - if(radius < minRadius) - return false; - auto centers_distance_squared = (x_center - region_origin_x)*(x_center - region_origin_x) + (y_center - region_origin_y)*(y_center - region_origin_y); - auto region_origin_radius_plus_tolerance = region_origin_radius + phiCut; - auto minimumOfIntersectionRange = (radius - region_origin_radius_plus_tolerance)*(radius - region_origin_radius_plus_tolerance); - - if (centers_distance_squared >= minimumOfIntersectionRange) { - auto maximumOfIntersectionRange = (radius + region_origin_radius_plus_tolerance)*(radius + region_origin_radius_plus_tolerance); - return centers_distance_squared <= maximumOfIntersectionRange; - } - - return false; - - - - } - - unsigned int getCAState() const { - return theCAState; - } - // if there is at least one left neighbor with the same state (friend), the state has to be increased by 1. - - void updateState() { - theCAState += hasSameStateNeighbors; - } - - bool isRootCell(const unsigned int minimumCAState) const { - return (theCAState >= minimumCAState); - } - - // trying to free the track building process from hardcoded layers, leaving the visit of the graph - // based on the neighborhood connections between cells. - - void findNtuplets(std::vector& foundNtuplets, CAntuplet& tmpNtuplet, const unsigned int minHitsPerNtuplet) const { - - // the building process for a track ends if: - // it has no outer neighbor - // it has no compatible neighbor - // the ntuplets is then saved if the number of hits it contains is greater than a threshold - - if (tmpNtuplet.size() == minHitsPerNtuplet - 1) - { - foundNtuplets.push_back(tmpNtuplet); - } - else - { - if(theOuterNeighbors.size() == 0) - { - return; - } - else - { - unsigned int numberOfOuterNeighbors = theOuterNeighbors.size(); - for (unsigned int i = 0; i < numberOfOuterNeighbors; ++i) { - tmpNtuplet.push_back((theOuterNeighbors[i])); - theOuterNeighbors[i]->findNtuplets(foundNtuplets, tmpNtuplet, minHitsPerNtuplet); - tmpNtuplet.pop_back(); - } - } - - } +public: + + unsigned char getCAState() const { + return theCAState; + } + + // if there is at least one left neighbor with the same state (friend), the state has to be increased by 1. + void updateState() { + theCAState += hasSameStateNeighbors; + } + + bool isRootCell(const unsigned int minimumCAState) const { + return (theCAState >= minimumCAState); + } + + public: + unsigned char theCAState=0; + unsigned char hasSameStateNeighbors=0; + +}; +class CACell { +public: + using Hit = RecHitsSortedInPhi::Hit; + using CAntuple = std::vector; + using CAntuplet = std::vector; + using CAColl = std::vector; + using CAStatusColl = std::vector; + + + CACell(const HitDoublets* doublets, int doubletId, const int innerHitId, const int outerHitId) : + theDoublets(doublets), theDoubletId(doubletId) + ,theInnerR(doublets->rv(doubletId, HitDoublets::inner)) + ,theInnerZ(doublets->z(doubletId, HitDoublets::inner)) + {} + + + + Hit const & getInnerHit() const { + return theDoublets->hit(theDoubletId, HitDoublets::inner); + } + + Hit const & getOuterHit() const { + return theDoublets->hit(theDoubletId, HitDoublets::outer); + } + + + float getInnerX() const { + return theDoublets->x(theDoubletId, HitDoublets::inner); + } + + float getOuterX() const { + return theDoublets->x(theDoubletId, HitDoublets::outer); + } + + float getInnerY() const { + return theDoublets->y(theDoubletId, HitDoublets::inner); + } + + float getOuterY() const { + return theDoublets->y(theDoubletId, HitDoublets::outer); + } + + float getInnerZ() const { + return theInnerZ; + } + + float getOuterZ() const { + return theDoublets->z(theDoubletId, HitDoublets::outer); + } + + float getInnerR() const { + return theInnerR; + } + + float getOuterR() const { + return theDoublets->rv(theDoubletId, HitDoublets::outer); + } + + float getInnerPhi() const { + return theDoublets->phi(theDoubletId, HitDoublets::inner); + } + + float getOuterPhi() const { + return theDoublets->phi(theDoubletId, HitDoublets::outer); + } + + void evolve(unsigned int me, CAStatusColl& allStatus) { + + allStatus[me].hasSameStateNeighbors = 0; + auto mystate = allStatus[me].theCAState; + + for (auto oc : theOuterNeighbors ) { + + if (allStatus[oc].getCAState() == mystate) { + + allStatus[me].hasSameStateNeighbors = 1; + + break; + } } + } + + + void checkAlignmentAndAct(CAColl& allCells, CAntuple & innerCells, const float ptmin, const float region_origin_x, + const float region_origin_y, const float region_origin_radius, const float thetaCut, + const float phiCut, const float hardPtCut, std::vector * foundTriplets) { + int ncells = innerCells.size(); + int constexpr VSIZE = 16; + int ok[VSIZE]; + float r1[VSIZE]; + float z1[VSIZE]; + auto ro = getOuterR(); + auto zo = getOuterZ(); + unsigned int cellId = this - &allCells.front(); + auto loop = [&](int i, int vs) { + for (int j=0;jemplace_back(CACell::CAntuplet{koc,cellId}); + else { + oc.tagAsOuterNeighbor(cellId); + } + } + } + }; + auto lim = VSIZE*(ncells/VSIZE); + for (int i=0; i& foundTriplets, + const float ptmin, const float region_origin_x, const float region_origin_y, + const float region_origin_radius, const float thetaCut, const float phiCut, + const float hardPtCut) { + checkAlignmentAndAct(allCells, innerCells, ptmin, region_origin_x, region_origin_y, region_origin_radius, thetaCut, + phiCut, hardPtCut, &foundTriplets); + } + + + int areAlignedRZ(float r1, float z1, float ro, float zo, const float ptmin, const float thetaCut) const + { + float radius_diff = std::abs(r1 - ro); + float distance_13_squared = radius_diff*radius_diff + (z1 - zo)*(z1 - zo); + + float pMin = ptmin*std::sqrt(distance_13_squared); //this needs to be divided by radius_diff later + + float tan_12_13_half_mul_distance_13_squared = fabs(z1 * (getInnerR() - ro) + getInnerZ() * (ro - r1) + zo * (r1 - getInnerR())) ; + return tan_12_13_half_mul_distance_13_squared * pMin <= thetaCut * distance_13_squared * radius_diff; + } + + + void tagAsOuterNeighbor(unsigned int otherCell) + { + theOuterNeighbors.push_back(otherCell); + } + + + bool haveSimilarCurvature(const CACell & otherCell, const float ptmin, + const float region_origin_x, const float region_origin_y, const float region_origin_radius, const float phiCut, const float hardPtCut) const + { + + + auto x1 = otherCell.getInnerX(); + auto y1 = otherCell.getInnerY(); + + auto x2 = getInnerX(); + auto y2 = getInnerY(); + + auto x3 = getOuterX(); + auto y3 = getOuterY(); + + float distance_13_squared = (x1 - x3)*(x1 - x3) + (y1 - y3)*(y1 - y3); + float tan_12_13_half_mul_distance_13_squared = std::abs(y1 * (x2 - x3) + y2 * (x3 - x1) + y3 * (x1 - x2)) ; + // high pt : just straight + if(tan_12_13_half_mul_distance_13_squared * ptmin <= 1.0e-4f*distance_13_squared) + { + + float distance_3_beamspot_squared = (x3-region_origin_x) * (x3-region_origin_x) + (y3-region_origin_y) * (y3-region_origin_y); + + float dot_bs3_13 = ((x1 - x3)*( region_origin_x - x3) + (y1 - y3) * (region_origin_y-y3)); + float proj_bs3_on_13_squared = dot_bs3_13*dot_bs3_13/distance_13_squared; + + float distance_13_beamspot_squared = distance_3_beamspot_squared - proj_bs3_on_13_squared; + + return distance_13_beamspot_squared < (region_origin_radius+phiCut)*(region_origin_radius+phiCut); + } + + //87 cm/GeV = 1/(3.8T * 0.3) + + //take less than radius given by the hardPtCut and reject everything below + float minRadius = hardPtCut*87.f; // FIXME move out and use real MagField + + auto det = (x1 - x2) * (y2 - y3) - (x2 - x3) * (y1 - y2); + + + auto offset = x2 * x2 + y2*y2; + + auto bc = (x1 * x1 + y1 * y1 - offset)*0.5f; + + auto cd = (offset - x3 * x3 - y3 * y3)*0.5f; + + + auto idet = 1.f / det; + + auto x_center = (bc * (y2 - y3) - cd * (y1 - y2)) * idet; + auto y_center = (cd * (x1 - x2) - bc * (x2 - x3)) * idet; + + auto radius = std::sqrt((x2 - x_center)*(x2 - x_center) + (y2 - y_center)*(y2 - y_center)); + + if(radius < minRadius) return false; // hard cut on pt + + auto centers_distance_squared = (x_center - region_origin_x)*(x_center - region_origin_x) + (y_center - region_origin_y)*(y_center - region_origin_y); + auto region_origin_radius_plus_tolerance = region_origin_radius + phiCut; + auto minimumOfIntersectionRange = (radius - region_origin_radius_plus_tolerance)*(radius - region_origin_radius_plus_tolerance); + + if (centers_distance_squared >= minimumOfIntersectionRange) { + auto maximumOfIntersectionRange = (radius + region_origin_radius_plus_tolerance)*(radius + region_origin_radius_plus_tolerance); + return centers_distance_squared <= maximumOfIntersectionRange; + } + + return false; + + } + + + // trying to free the track building process from hardcoded layers, leaving the visit of the graph + // based on the neighborhood connections between cells. + + void findNtuplets(CAColl& allCells, std::vector& foundNtuplets, CAntuplet& tmpNtuplet, const unsigned int minHitsPerNtuplet) const { + + // the building process for a track ends if: + // it has no outer neighbor + // it has no compatible neighbor + // the ntuplets is then saved if the number of hits it contains is greater than a threshold + + if (tmpNtuplet.size() == minHitsPerNtuplet - 1) + { + foundNtuplets.push_back(tmpNtuplet); + } + else + { + unsigned int numberOfOuterNeighbors = theOuterNeighbors.size(); + for (unsigned int i = 0; i < numberOfOuterNeighbors; ++i) { + tmpNtuplet.push_back((theOuterNeighbors[i])); + allCells[theOuterNeighbors[i]].findNtuplets(allCells,foundNtuplets, tmpNtuplet, minHitsPerNtuplet); + tmpNtuplet.pop_back(); + } + } + + } + + private: - std::vector theInnerNeighbors; - std::vector theOuterNeighbors; - - unsigned int theCAState; - - const unsigned int theInnerHitId; - const unsigned int theOuterHitId; - const unsigned int theCellId; - unsigned int hasSameStateNeighbors; - - const HitDoublets* theDoublets; - const int theDoubletId; - - const float theInnerR; - const float theOuterR; - const float theInnerZ; - const float theOuterZ; - + + CAntuple theOuterNeighbors; + + const HitDoublets* theDoublets; + const int theDoubletId; + + const float theInnerR; + const float theInnerZ; + }; diff --git a/RecoPixelVertexing/PixelTriplets/plugins/CAGraph.h b/RecoPixelVertexing/PixelTriplets/plugins/CAGraph.h index 610ca943321c9..09658241d0f60 100644 --- a/RecoPixelVertexing/PixelTriplets/plugins/CAGraph.h +++ b/RecoPixelVertexing/PixelTriplets/plugins/CAGraph.h @@ -4,9 +4,6 @@ #include #include #include -#include -#include -#include "CACell.h" struct CALayer { @@ -31,7 +28,7 @@ struct CALayer std::vector theOuterLayers; std::vector theInnerLayers; - std::vector< std::vector > isOuterHitOfCell; + std::vector< std::vector > isOuterHitOfCell; private: @@ -58,8 +55,7 @@ struct CALayerPair } std::array theLayers; - std::vector theFoundCells; - + std::array theFoundCells= {{0,0}}; }; struct CAGraph @@ -67,7 +63,6 @@ struct CAGraph std::vector theLayers; std::vector theLayerPairs; std::vector theRootLayers; - }; #endif /* CAGRAPH_H_ */ diff --git a/RecoPixelVertexing/PixelTriplets/plugins/CAHitQuadrupletGenerator.cc b/RecoPixelVertexing/PixelTriplets/plugins/CAHitQuadrupletGenerator.cc index 078483c8d7965..496dbeba9ea3c 100644 --- a/RecoPixelVertexing/PixelTriplets/plugins/CAHitQuadrupletGenerator.cc +++ b/RecoPixelVertexing/PixelTriplets/plugins/CAHitQuadrupletGenerator.cc @@ -138,7 +138,7 @@ namespace { g.theLayers[i].theInnerLayerPairs.clear(); g.theLayers[i].theOuterLayers.clear(); g.theLayers[i].theOuterLayerPairs.clear(); - + for (auto & v : g.theLayers[i].isOuterHitOfCell) v.clear(); } } @@ -234,6 +234,7 @@ void CAHitQuadrupletGenerator::hitQuadruplets(const TrackingRegion& region, hitQuadruplets(region, result, hitDoubletsPtr, g, es); theLayerCache.clear(); } + void CAHitQuadrupletGenerator::hitNtuplets(const IntermediateHitDoublets& regionDoublets, std::vector& result, const edm::EventSetup& es, @@ -292,7 +293,8 @@ void CAHitQuadrupletGenerator::hitNtuplets(const IntermediateHitDoublets& region ca.findNtuplets(foundQuadruplets, numberOfHitsInNtuplet); - + auto & allCells = ca.getAllCells(); + const QuantityDependsPtEval maxChi2Eval = maxChi2.evaluator(es); // re-used thoughout @@ -321,13 +323,13 @@ void CAHitQuadrupletGenerator::hitNtuplets(const IntermediateHitDoublets& region }; for(unsigned int i = 0; i< 3; ++i) { - auto const& ahit = foundQuadruplets[quadId][i]->getInnerHit(); + auto const& ahit = allCells[foundQuadruplets[quadId][i]].getInnerHit(); gps[i] = ahit->globalPosition(); ges[i] = ahit->globalPositionError(); barrels[i] = isBarrel(ahit->geographicalId().subdetId()); } - auto const& ahit = foundQuadruplets[quadId][2]->getOuterHit(); + auto const& ahit = allCells[foundQuadruplets[quadId][2]].getOuterHit(); gps[3] = ahit->globalPosition(); ges[3] = ahit->globalPositionError(); barrels[3] = isBarrel(ahit->geographicalId().subdetId()); @@ -336,10 +338,10 @@ void CAHitQuadrupletGenerator::hitNtuplets(const IntermediateHitDoublets& region const auto fourthLayerId = tTopo->layer(ahit->geographicalId()); const auto sideId = tTopo->side(ahit->geographicalId()); const auto subDetId = ahit->geographicalId().subdetId(); - const auto isTheSameTriplet = (quadId != 0) && (foundQuadruplets[quadId][0]->getCellId() == previousCellIds[0]) && (foundQuadruplets[quadId][1]->getCellId() == previousCellIds[1]); + const auto isTheSameTriplet = (quadId != 0) && (foundQuadruplets[quadId][0] == previousCellIds[0]) && (foundQuadruplets[quadId][1] == previousCellIds[1]); const auto isTheSameFourthLayer = (quadId != 0) && (fourthLayerId == previousfourthLayerId) && (subDetId == previousSubDetId) && (sideId == previousSideId); - previousCellIds = {{foundQuadruplets[quadId][0]->getCellId(), foundQuadruplets[quadId][1]->getCellId()}}; + previousCellIds = {{foundQuadruplets[quadId][0], foundQuadruplets[quadId][1]}}; previousfourthLayerId = fourthLayerId; @@ -362,7 +364,9 @@ void CAHitQuadrupletGenerator::hitNtuplets(const IntermediateHitDoublets& region const float thisMaxChi2 = maxChi2Eval.value(abscurv); if (theComparitor) { - SeedingHitSet tmpTriplet(foundQuadruplets[quadId][0]->getInnerHit(), foundQuadruplets[quadId][2]->getInnerHit(), foundQuadruplets[quadId][2]->getOuterHit()); + SeedingHitSet tmpTriplet(allCells[foundQuadruplets[quadId][0]].getInnerHit(), + allCells[foundQuadruplets[quadId][2]].getInnerHit(), + allCells[foundQuadruplets[quadId][2]].getOuterHit()); if (!theComparitor->compatible(tmpTriplet) ) @@ -417,13 +421,19 @@ void CAHitQuadrupletGenerator::hitNtuplets(const IntermediateHitDoublets& region { result[index].pop_back(); } - result[index].emplace_back(foundQuadruplets[quadId][0]->getInnerHit(), foundQuadruplets[quadId][1]->getInnerHit(),foundQuadruplets[quadId][2]->getInnerHit(), foundQuadruplets[quadId][2]->getOuterHit()); + result[index].emplace_back(allCells[foundQuadruplets[quadId][0]].getInnerHit(), + allCells[foundQuadruplets[quadId][1]].getInnerHit(), + allCells[foundQuadruplets[quadId][2]].getInnerHit(), + allCells[foundQuadruplets[quadId][2]].getOuterHit()); hasAlreadyPushedACandidate = true; } } else { - result[index].emplace_back(foundQuadruplets[quadId][0]->getInnerHit(), foundQuadruplets[quadId][1]->getInnerHit(), foundQuadruplets[quadId][2]->getInnerHit(), foundQuadruplets[quadId][2]->getOuterHit()); + result[index].emplace_back(allCells[foundQuadruplets[quadId][0]].getInnerHit(), + allCells[foundQuadruplets[quadId][1]].getInnerHit(), + allCells[foundQuadruplets[quadId][2]].getInnerHit(), + allCells[foundQuadruplets[quadId][2]].getOuterHit()); } } index++; @@ -433,7 +443,7 @@ void CAHitQuadrupletGenerator::hitNtuplets(const IntermediateHitDoublets& region void CAHitQuadrupletGenerator::hitQuadruplets(const TrackingRegion& region, OrderedHitSeeds & result, - std::vector& hitDoublets, const CAGraph& g, + std::vector& hitDoublets, CAGraph& g, const edm::EventSetup& es) { //Retrieve tracker topology from geometry edm::ESHandle tTopoHand; @@ -452,6 +462,7 @@ void CAHitQuadrupletGenerator::hitQuadruplets(const TrackingRegion& region, ca.findNtuplets(foundQuadruplets, numberOfHitsInNtuplet); + auto & allCells = ca.getAllCells(); const QuantityDependsPtEval maxChi2Eval = maxChi2.evaluator(es); @@ -481,13 +492,13 @@ void CAHitQuadrupletGenerator::hitQuadruplets(const TrackingRegion& region, }; for(unsigned int i = 0; i< 3; ++i) { - auto const& ahit = foundQuadruplets[quadId][i]->getInnerHit(); + auto const& ahit = allCells[foundQuadruplets[quadId][i]].getInnerHit(); gps[i] = ahit->globalPosition(); ges[i] = ahit->globalPositionError(); barrels[i] = isBarrel(ahit->geographicalId().subdetId()); } - auto const& ahit = foundQuadruplets[quadId][2]->getOuterHit(); + auto const& ahit = allCells[foundQuadruplets[quadId][2]].getOuterHit(); gps[3] = ahit->globalPosition(); ges[3] = ahit->globalPositionError(); barrels[3] = isBarrel(ahit->geographicalId().subdetId()); @@ -497,10 +508,10 @@ void CAHitQuadrupletGenerator::hitQuadruplets(const TrackingRegion& region, const auto fourthLayerId = tTopo->layer(ahit->geographicalId()); const auto sideId = tTopo->side(ahit->geographicalId()); const auto subDetId = ahit->geographicalId().subdetId(); - const auto isTheSameTriplet = (quadId != 0) && (foundQuadruplets[quadId][0]->getCellId() == previousCellIds[0]) && (foundQuadruplets[quadId][1]->getCellId() == previousCellIds[1]); + const auto isTheSameTriplet = (quadId != 0) && (foundQuadruplets[quadId][0] == previousCellIds[0]) && (foundQuadruplets[quadId][1] == previousCellIds[1]); const auto isTheSameFourthLayer = (quadId != 0) && (fourthLayerId == previousfourthLayerId) && (subDetId == previousSubDetId) && (sideId == previousSideId); - previousCellIds = {{foundQuadruplets[quadId][0]->getCellId(), foundQuadruplets[quadId][1]->getCellId()}}; + previousCellIds = {{foundQuadruplets[quadId][0], foundQuadruplets[quadId][1]}}; previousfourthLayerId = fourthLayerId; @@ -526,7 +537,9 @@ void CAHitQuadrupletGenerator::hitQuadruplets(const TrackingRegion& region, if (theComparitor) { - SeedingHitSet tmpTriplet(foundQuadruplets[quadId][0]->getInnerHit(), foundQuadruplets[quadId][2]->getInnerHit(), foundQuadruplets[quadId][2]->getOuterHit()); + SeedingHitSet tmpTriplet(allCells[foundQuadruplets[quadId][0]].getInnerHit(), + allCells[foundQuadruplets[quadId][2]].getInnerHit(), + allCells[foundQuadruplets[quadId][2]].getOuterHit()); if (!theComparitor->compatible(tmpTriplet) ) @@ -587,16 +600,21 @@ void CAHitQuadrupletGenerator::hitQuadruplets(const TrackingRegion& region, result.pop_back(); } - result.emplace_back(foundQuadruplets[quadId][0]->getInnerHit(), foundQuadruplets[quadId][1]->getInnerHit(), - foundQuadruplets[quadId][2]->getInnerHit(), foundQuadruplets[quadId][2]->getOuterHit()); + + result.emplace_back(allCells[foundQuadruplets[quadId][0]].getInnerHit(), + allCells[foundQuadruplets[quadId][1]].getInnerHit(), + allCells[foundQuadruplets[quadId][2]].getInnerHit(), + allCells[foundQuadruplets[quadId][2]].getOuterHit()); hasAlreadyPushedACandidate = true; - + } } else { - - result.emplace_back(foundQuadruplets[quadId][0]->getInnerHit(), foundQuadruplets[quadId][1]->getInnerHit(), foundQuadruplets[quadId][2]->getInnerHit(), foundQuadruplets[quadId][2]->getOuterHit()); + result.emplace_back(allCells[foundQuadruplets[quadId][0]].getInnerHit(), + allCells[foundQuadruplets[quadId][1]].getInnerHit(), + allCells[foundQuadruplets[quadId][2]].getInnerHit(), + allCells[foundQuadruplets[quadId][2]].getOuterHit()); } } diff --git a/RecoPixelVertexing/PixelTriplets/plugins/CAHitQuadrupletGenerator.h b/RecoPixelVertexing/PixelTriplets/plugins/CAHitQuadrupletGenerator.h index c98feb1709db9..e3878b92cb648 100644 --- a/RecoPixelVertexing/PixelTriplets/plugins/CAHitQuadrupletGenerator.h +++ b/RecoPixelVertexing/PixelTriplets/plugins/CAHitQuadrupletGenerator.h @@ -63,7 +63,7 @@ class CAHitQuadrupletGenerator : public HitQuadrupletGenerator { // actual work void hitQuadruplets(const TrackingRegion& reg, OrderedHitSeeds& result, std::vector& hitDoublets, - const CAGraph& g, + CAGraph& g, const edm::EventSetup& es); edm::EDGetTokenT theSeedingLayerToken; diff --git a/RecoPixelVertexing/PixelTriplets/plugins/CAHitTripletGenerator.cc b/RecoPixelVertexing/PixelTriplets/plugins/CAHitTripletGenerator.cc index 4e9b961025f4b..8100395b625b6 100644 --- a/RecoPixelVertexing/PixelTriplets/plugins/CAHitTripletGenerator.cc +++ b/RecoPixelVertexing/PixelTriplets/plugins/CAHitTripletGenerator.cc @@ -132,7 +132,7 @@ namespace { g.theLayers[i].theInnerLayerPairs.clear(); g.theLayers[i].theOuterLayers.clear(); g.theLayers[i].theOuterLayerPairs.clear(); - + for (auto & v : g.theLayers[i].isOuterHitOfCell) v.clear(); } } @@ -273,7 +273,7 @@ void CAHitTripletGenerator::hitNtuplets(const IntermediateHitDoublets& regionDou ca.findTriplets(hitDoublets, foundTriplets, region, caThetaCut, caPhiCut, caHardPtCut); - + auto & allCells = ca.getAllCells(); const QuantityDependsPtEval maxChi2Eval = maxChi2.evaluator(es); @@ -290,9 +290,9 @@ void CAHitTripletGenerator::hitNtuplets(const IntermediateHitDoublets& regionDou ++tripletId) { - OrderedHitTriplet tmpTriplet(foundTriplets[tripletId][0]->getInnerHit(), - foundTriplets[tripletId][0]->getOuterHit(), - foundTriplets[tripletId][1]->getOuterHit()); + OrderedHitTriplet tmpTriplet(allCells[foundTriplets[tripletId][0]].getInnerHit(), + allCells[foundTriplets[tripletId][0]].getOuterHit(), + allCells[foundTriplets[tripletId][1]].getOuterHit()); auto isBarrel = [](const unsigned id) -> bool { @@ -300,13 +300,13 @@ void CAHitTripletGenerator::hitNtuplets(const IntermediateHitDoublets& regionDou }; for (unsigned int i = 0; i < 2; ++i) { - auto const& ahit = foundTriplets[tripletId][i]->getInnerHit(); + auto const& ahit = allCells[foundTriplets[tripletId][i]].getInnerHit(); gps[i] = ahit->globalPosition(); ges[i] = ahit->globalPositionError(); barrels[i] = isBarrel(ahit->geographicalId().subdetId()); } - auto const& ahit = foundTriplets[tripletId][1]->getOuterHit(); + auto const& ahit = allCells[foundTriplets[tripletId][1]].getOuterHit(); gps[2] = ahit->globalPosition(); ges[2] = ahit->globalPositionError(); barrels[2] = isBarrel(ahit->geographicalId().subdetId()); @@ -374,13 +374,15 @@ void CAHitTripletGenerator::hitNtuplets(const IntermediateHitDoublets& regionDou void CAHitTripletGenerator::hitTriplets(const TrackingRegion& region, OrderedHitTriplets & result, - std::vector& hitDoublets, const CAGraph& g, + std::vector& hitDoublets, CAGraph& g, const edm::EventSetup& es) { std::vector foundTriplets; CellularAutomaton ca(g); ca.findTriplets(hitDoublets, foundTriplets, region, caThetaCut, caPhiCut, caHardPtCut); + + auto & allCells = ca.getAllCells(); unsigned int numberOfFoundTriplets = foundTriplets.size(); @@ -398,9 +400,9 @@ void CAHitTripletGenerator::hitTriplets(const TrackingRegion& region, ++tripletId) { - OrderedHitTriplet tmpTriplet(foundTriplets[tripletId][0]->getInnerHit(), - foundTriplets[tripletId][0]->getOuterHit(), - foundTriplets[tripletId][1]->getOuterHit()); + OrderedHitTriplet tmpTriplet(allCells[foundTriplets[tripletId][0]].getInnerHit(), + allCells[foundTriplets[tripletId][0]].getOuterHit(), + allCells[foundTriplets[tripletId][1]].getOuterHit()); auto isBarrel = [](const unsigned id) -> bool { @@ -408,13 +410,13 @@ void CAHitTripletGenerator::hitTriplets(const TrackingRegion& region, }; for (unsigned int i = 0; i < 2; ++i) { - auto const& ahit = foundTriplets[tripletId][i]->getInnerHit(); + auto const& ahit = allCells[foundTriplets[tripletId][i]].getInnerHit(); gps[i] = ahit->globalPosition(); ges[i] = ahit->globalPositionError(); barrels[i] = isBarrel(ahit->geographicalId().subdetId()); } - auto const& ahit = foundTriplets[tripletId][1]->getOuterHit(); + auto const& ahit = allCells[foundTriplets[tripletId][1]].getOuterHit(); gps[2] = ahit->globalPosition(); ges[2] = ahit->globalPositionError(); barrels[2] = isBarrel(ahit->geographicalId().subdetId()); diff --git a/RecoPixelVertexing/PixelTriplets/plugins/CAHitTripletGenerator.h b/RecoPixelVertexing/PixelTriplets/plugins/CAHitTripletGenerator.h index e71acec54cad3..4754a4e9d3ffd 100644 --- a/RecoPixelVertexing/PixelTriplets/plugins/CAHitTripletGenerator.h +++ b/RecoPixelVertexing/PixelTriplets/plugins/CAHitTripletGenerator.h @@ -60,7 +60,7 @@ class CAHitTripletGenerator : public HitTripletGenerator { // actual work void hitTriplets(const TrackingRegion& reg, OrderedHitTriplets& result, std::vector& hitDoublets, - const CAGraph& g, + CAGraph& g, const edm::EventSetup& es); edm::EDGetTokenT theSeedingLayerToken; diff --git a/RecoPixelVertexing/PixelTriplets/plugins/CellularAutomaton.cc b/RecoPixelVertexing/PixelTriplets/plugins/CellularAutomaton.cc index 7a3f7aef73cb1..d76a6793b36ae 100644 --- a/RecoPixelVertexing/PixelTriplets/plugins/CellularAutomaton.cc +++ b/RecoPixelVertexing/PixelTriplets/plugins/CellularAutomaton.cc @@ -1,9 +1,15 @@ #include "CellularAutomaton.h" +#include + + void CellularAutomaton::createAndConnectCells(const std::vector& hitDoublets, const TrackingRegion& region, const float thetaCut, const float phiCut, const float hardPtCut) { - unsigned int cellId = 0; + int tsize=0; + for ( auto hd : hitDoublets) tsize+=hd->size(); + allCells.reserve(tsize); + unsigned int cellId = 0; float ptmin = region.ptMin(); float region_origin_x = region.origin().x(); float region_origin_y = region.origin().y(); @@ -50,27 +56,29 @@ void CellularAutomaton::createAndConnectCells(const std::vectorsize(); - currentLayerPairRef.theFoundCells.reserve(numberOfDoublets); + currentLayerPairRef.theFoundCells[0] = cellId; + currentLayerPairRef.theFoundCells[1] = cellId+numberOfDoublets; for (unsigned int i = 0; i < numberOfDoublets; ++i) { - currentLayerPairRef.theFoundCells.emplace_back( - doubletLayerPairId, i, cellId, + allCells.emplace_back(doubletLayerPairId, i, doubletLayerPairId->innerHitId(i), doubletLayerPairId->outerHitId(i)); - currentOuterLayerRef.isOuterHitOfCell[doubletLayerPairId->outerHitId(i)].push_back( - &(currentLayerPairRef.theFoundCells[i])); - cellId++; - - for (auto neigCell : currentInnerLayerRef.isOuterHitOfCell[doubletLayerPairId->innerHitId(i)]) - { - currentLayerPairRef.theFoundCells[i].checkAlignmentAndTag( - neigCell, ptmin, region_origin_x, - region_origin_y, region_origin_radius, thetaCut, - phiCut, hardPtCut); - } + + + currentOuterLayerRef.isOuterHitOfCell[doubletLayerPairId->outerHitId(i)].push_back(cellId); + + cellId++; + + auto & neigCells = currentInnerLayerRef.isOuterHitOfCell[doubletLayerPairId->innerHitId(i)]; + allCells.back().checkAlignmentAndTag(allCells, + neigCells, ptmin, region_origin_x, + region_origin_y, region_origin_radius, thetaCut, + phiCut, hardPtCut + ); + } - + assert(cellId==currentLayerPairRef.theFoundCells[1]); for (auto outerLayerPair : currentOuterLayerRef.theOuterLayerPairs) { LayerPairsToVisit.push(outerLayerPair); @@ -93,62 +101,67 @@ void CellularAutomaton::createAndConnectCells(const std::vector& foundNtuplets, const unsigned int minHitsPerNtuplet) { - std::vector tmpNtuplet; + CACell::CAntuple tmpNtuplet; tmpNtuplet.reserve(minHitsPerNtuplet); - for (CACell* root_cell : theRootCells) + for (auto root_cell : theRootCells) { - tmpNtuplet.clear(); - tmpNtuplet.push_back(root_cell); - root_cell->findNtuplets(foundNtuplets, tmpNtuplet, minHitsPerNtuplet); + tmpNtuplet.clear(); + tmpNtuplet.push_back(root_cell); + allCells[root_cell].findNtuplets(allCells,foundNtuplets, tmpNtuplet, minHitsPerNtuplet); } } @@ -157,7 +170,11 @@ void CellularAutomaton::findNtuplets( void CellularAutomaton::findTriplets(const std::vector& hitDoublets,std::vector& foundTriplets, const TrackingRegion& region, const float thetaCut, const float phiCut, const float hardPtCut) { - unsigned int cellId = 0; + int tsize=0; + for ( auto hd : hitDoublets) tsize+=hd->size(); + allCells.reserve(tsize); + + unsigned int cellId = 0; float ptmin = region.ptMin(); float region_origin_x = region.origin().x(); float region_origin_y = region.origin().y(); @@ -204,26 +221,28 @@ void CellularAutomaton::findTriplets(const std::vector& hitD const HitDoublets* doubletLayerPairId = hitDoublets[currentLayerPair]; auto numberOfDoublets = doubletLayerPairId->size(); - currentLayerPairRef.theFoundCells.reserve(numberOfDoublets); + currentLayerPairRef.theFoundCells[0] = cellId; + currentLayerPairRef.theFoundCells[1] = cellId+numberOfDoublets; for (unsigned int i = 0; i < numberOfDoublets; ++i) { - currentLayerPairRef.theFoundCells.emplace_back( - doubletLayerPairId, i, cellId, + allCells.emplace_back(doubletLayerPairId, i, doubletLayerPairId->innerHitId(i), doubletLayerPairId->outerHitId(i)); - currentOuterLayerRef.isOuterHitOfCell[doubletLayerPairId->outerHitId(i)].push_back( - &(currentLayerPairRef.theFoundCells[i])); - cellId++; - - for (auto neigCell : currentInnerLayerRef.isOuterHitOfCell[doubletLayerPairId->innerHitId(i)]) - { - currentLayerPairRef.theFoundCells[i].checkAlignmentAndPushTriplet( - neigCell, foundTriplets, ptmin, region_origin_x, - region_origin_y, region_origin_radius, thetaCut, - phiCut, hardPtCut); - } - + + + currentOuterLayerRef.isOuterHitOfCell[doubletLayerPairId->outerHitId(i)].push_back(cellId); + + cellId++; + + auto & neigCells = currentInnerLayerRef.isOuterHitOfCell[doubletLayerPairId->innerHitId(i)]; + allCells.back().checkAlignmentAndPushTriplet(allCells, + neigCells, foundTriplets, ptmin, region_origin_x, + region_origin_y, region_origin_radius, thetaCut, + phiCut, hardPtCut + ); } + assert(cellId==currentLayerPairRef.theFoundCells[1]); + for (auto outerLayerPair : currentOuterLayerRef.theOuterLayerPairs) { diff --git a/RecoPixelVertexing/PixelTriplets/plugins/CellularAutomaton.h b/RecoPixelVertexing/PixelTriplets/plugins/CellularAutomaton.h index ec144f00985a9..2166e271cc66c 100644 --- a/RecoPixelVertexing/PixelTriplets/plugins/CellularAutomaton.h +++ b/RecoPixelVertexing/PixelTriplets/plugins/CellularAutomaton.h @@ -8,25 +8,31 @@ class CellularAutomaton { public: - CellularAutomaton(const CAGraph& graph) - : theLayerGraph(graph) - { - - } - - void createAndConnectCells(const std::vector&, - const TrackingRegion&, const float, const float, const float); - - void evolve(const unsigned int); - void findNtuplets(std::vector&, const unsigned int); - void findTriplets(const std::vector& hitDoublets,std::vector& foundTriplets, const TrackingRegion& region, - const float thetaCut, const float phiCut, const float hardPtCut); - + CellularAutomaton(CAGraph& graph) + : theLayerGraph(graph) + { + + } + + std::vector & getAllCells() { return allCells;} + + void createAndConnectCells(const std::vector&, + const TrackingRegion&, const float, const float, const float); + + void evolve(const unsigned int); + void findNtuplets(std::vector&, const unsigned int); + void findTriplets(const std::vector& hitDoublets,std::vector& foundTriplets, const TrackingRegion& region, + const float thetaCut, const float phiCut, const float hardPtCut); + private: - CAGraph theLayerGraph; - std::vector theRootCells; - std::vector > theNtuplets; + CAGraph & theLayerGraph; + + std::vector allCells; + std::vector allStatus; + std::vector theRootCells; + std::vector > theNtuplets; + }; #endif diff --git a/RecoTracker/TkHitPairs/interface/RecHitsSortedInPhi.h b/RecoTracker/TkHitPairs/interface/RecHitsSortedInPhi.h index c61736a94db21..70f7db42600b9 100644 --- a/RecoTracker/TkHitPairs/interface/RecHitsSortedInPhi.h +++ b/RecoTracker/TkHitPairs/interface/RecHitsSortedInPhi.h @@ -7,6 +7,8 @@ #include #include +#include + /** A RecHit container sorted in phi. * Provides fast access for hits in a given phi window * using binary search. @@ -85,8 +87,8 @@ class RecHitsSortedInPhi { public: float phi(int i) const { return theHits[i].phi();} - float gv(int i) const { return isBarrel ? z[i] : gp(i).perp();} // global v - float rv(int i) const { return isBarrel ? u[i] : v[i];} // dispaced r + float gv(int i) const { return isBarrel ? z[i] : gp(i).perp();} // global v + float rv(int i) const { return isBarrel ? u[i] : v[i];} // dispaced r GlobalPoint gp(int i) const { return GlobalPoint(x[i],y[i],z[i]);} public: @@ -128,33 +130,41 @@ class HitDoublets { public: enum layer { inner=0, outer=1}; + using HitLayer = RecHitsSortedInPhi; using Hit=RecHitsSortedInPhi::Hit; - - + using ADoublet = std::pair; + HitDoublets( RecHitsSortedInPhi const & in, RecHitsSortedInPhi const & out) : layers{{&in,&out}}{} HitDoublets(HitDoublets && rh) : layers(std::move(rh.layers)), indeces(std::move(rh.indeces)){} - - void reserve(std::size_t s) { indeces.reserve(2*s);} - std::size_t size() const { return indeces.size()/2;} + + void reserve(std::size_t s) { indeces.reserve(s);} + std::size_t size() const { return indeces.size();} bool empty() const { return indeces.empty();} void clear() { indeces.clear();} - void shrink_to_fit() { indeces.shrink_to_fit();} - - void add (int il, int ol) { indeces.push_back(il);indeces.push_back(ol);} + void shrink_to_fit() { + indeces.shrink_to_fit(); + } + + void add (int il, int ol) { + indeces.emplace_back(il,ol); + } + int index(int i, layer l) const { return l==inner ? innerHitId(i) : outerHitId(i);} DetLayer const * detLayer(layer l) const { return layers[l]->layer; } - int innerHitId(int i) const {return indeces[2*i];} - int outerHitId(int i) const {return indeces[2*i+1];} - Hit const & hit(int i, layer l) const { return layers[l]->theHits[indeces[2*i+l]].hit();} - float phi(int i, layer l) const { return layers[l]->phi(indeces[2*i+l]);} - float rv(int i, layer l) const { return layers[l]->rv(indeces[2*i+l]);} - float r(int i, layer l) const { float xp = x(i,l); float yp = y(i,l); return sqrt (xp*xp + yp*yp);} - float z(int i, layer l) const { return layers[l]->z[indeces[2*i+l]];} - float x(int i, layer l) const { return layers[l]->x[indeces[2*i+l]];} - float y(int i, layer l) const { return layers[l]->y[indeces[2*i+l]];} + HitLayer const & innerLayer() const { return *layers[inner];} + HitLayer const & outerLayer() const { return *layers[outer];} + int innerHitId(int i) const {return indeces[i].first;} + int outerHitId(int i) const {return indeces[i].second;} + Hit const & hit(int i, layer l) const { return layers[l]->theHits[index(i,l)].hit();} + float phi(int i, layer l) const { return layers[l]->phi(index(i,l));} + float rv(int i, layer l) const { return layers[l]->rv(index(i,l));} + float r(int i, layer l) const { float xp = x(i,l); float yp = y(i,l); return std::sqrt (xp*xp + yp*yp);} + float z(int i, layer l) const { return layers[l]->z[index(i,l)];} + float x(int i, layer l) const { return layers[l]->x[index(i,l)];} + float y(int i, layer l) const { return layers[l]->y[index(i,l)];} GlobalPoint gp(int i, layer l) const { return GlobalPoint(x(i,l),y(i,l),z(i,l));} private: @@ -162,7 +172,7 @@ class HitDoublets { std::array layers; - std::vector indeces; + std::vector indeces; // naturally sorted by outerId };