Skip to content

Commit

Permalink
Merge pull request #7896 from VinInn/Loopers
Browse files Browse the repository at this point in the history 
factorize out choice of plane, use it for loopers as well
  • Loading branch information
@cmsbuild
cmsbuild committed Feb 26, 2015
2 parents e39920c + 18b28e8 commit 764ccc0
Show file tree
Hide file tree
Showing 3 changed files with 150 additions and 31 deletions.
63 changes: 63 additions & 0 deletions TrackingTools/GeomPropagators/interface/OptimalHelixPlaneCrossing.h
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,63 @@
#ifndef OptimalHelixPlaneCrossing_H
#define OptimalHelixPlaneCrossing_H


#include "TrackingTools/GeomPropagators/interface/HelixBarrelPlaneCrossingByCircle.h"
#include "TrackingTools/GeomPropagators/interface/HelixForwardPlaneCrossing.h"
#include "TrackingTools/GeomPropagators/interface/HelixArbitraryPlaneCrossing.h"

#include<memory>


class OptimalHelixPlaneCrossing {
public:
using Base = HelixPlaneCrossing;


template<typename ... Args>
OptimalHelixPlaneCrossing(Plane const & plane, Args&&... args) {

GlobalVector u = plane.normalVector();
constexpr float small = 1.e-6; // for orientation of planes


if (std::abs(u.z()) < small) {
// barrel plane:
// instantiate HelixBarrelPlaneCrossing,
new(get()) HelixBarrelPlaneCrossingByCircle(args...);
}
else if ( (std::abs(u.x()) < small) & (std::abs(u.y()) < small) ) {
// forward plane:
// instantiate HelixForwardPlaneCrossing
new(get()) HelixForwardPlaneCrossing(args...);
} else {
// arbitrary plane:
// instantiate HelixArbitraryPlaneCrossing
new(get()) HelixArbitraryPlaneCrossing(args...);
}

}


~OptimalHelixPlaneCrossing() { get()->~Base(); }

Base & operator*() { return *get(); }
Base const & operator*() const { return *get();}



private:


Base * get() { return (Base*)&mem;}
Base const * get() const { return (Base const *)&mem;}

union Tmp { HelixBarrelPlaneCrossingByCircle a; HelixForwardPlaneCrossing b; HelixArbitraryPlaneCrossing c;};
using aligned_union_t = typename std::aligned_storage<sizeof(Tmp),alignof(Tmp)>::type;
aligned_union_t mem;

};



#endif
37 changes: 10 additions & 27 deletions TrackingTools/GeomPropagators/src/AnalyticalPropagator.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,9 +9,7 @@
#include "TrackingTools/GeomPropagators/interface/PropagationExceptions.h"
#include "TrackingTools/GeomPropagators/interface/StraightLinePlaneCrossing.h"
#include "TrackingTools/GeomPropagators/interface/StraightLineBarrelCylinderCrossing.h"
#include "TrackingTools/GeomPropagators/interface/HelixBarrelPlaneCrossingByCircle.h"
#include "TrackingTools/GeomPropagators/interface/HelixForwardPlaneCrossing.h"
#include "TrackingTools/GeomPropagators/interface/HelixArbitraryPlaneCrossing.h"
#include "TrackingTools/GeomPropagators/interface/OptimalHelixPlaneCrossing.h"
#include "TrackingTools/GeomPropagators/interface/HelixBarrelCylinderCrossing.h"
#include "TrackingTools/AnalyticalJacobians/interface/AnalyticalCurvilinearJacobian.h"
#include "TrackingTools/GeomPropagators/interface/PropagationDirectionFromPath.h"
Expand Down Expand Up @@ -207,8 +205,7 @@ AnalyticalPropagator::propagateParametersOnPlane(const FreeTrajectoryState& fts,
//
// Helix case
//
GlobalVector u = plane.normalVector();
constexpr float small = 1.e-6; // for orientation of planes

//
// Frame-independant point and vector are created explicitely to
// avoid confusing gcc (refuses to compile with temporary objects
Expand All @@ -217,22 +214,8 @@ AnalyticalPropagator::propagateParametersOnPlane(const FreeTrajectoryState& fts,
HelixPlaneCrossing::PositionType helixPos(x);
HelixPlaneCrossing::DirectionType helixDir(p);
if likely(isOldPropagationType) {
if (std::abs(u.z()) < small) {
// barrel plane:
// instantiate HelixBarrelPlaneCrossing, get vector of solutions and check for existance
HelixBarrelPlaneCrossingByCircle planeCrossing(helixPos,helixDir,rho,propagationDirection());
return propagateWithHelixCrossing(planeCrossing,plane,fts.momentum().mag(),x,p,s);
}
if ( (std::abs(u.x()) < small) & (std::abs(u.y()) < small) ) {
// forward plane:
// instantiate HelixForwardPlaneCrossing, get vector of solutions and check for existance
HelixForwardPlaneCrossing planeCrossing(helixPos,helixDir,rho,propagationDirection());
return propagateWithHelixCrossing(planeCrossing,plane,fts.momentum().mag(),x,p,s);
}
// arbitrary plane:
// instantiate HelixArbitraryPlaneCrossing, get vector of solutions and check for existance
HelixArbitraryPlaneCrossing planeCrossing(helixPos,helixDir,rho,propagationDirection());
return propagateWithHelixCrossing(planeCrossing,plane,fts.momentum().mag(),x,p,s);
OptimalHelixPlaneCrossing planeCrossing(plane,helixPos,helixDir,rho,propagationDirection());
return propagateWithHelixCrossing(*planeCrossing,plane,fts.momentum().mag(),x,p,s);
}


Expand All @@ -257,7 +240,7 @@ AnalyticalPropagator::propagateParametersOnPlane(const FreeTrajectoryState& fts,
HelixPlaneCrossing::PositionType helixPos1(gp1);
HelixPlaneCrossing::DirectionType helixDir1(gm1);
LogDebug("AnalyticalPropagator") << "gp1 before calling planeCrossing1: " << gp1 << "\n";
HelixArbitraryPlaneCrossing planeCrossing1(helixPos1,helixDir1,rho1,propagationDirection());
OptimalHelixPlaneCrossing planeCrossing1(plane,helixPos1,helixDir1,rho1,propagationDirection());

HelixPlaneCrossing::PositionType xGen;
HelixPlaneCrossing::DirectionType pGen;
Expand All @@ -266,7 +249,7 @@ AnalyticalPropagator::propagateParametersOnPlane(const FreeTrajectoryState& fts,
if(propagationDirection()==oppositeToMomentum)
tolerance *=-1;

bool check1 = propagateWithHelixCrossing(planeCrossing1,plane,fts.momentum().mag(),gp1,gm1,s1);
bool check1 = propagateWithHelixCrossing(*planeCrossing1,plane,fts.momentum().mag(),gp1,gm1,s1);
double dphi1 = fabs(fts.momentum().phi()-gm1.phi());
LogDebug("AnalyticalPropagator") << "check1, s1, dphi, gp1: "
<< check1 << " , "
Expand All @@ -276,8 +259,8 @@ AnalyticalPropagator::propagateParametersOnPlane(const FreeTrajectoryState& fts,

//move forward a bit to avoid that the propagator doesn't propagate because the state is already on surface.
//we want to go to the other point of intersection between the helix and the plane
xGen = planeCrossing1.position(s1+tolerance);
pGen = planeCrossing1.direction(s1+tolerance);
xGen = (*planeCrossing1).position(s1+tolerance);
pGen = (*planeCrossing1).direction(s1+tolerance);

/*
if(!check1 || s1>170 ){
Expand Down Expand Up @@ -306,9 +289,9 @@ AnalyticalPropagator::propagateParametersOnPlane(const FreeTrajectoryState& fts,
double rho2 = rho1;
HelixPlaneCrossing::PositionType helixPos2(gp2);
HelixPlaneCrossing::DirectionType helixDir2(gm2);
HelixArbitraryPlaneCrossing planeCrossing2(helixPos2,helixDir2,rho2,propagationDirection());
OptimalHelixPlaneCrossing planeCrossing2(plane,helixPos2,helixDir2,rho2,propagationDirection());

bool check2 = propagateWithHelixCrossing(planeCrossing2,plane,gm2.mag(),gp2,gm2,s2);
bool check2 = propagateWithHelixCrossing(*planeCrossing2,plane,gm2.mag(),gp2,gm2,s2);

if(!check2){
x = gp1;
Expand Down
81 changes: 77 additions & 4 deletions TrackingTools/GeomPropagators/test/HelixPropagators_t.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,8 @@
#include "DataFormats/GeometrySurface/interface/Plane.h"
#include "TrackingTools/GeomPropagators/interface/HelixBarrelPlaneCrossing2OrderLocal.h"
#include "TrackingTools/GeomPropagators/interface/HelixBarrelPlaneCrossingByCircle.h"
#include "TrackingTools/GeomPropagators/interface/OptimalHelixPlaneCrossing.h"


#include <algorithm>
#include <cmath>
Expand Down Expand Up @@ -67,8 +69,11 @@ void crossing1() {
-0.0969205,0.995289,0.00253912);

std::cout << rot << std::endl;

Plane plane(pos,rot);
GlobalVector u = plane.normalVector();
std::cout << "norm " << u << std::endl;


HelixBarrelPlaneCrossingByCircle precise(startingPos, startingDir, rho,alongMomentum);
bool cross; double s;
Expand Down Expand Up @@ -97,18 +102,40 @@ void crossing2() {
0.000108324,-0.00233467,0.999997,
0.0994701,-0.995038,-0.00233387);
std::cout << rot << std::endl;



Plane plane(pos,rot);

GlobalVector u = plane.normalVector();
std::cout << "norm " << u << std::endl;



double rho = 0.00223254;

HelixBarrelPlaneCrossingByCircle precise(startingPos, startingDir, rho,alongMomentum);
bool cross; double s;
std::tie(cross,s) = precise.pathLength(plane);

HelixBarrelPlaneCrossing2OrderLocal crossing(startingPos, startingDir, rho, plane);

std::cout << s << ' ' << precise.position(s) << ' ' << precise.direction(s) << std::endl;
std::cout << plane.toLocal(GlobalPoint(precise.position(s))) << " " << plane.toLocal(GlobalVector(precise.direction(s))) << std::endl;

{

OptimalHelixPlaneCrossing ocrossing(plane, HelixPlaneCrossing::PositionType(startingPos), HelixPlaneCrossing::DirectionType(startingDir), rho, alongMomentum);
auto & crossing = *ocrossing;
bool cross; double s;
std::tie(cross,s) = crossing.pathLength(plane);
if (!cross) std::cout << "missed!" << std::endl;
else {
std::cout << s << ' ' << crossing.position(s) << ' ' << crossing.direction(s) << std::endl;
std::cout << plane.toLocal(GlobalPoint(crossing.position(s))) << " "
<< plane.toLocal(GlobalVector(crossing.direction(s))) << std::endl;

}
}

HelixBarrelPlaneCrossing2OrderLocal crossing(startingPos, startingDir, rho, plane);
std::cout << HelixBarrelPlaneCrossing2OrderLocal::positionOnly(startingPos, startingDir, rho, plane) << ' ';
std::cout << crossing.position() << ' ' << crossing.direction() << std::endl;

Expand All @@ -121,12 +148,58 @@ void crossing2() {

}



void crossing3() {

GlobalPoint startingPos(-8.12604,-50.829,9.82116);
GlobalVector startingDir(-0.517536,-5.09581,1.02212);

double rho = 0.00223254;


GlobalPoint pos(26.2991,36.199,100.263);
Surface::RotationType rot(-0.808996,0.587813,1.16466e-16,
0.587813,0.808996,-3.78444e-17,
-1.16466e-16,3.78444e-17,-1);

std::cout << rot << std::endl;

Plane plane(pos,rot);
GlobalVector u = plane.normalVector();
std::cout << "norm " << u << std::endl;


OptimalHelixPlaneCrossing ocrossing(plane, HelixPlaneCrossing::PositionType(startingPos), HelixPlaneCrossing::DirectionType(startingDir), rho, alongMomentum);
auto & crossing = *ocrossing;
bool cross; double s;
std::tie(cross,s) = crossing.pathLength(plane);
if (!cross) std::cout << "missed!" << std::endl;
else {
std::cout << crossing.position(s) << ' ' << crossing.direction(s) << std::endl;
std::cout << plane.toLocal(GlobalPoint(crossing.position(s))) << " "
<< plane.toLocal(GlobalVector(crossing.direction(s))) << std::endl;

}

LocalPoint thePos; LocalVector theDir;
std::tie(thePos,theDir) = secondOrderAccurate(startingPos, startingDir, rho, plane);

std::cout << thePos << ' ' << theDir << std::endl;

}



void testHelixBarrelPlaneCrossing2OrderLocal() {

crossing1();
std::cout << std::endl;

crossing2();
std::cout << std::endl;

crossing3();


}
Expand Down

0 comments on commit 764ccc0

Please sign in to comment.