CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
List of all members | Classes | Public Types | Public Member Functions | Protected Types | Protected Member Functions | Static Protected Attributes | Private Types | Private Attributes | Static Private Attributes
SteppingHelixPropagator Class Referencefinal

#include <SteppingHelixPropagator.h>

Inheritance diagram for SteppingHelixPropagator:
Propagator

Classes

struct  Basis
 

Public Types

enum  DestType {
  RADIUS_DT = 0, Z_DT, PLANE_DT, CONE_DT,
  CYLINDER_DT, PATHL_DT, POINT_PCA_DT, LINE_PCA_DT,
  UNDEFINED_DT
}
 
enum  Fancy {
  HEL_AS_F = 0, HEL_ALL_F, POL_1_F, POL_2_F,
  POL_M_F
}
 
enum  Pars { RADIUS_P = 0, Z_P = 0, PATHL_P = 0 }
 
typedef CLHEP::Hep3Vector Point
 
typedef
SteppingHelixStateInfo::Result 		Result
 
typedef SteppingHelixStateInfo StateInfo
 
typedef CLHEP::Hep3Vector Vector
 

Public Member Functions

void applyRadX0Correction (bool applyRadX0Correction)
 
SteppingHelixPropagator * clone () const override
 
const MagneticField * magneticField () const override
 
void propagate (const SteppingHelixStateInfo &ftsStart, const Surface &sDest, SteppingHelixStateInfo &out) const
 Propagate to Plane given a starting point. More...
 
void propagate (const SteppingHelixStateInfo &ftsStart, const Plane &pDest, SteppingHelixStateInfo &out) const
 
void propagate (const SteppingHelixStateInfo &ftsStart, const Cylinder &cDest, SteppingHelixStateInfo &out) const
 Propagate to Cylinder given a starting point (a Cylinder is assumed to be positioned at 0,0,0) More...
 
void propagate (const SteppingHelixStateInfo &ftsStart, const GlobalPoint &pDest, SteppingHelixStateInfo &out) const
 Propagate to PCA to point given a starting point. More...
 
void propagate (const SteppingHelixStateInfo &ftsStart, const GlobalPoint &pDest1, const GlobalPoint &pDest2, SteppingHelixStateInfo &out) const
 Propagate to PCA to a line (given by 2 points) given a starting point. More...
 
std::pair
< TrajectoryStateOnSurface,
double > 		propagateWithPath (const FreeTrajectoryState &ftsStart, const Plane &pDest) const override
 
std::pair
< TrajectoryStateOnSurface,
double > 		propagateWithPath (const FreeTrajectoryState &ftsStart, const Cylinder &cDest) const override
 
std::pair< FreeTrajectoryState,
double > 		propagateWithPath (const FreeTrajectoryState &ftsStart, const GlobalPoint &pDest) const override
 Propagate to PCA to point given a starting point. More...
 
std::pair< FreeTrajectoryState,
double > 		propagateWithPath (const FreeTrajectoryState &ftsStart, const GlobalPoint &pDest1, const GlobalPoint &pDest2) const override
 Propagate to PCA to a line (given by 2 points) given a starting point. More...
 
std::pair< FreeTrajectoryState,
double > 		propagateWithPath (const FreeTrajectoryState &ftsStart, const reco::BeamSpot &beamSpot) const override
 Propagate to PCA to a line (given by beamSpot position and slope) given a starting point. More...
 
void setDebug (bool debug)
 Switch debug printouts (to cout) .. very verbose. More...
 
void setEndcapShiftsInZPosNeg (double valPos, double valNeg)
 set shifts in Z for endcap pieces (includes EE, HE, ME, YE) More...
 
void setMaterialMode (bool noMaterial)
 Switch for material effects mode: no material effects if true. More...
 
void setNoErrorPropagation (bool noErrorPropagation)
 Force no error propagation. More...
 
void setReturnTangentPlane (bool val)
 flag to return tangent plane for non-plane input More...
 
void setSendLogWarning (bool val)
 flag to send LogWarning on failures More...
 
void setUseInTeslaFromMagField (bool val)
 force getting field value from MagneticField, not the geometric one More...
 
void setUseIsYokeFlag (bool val)
 
void setUseMagVolumes (bool val)
 Switch to using MagneticField Volumes .. as in VolumeBasedMagneticField. More...
 
void setUseMatVolumes (bool val)
 Switch to using Material Volumes .. internally defined for now. More...
 
void setUseTuningForL2Speed (bool val)
 
void setVBFPointer (const VolumeBasedMagneticField *val)
 
 SteppingHelixPropagator ()
 Constructors. More...
 
 SteppingHelixPropagator (const MagneticField *field, PropagationDirection dir=alongMomentum)
 
 ~SteppingHelixPropagator () override
 Destructor. More...
 
- Public Member Functions inherited from Propagator
template<typename STA , typename SUR >
TrajectoryStateOnSurface propagate (STA const &state, SUR const &surface) const
 
virtual FreeTrajectoryState propagate (const FreeTrajectoryState &ftsStart, const GlobalPoint &pDest) const final
 
virtual FreeTrajectoryState propagate (const FreeTrajectoryState &ftsStart, const GlobalPoint &pDest1, const GlobalPoint &pDest2) const final
 
virtual FreeTrajectoryState propagate (const FreeTrajectoryState &ftsStart, const reco::BeamSpot &beamSpot) const final
 
virtual std::pair
< TrajectoryStateOnSurface,
double > 		propagateWithPath (const FreeTrajectoryState &, const Surface &) const final
 
virtual std::pair
< TrajectoryStateOnSurface,
double > 		propagateWithPath (const TrajectoryStateOnSurface &tsos, const Surface &sur) const final
 
virtual std::pair
< TrajectoryStateOnSurface,
double > 		propagateWithPath (const TrajectoryStateOnSurface &tsos, const Plane &sur) const
 
virtual std::pair
< TrajectoryStateOnSurface,
double > 		propagateWithPath (const TrajectoryStateOnSurface &tsos, const Cylinder &sur) const
 
virtual PropagationDirection propagationDirection () const final
 
 Propagator (PropagationDirection dir=alongMomentum)
 
virtual bool setMaxDirectionChange (float phiMax)
 
virtual void setPropagationDirection (PropagationDirection dir)
 
virtual ~Propagator ()
 

Protected Types

typedef
SteppingHelixStateInfo::VolumeBounds 		MatBounds
 
typedef StateInfo StateArray [MAX_POINTS+1]
 

Protected Member Functions

int cIndex_ (int ind) const
 (Internals) circular index for array of transients More...
 
double getDeDx (const SteppingHelixPropagator::StateInfo &sv, double &dEdXPrime, double &radX0, MatBounds &rzLims) const
 estimate average (in fact smth. close to MPV and median) energy loss per unit path length More...
 
void getNextState (const SteppingHelixPropagator::StateInfo &svPrevious, SteppingHelixPropagator::StateInfo &svNext, double dP, const SteppingHelixPropagator::Vector &tau, const SteppingHelixPropagator::Vector &drVec, double dS, double dX0, const AlgebraicMatrix55 &dCovCurv) const
 
void initStateArraySHPSpecific (StateArray &svBuf, bool flagsOnly) const
 (Internals) set defaults needed for states used inside propagate methods More...
 
bool isYokeVolume (const MagVolume *vol) const
 check if it's a yoke/iron based on this MagVol internals More...
 
void loadState (SteppingHelixPropagator::StateInfo &svCurrent, const SteppingHelixPropagator::Vector &p3, const SteppingHelixPropagator::Point &r3, int charge, PropagationDirection dir, const AlgebraicSymMatrix55 &covCurv) const
 
bool makeAtomStep (SteppingHelixPropagator::StateInfo &svCurrent, SteppingHelixPropagator::StateInfo &svNext, double dS, PropagationDirection dir, SteppingHelixPropagator::Fancy fancy) const
 main stepping function: compute next state vector after a step of length dS More...
 
Result propagate (StateArray &svBuff, int &nPoints, SteppingHelixPropagator::DestType type, const double pars[6], double epsilon=1e-3) const
 
Result refToDest (DestType dest, const SteppingHelixPropagator::StateInfo &sv, const double pars[6], double &dist, double &tanDist, PropagationDirection &refDirection, double fastSkipDist=1e12) const
 (Internals) determine distance and direction from the current position to the plane More...
 
Result refToMagVolume (const SteppingHelixPropagator::StateInfo &sv, PropagationDirection dir, double &dist, double &tanDist, double fastSkipDist=1e12, bool expectNewMagVolume=false, double maxStep=1e12) const
 
Result refToMatVolume (const SteppingHelixPropagator::StateInfo &sv, PropagationDirection dir, double &dist, double &tanDist, double fastSkipDist=1e12) const
 
void setIState (const SteppingHelixStateInfo &sStart, StateArray &svBuff, int &nPoints) const
 (Internals) Init starting point More...
 
void setRep (SteppingHelixPropagator::Basis &rep, const SteppingHelixPropagator::Vector &tau) const
 Set/compute basis vectors for local coordinates at current step (called by incrementState) More...
 

Static Protected Attributes

static const int MAX_POINTS = 7
 

Private Types

typedef std::pair
< FreeTrajectoryState, double > 		FtsPP
 
typedef std::pair
< TrajectoryStateOnSurface,
double > 		TsosPP
 

Private Attributes

bool applyRadX0Correction_
 
bool debug_
 
double defaultStep_
 
double ecShiftNeg_
 
double ecShiftPos_
 
const MagneticField * field_
 
StateInfo invalidState_
 
bool noErrorPropagation_
 
bool noMaterialMode_
 
bool returnTangentPlane_
 
bool sendLogWarning_
 
const AlgebraicSymMatrix55 unit55_
 
bool useInTeslaFromMagField_
 
bool useIsYokeFlag_
 
bool useMagVolumes_
 
bool useMatVolumes_
 
bool useTuningForL2Speed_
 
const VolumeBasedMagneticField * vbField_
 

Static Private Attributes

static const int MAX_STEPS = 10000
 

Detailed Description

Propagator implementation using steps along a helix. Minimal geometry navigation. Material effects (multiple scattering and energy loss) are based on tuning to MC and (eventually) data.

Author
Vyacheslav Krutelyov (slava77)

Propagator implementation using steps along a helix. Minimal geometry navigation. Material effects (multiple scattering and energy loss) are based on tuning to MC and (eventually) data. Implementation file contents follow.

Author
Vyacheslav Krutelyov (slava77)

Definition at line 36 of file SteppingHelixPropagator.h.

Member Typedef Documentation

typedef std::pair<FreeTrajectoryState, double> SteppingHelixPropagator::FtsPP
private

Definition at line 258 of file SteppingHelixPropagator.h.

typedef SteppingHelixStateInfo::VolumeBounds SteppingHelixPropagator::MatBounds
protected

Definition at line 170 of file SteppingHelixPropagator.h.

typedef CLHEP::Hep3Vector SteppingHelixPropagator::Point

Definition at line 39 of file SteppingHelixPropagator.h.

typedef SteppingHelixStateInfo::Result SteppingHelixPropagator::Result

Definition at line 42 of file SteppingHelixPropagator.h.

typedef StateInfo SteppingHelixPropagator::StateArray[MAX_POINTS+1]
protected

Definition at line 172 of file SteppingHelixPropagator.h.

typedef SteppingHelixStateInfo SteppingHelixPropagator::StateInfo

Definition at line 41 of file SteppingHelixPropagator.h.

typedef std::pair<TrajectoryStateOnSurface, double> SteppingHelixPropagator::TsosPP
private

Definition at line 257 of file SteppingHelixPropagator.h.

typedef CLHEP::Hep3Vector SteppingHelixPropagator::Vector

Definition at line 38 of file SteppingHelixPropagator.h.

Member Enumeration Documentation

enum SteppingHelixPropagator::DestType
Enumerator
RADIUS_DT 
Z_DT 
PLANE_DT 
CONE_DT 
CYLINDER_DT 
PATHL_DT 
POINT_PCA_DT 
LINE_PCA_DT 
UNDEFINED_DT 

Definition at line 52 of file SteppingHelixPropagator.h.

52  {
53  RADIUS_DT = 0,
54  Z_DT,
55  PLANE_DT,
56  CONE_DT,
57  CYLINDER_DT,
58  PATHL_DT,
59  POINT_PCA_DT,
60  LINE_PCA_DT,
61  UNDEFINED_DT
62  };
enum SteppingHelixPropagator::Fancy
Enumerator
HEL_AS_F 
HEL_ALL_F 
POL_1_F 
POL_2_F 
POL_M_F 

Definition at line 64 of file SteppingHelixPropagator.h.

64  {
65  HEL_AS_F = 0, //simple analytical helix, eloss at end of step
66  HEL_ALL_F, //analytical helix with linear eloss
67  POL_1_F, //1st order approximation, straight line
68  POL_2_F, //2nd order
69  POL_M_F //highest available
70  };
enum SteppingHelixPropagator::Pars
Enumerator
RADIUS_P 
Z_P 
PATHL_P 

Definition at line 50 of file SteppingHelixPropagator.h.

50 { RADIUS_P = 0, Z_P = 0, PATHL_P = 0 };

Constructor & Destructor Documentation

SteppingHelixPropagator::SteppingHelixPropagator ( )

Constructors.

Definition at line 56 of file SteppingHelixPropagator.cc.

References field_.

Referenced by clone().

56 : Propagator(anyDirection) { field_ = nullptr; }
Propagator::Propagator
Propagator(PropagationDirection dir=alongMomentum)
Definition: Propagator.h:46
SteppingHelixPropagator::field_
const MagneticField * field_
Definition: SteppingHelixPropagator.h:265
SteppingHelixPropagator::SteppingHelixPropagator ( const MagneticField *  field,
PropagationDirection  dir = alongMomentum 
)

Definition at line 58 of file SteppingHelixPropagator.cc.

References applyRadX0Correction_, debug_, defaultStep_, ecShiftNeg_, ecShiftPos_, field_, noErrorPropagation_, noMaterialMode_, returnTangentPlane_, sendLogWarning_, useInTeslaFromMagField_, useIsYokeFlag_, useMagVolumes_, useMatVolumes_, useTuningForL2Speed_, and vbField_.

59  : Propagator(dir), unit55_(AlgebraicMatrixID()) {
60  field_ = field;
61  vbField_ = dynamic_cast<const VolumeBasedMagneticField*>(field_);
62  debug_ = false;
63  noMaterialMode_ = false;
64  noErrorPropagation_ = false;
65  applyRadX0Correction_ = true;
66  useMagVolumes_ = true;
67  useIsYokeFlag_ = true;
68  useMatVolumes_ = true;
69  useInTeslaFromMagField_ = false; //overrides behavior only if true
70  returnTangentPlane_ = true;
71  sendLogWarning_ = false;
72  useTuningForL2Speed_ = false;
73  defaultStep_ = 5.;
74 
75  ecShiftPos_ = 0;
76  ecShiftNeg_ = 0;
77 }
Propagator::Propagator
Propagator(PropagationDirection dir=alongMomentum)
Definition: Propagator.h:46
AlgebraicMatrixID
ROOT::Math::SMatrixIdentity AlgebraicMatrixID
Definition: AlgebraicROOTObjects.h:72
DeadROC_duringRun.dir
string dir
Definition: DeadROC_duringRun.py:23
SteppingHelixPropagator::vbField_
const VolumeBasedMagneticField * vbField_
Definition: SteppingHelixPropagator.h:266
SteppingHelixPropagator::unit55_
const AlgebraicSymMatrix55 unit55_
Definition: SteppingHelixPropagator.h:267
SteppingHelixPropagator::field_
const MagneticField * field_
Definition: SteppingHelixPropagator.h:265
SteppingHelixPropagator::~SteppingHelixPropagator ( )
override

Destructor.

Definition at line 54 of file SteppingHelixPropagator.cc.

54 {}

Member Function Documentation

void SteppingHelixPropagator::applyRadX0Correction ( bool  applyRadX0Correction)
inline

Apply radLength correction (1+0.036*ln(radX0+1)) to covariance matrix +1 is added for IR-safety Should be done with care: it's easy to make the end-point result dependent on the intermediate stop points

Definition at line 133 of file SteppingHelixPropagator.h.

References applyRadX0Correction(), and applyRadX0Correction_.

Referenced by applyRadX0Correction(), and AlCaHOCalibProducer::fillHOStore().

133 { applyRadX0Correction_ = applyRadX0Correction; }
SteppingHelixPropagator::applyRadX0Correction
void applyRadX0Correction(bool applyRadX0Correction)
Definition: SteppingHelixPropagator.h:133
int SteppingHelixPropagator::cIndex_ ( int  ind) const
protected

(Internals) circular index for array of transients

Definition at line 1773 of file SteppingHelixPropagator.cc.

References MAX_POINTS, and mps_fire::result.

Referenced by propagate(), and setIState().

1773  {
1774  int result = ind % MAX_POINTS;
1775  if (ind != 0 && result == 0) {
1776  result = MAX_POINTS;
1777  }
1778  return result;
1779 }
mps_fire.result
tuple result
Definition: mps_fire.py:311
SteppingHelixPropagator* SteppingHelixPropagator::clone ( void  ) const
inlineoverridevirtual

Implements Propagator.

Definition at line 76 of file SteppingHelixPropagator.h.

References SteppingHelixPropagator().

76 { return new SteppingHelixPropagator(*this); }
double SteppingHelixPropagator::getDeDx ( const SteppingHelixPropagator::StateInfo &  sv,
double &  dEdXPrime,
double &  radX0,
MatBounds &  rzLims 
) const
protected

estimate average (in fact smth. close to MPV and median) energy loss per unit path length

Definition at line 1317 of file SteppingHelixPropagator.cc.

References SteppingHelixStateInfo::bf, plot_hgcal_utils::dEdx, alignCSCRings::e, ecShiftNeg_, ecShiftPos_, SteppingHelixStateInfo::isYokeVol, log, SteppingHelixStateInfo::magVol, noMaterialMode_, SteppingHelixStateInfo::p3, SteppingHelixStateInfo::r3, mathSSE::sqrt(), useIsYokeFlag_, and useMagVolumes_.

Referenced by getNextState(), loadState(), and makeAtomStep().

1320  {
1321  radX0 = 1.e24;
1322  dEdXPrime = 0.;
1323  rzLims = MatBounds();
1324  if (noMaterialMode_)
1325  return 0;
1326 
1327  double dEdx = 0.;
1328 
1329  double lR = sv.r3.perp();
1330  double lZ = fabs(sv.r3.z());
1331 
1332  //assume "Iron" .. seems to be quite the same for brass/iron/PbW04
1333  //good for Fe within 3% for 0.2 GeV to 10PeV
1334 
1335  double dEdX_HCal = 0.95; //extracted from sim
1336  double dEdX_ECal = 0.45;
1337  double dEdX_coil = 0.35; //extracted from sim .. closer to 40% in fact
1338  double dEdX_Fe = 1;
1339  double dEdX_MCh = 0.053; //chambers on average
1340  double dEdX_Trk = 0.0114;
1341  double dEdX_Air = 2E-4;
1342  double dEdX_Vac = 0.0;
1343 
1344  double radX0_HCal = 1.44 / 0.8; //guessing
1345  double radX0_ECal = 0.89 / 0.7;
1346  double radX0_coil = 4.; //
1347  double radX0_Fe = 1.76;
1348  double radX0_MCh = 1e3; //
1349  double radX0_Trk = 320.;
1350  double radX0_Air = 3.e4;
1351  double radX0_Vac = 3.e9; //"big" number for vacuum
1352 
1353  //not all the boundaries are set below: this will be a default
1354  if (!(lR < 380 && lZ < 785)) {
1355  if (lZ > 785)
1356  rzLims = MatBounds(0, 1e4, 785, 1e4);
1357  if (lZ < 785)
1358  rzLims = MatBounds(380, 1e4, 0, 785);
1359  }
1360 
1361  //this def makes sense assuming we don't jump into endcap volume from the other z-side in one step
1362  //also, it is a positive shift only (at least for now): can't move ec further into the detector
1363  double ecShift = sv.r3.z() > 0 ? fabs(ecShiftPos_) : fabs(ecShiftNeg_);
1364 
1365  //this should roughly figure out where things are
1366  //(numbers taken from Fig1.1.2 TDR and from geom xmls)
1367  if (lR < 2.9) { //inside beampipe
1368  dEdx = dEdX_Vac;
1369  radX0 = radX0_Vac;
1370  rzLims = MatBounds(0, 2.9, 0, 1E4);
1371  } else if (lR < 129) {
1372  if (lZ < 294) {
1373  rzLims = MatBounds(2.9, 129, 0, 294); //tracker boundaries
1374  dEdx = dEdX_Trk;
1375  radX0 = radX0_Trk;
1376  //somewhat empirical formula that ~ matches the average if going from 0,0,0
1377  //assuming "uniform" tracker material
1378  //doesn't really track material layer to layer
1379  double lEtaDet = fabs(sv.r3.eta());
1380  double scaleRadX = lEtaDet > 1.5 ? 0.7724 : 1. / cosh(0.5 * lEtaDet); //sin(2.*atan(exp(-0.5*lEtaDet)));
1381  scaleRadX *= scaleRadX;
1382  if (lEtaDet > 2 && lZ > 20)
1383  scaleRadX *= (lEtaDet - 1.);
1384  if (lEtaDet > 2.5 && lZ > 20)
1385  scaleRadX *= (lEtaDet - 1.);
1386  radX0 *= scaleRadX;
1387  }
1388  //endcap part begins here
1389  else if (lZ < 294 + ecShift) {
1390  //gap in front of EE here, piece inside 2.9<R<129
1391  rzLims = MatBounds(2.9, 129, 294, 294 + ecShift);
1392  dEdx = dEdX_Air;
1393  radX0 = radX0_Air;
1394  } else if (lZ < 372 + ecShift) {
1395  rzLims = MatBounds(2.9, 129, 294 + ecShift, 372 + ecShift); //EE here, piece inside 2.9<R<129
1396  dEdx = dEdX_ECal;
1397  radX0 = radX0_ECal;
1398  } //EE averaged out over a larger space
1399  else if (lZ < 398 + ecShift) {
1400  rzLims =
1401  MatBounds(2.9, 129, 372 + ecShift, 398 + ecShift); //whatever goes behind EE 2.9<R<129 is air up to Z=398
1402  dEdx = dEdX_HCal * 0.05;
1403  radX0 = radX0_Air;
1404  } //betw EE and HE
1405  else if (lZ < 555 + ecShift) {
1406  rzLims = MatBounds(2.9, 129, 398 + ecShift, 555 + ecShift); //HE piece 2.9<R<129;
1407  dEdx = dEdX_HCal * 0.96;
1408  radX0 = radX0_HCal / 0.96;
1409  } //HE calor abit less dense
1410  else {
1411  // rzLims = MatBounds(2.9, 129, 555, 785);
1412  // set the boundaries first: they serve as stop-points here
1413  // the material is set below
1414  if (lZ < 568 + ecShift)
1415  rzLims = MatBounds(2.9, 129, 555 + ecShift, 568 + ecShift); //a piece of HE support R<129, 555<Z<568
1416  else if (lZ < 625 + ecShift) {
1417  if (lR < 85 + ecShift)
1418  rzLims = MatBounds(2.9, 85, 568 + ecShift, 625 + ecShift);
1419  else
1420  rzLims = MatBounds(85, 129, 568 + ecShift, 625 + ecShift);
1421  } else if (lZ < 785 + ecShift)
1422  rzLims = MatBounds(2.9, 129, 625 + ecShift, 785 + ecShift);
1423  else if (lZ < 1500 + ecShift)
1424  rzLims = MatBounds(2.9, 129, 785 + ecShift, 1500 + ecShift);
1425  else
1426  rzLims = MatBounds(2.9, 129, 1500 + ecShift, 1E4);
1427 
1428  //iron .. don't care about no material in front of HF (too forward)
1429  if (!(lZ > 568 + ecShift && lZ < 625 + ecShift && lR > 85) // HE support
1430  && !(lZ > 785 + ecShift && lZ < 850 + ecShift && lR > 118)) {
1431  dEdx = dEdX_Fe;
1432  radX0 = radX0_Fe;
1433  } else {
1434  dEdx = dEdX_MCh;
1435  radX0 = radX0_MCh;
1436  } //ME at eta > 2.2
1437  }
1438  } else if (lR < 287) {
1439  if (lZ < 372 + ecShift && lR < 177) { // 129<<R<177
1440  if (lZ < 304)
1441  rzLims = MatBounds(129, 177, 0, 304); //EB 129<<R<177 0<Z<304
1442  else if (lZ < 343) { // 129<<R<177 304<Z<343
1443  if (lR < 135)
1444  rzLims = MatBounds(129, 135, 304, 343); // tk piece 129<<R<135 304<Z<343
1445  else if (lR < 172) { //
1446  if (lZ < 311)
1447  rzLims = MatBounds(135, 172, 304, 311);
1448  else
1449  rzLims = MatBounds(135, 172, 311, 343);
1450  } else {
1451  if (lZ < 328)
1452  rzLims = MatBounds(172, 177, 304, 328);
1453  else
1454  rzLims = MatBounds(172, 177, 328, 343);
1455  }
1456  } else if (lZ < 343 + ecShift) {
1457  rzLims = MatBounds(129, 177, 343, 343 + ecShift); //gap
1458  } else {
1459  if (lR < 156)
1460  rzLims = MatBounds(129, 156, 343 + ecShift, 372 + ecShift);
1461  else if ((lZ - 343 - ecShift) > (lR - 156) * 1.38)
1462  rzLims = MatBounds(156, 177, 127.72 + ecShift, 372 + ecShift, atan(1.38), 0);
1463  else
1464  rzLims = MatBounds(156, 177, 343 + ecShift, 127.72 + ecShift, 0, atan(1.38));
1465  }
1466 
1467  if (!(lR > 135 && lZ < 343 + ecShift && lZ > 304) &&
1468  !(lR > 156 && lZ < 372 + ecShift && lZ > 343 + ecShift && ((lZ - 343. - ecShift) < (lR - 156.) * 1.38))) {
1469  //the crystals are the same length, but they are not 100% of material
1470  double cosThetaEquiv = 0.8 / sqrt(1. + lZ * lZ / lR / lR) + 0.2;
1471  if (lZ > 343)
1472  cosThetaEquiv = 1.;
1473  dEdx = dEdX_ECal * cosThetaEquiv;
1474  radX0 = radX0_ECal / cosThetaEquiv;
1475  } //EB
1476  else {
1477  if ((lZ > 304 && lZ < 328 && lR < 177 && lR > 135) && !(lZ > 311 && lR < 172)) {
1478  dEdx = dEdX_Fe;
1479  radX0 = radX0_Fe;
1480  } //Tk_Support
1481  else if (lZ > 343 && lZ < 343 + ecShift) {
1482  dEdx = dEdX_Air;
1483  radX0 = radX0_Air;
1484  } else {
1485  dEdx = dEdX_ECal * 0.2;
1486  radX0 = radX0_Air;
1487  } //cables go here <-- will be abit too dense for ecShift > 0
1488  }
1489  } else if (lZ < 554 + ecShift) { // 129<R<177 372<Z<554 AND 177<R<287 0<Z<554
1490  if (lR < 177) { // 129<R<177 372<Z<554
1491  if (lZ > 372 + ecShift && lZ < 398 + ecShift)
1492  rzLims = MatBounds(129, 177, 372 + ecShift, 398 + ecShift); // EE 129<R<177 372<Z<398
1493  else if (lZ < 548 + ecShift)
1494  rzLims = MatBounds(129, 177, 398 + ecShift, 548 + ecShift); // HE 129<R<177 398<Z<548
1495  else
1496  rzLims = MatBounds(129, 177, 548 + ecShift, 554 + ecShift); // HE gap 129<R<177 548<Z<554
1497  } else if (lR < 193) { // 177<R<193 0<Z<554
1498  if ((lZ - 307) < (lR - 177.) * 1.739)
1499  rzLims = MatBounds(177, 193, 0, -0.803, 0, atan(1.739));
1500  else if (lZ < 389)
1501  rzLims = MatBounds(177, 193, -0.803, 389, atan(1.739), 0.);
1502  else if (lZ < 389 + ecShift)
1503  rzLims = MatBounds(177, 193, 389, 389 + ecShift); // air insert
1504  else if (lZ < 548 + ecShift)
1505  rzLims = MatBounds(177, 193, 389 + ecShift, 548 + ecShift); // HE 177<R<193 389<Z<548
1506  else
1507  rzLims = MatBounds(177, 193, 548 + ecShift, 554 + ecShift); // HE gap 177<R<193 548<Z<554
1508  } else if (lR < 264) { // 193<R<264 0<Z<554
1509  double anApex = 375.7278 - 193. / 1.327; // 230.28695599096
1510  if ((lZ - 375.7278) < (lR - 193.) / 1.327)
1511  rzLims = MatBounds(193, 264, 0, anApex, 0, atan(1. / 1.327)); //HB
1512  else if ((lZ - 392.7278) < (lR - 193.) / 1.327)
1513  rzLims = MatBounds(193, 264, anApex, anApex + 17., atan(1. / 1.327), atan(1. / 1.327)); // HB-HE gap
1514  else if ((lZ - 392.7278 - ecShift) < (lR - 193.) / 1.327)
1515  rzLims = MatBounds(193,
1516  264,
1517  anApex + 17.,
1518  anApex + 17. + ecShift,
1519  atan(1. / 1.327),
1520  atan(1. / 1.327)); // HB-HE gap air insert
1521  // HE (372,193)-(517,193)-(517,264)-(417.8,264)
1522  else if (lZ < 517 + ecShift)
1523  rzLims = MatBounds(193, 264, anApex + 17. + ecShift, 517 + ecShift, atan(1. / 1.327), 0);
1524  else if (lZ < 548 + ecShift) { // HE+gap 193<R<264 517<Z<548
1525  if (lR < 246)
1526  rzLims = MatBounds(193, 246, 517 + ecShift, 548 + ecShift); // HE 193<R<246 517<Z<548
1527  else
1528  rzLims = MatBounds(246, 264, 517 + ecShift, 548 + ecShift); // HE gap 246<R<264 517<Z<548
1529  } else
1530  rzLims = MatBounds(193, 264, 548 + ecShift, 554 + ecShift); // HE gap 193<R<246 548<Z<554
1531  } else if (lR < 275) { // 264<R<275 0<Z<554
1532  if (lZ < 433)
1533  rzLims = MatBounds(264, 275, 0, 433); //HB
1534  else if (lZ < 554)
1535  rzLims = MatBounds(264, 275, 433, 554); // HE gap 264<R<275 433<Z<554
1536  else
1537  rzLims = MatBounds(264, 275, 554, 554 + ecShift); // HE gap 264<R<275 554<Z<554 air insert
1538  } else { // 275<R<287 0<Z<554
1539  if (lZ < 402)
1540  rzLims = MatBounds(275, 287, 0, 402); // HB 275<R<287 0<Z<402
1541  else if (lZ < 554)
1542  rzLims = MatBounds(275, 287, 402, 554); // //HE gap 275<R<287 402<Z<554
1543  else
1544  rzLims = MatBounds(275, 287, 554, 554 + ecShift); //HE gap 275<R<287 554<Z<554 air insert
1545  }
1546 
1547  if ((lZ < 433 || lR < 264) && (lZ < 402 || lR < 275) &&
1548  (lZ < 517 + ecShift || lR < 246) //notches
1549  //I should've made HE and HF different .. now need to shorten HE to match
1550  && lZ < 548 + ecShift && !(lZ < 389 + ecShift && lZ > 335 && lR < 193) //not a gap EB-EE 129<R<193
1551  && !(lZ > 307 && lZ < 335 && lR < 193 && ((lZ - 307) > (lR - 177.) * 1.739)) //not a gap
1552  && !(lR < 177 && lZ < 398 + ecShift) //under the HE nose
1553  && !(lR < 264 && lR > 193 &&
1554  fabs(441.5 + 0.5 * ecShift - lZ + (lR - 269.) / 1.327) < (8.5 + ecShift * 0.5))) //not a gap
1555  {
1556  dEdx = dEdX_HCal;
1557  radX0 = radX0_HCal; //hcal
1558  } else if ((lR < 193 && lZ > 389 && lZ < 389 + ecShift) ||
1559  (lR > 264 && lR < 287 && lZ > 554 && lZ < 554 + ecShift) ||
1560  (lR < 264 && lR > 193 &&
1561  fabs(441.5 + 8.5 + 0.5 * ecShift - lZ + (lR - 269.) / 1.327) < ecShift * 0.5)) {
1562  dEdx = dEdX_Air;
1563  radX0 = radX0_Air;
1564  } else {
1565  dEdx = dEdX_HCal * 0.05;
1566  radX0 = radX0_Air;
1567  } //endcap gap
1568  }
1569  //the rest is a tube of endcap volume 129<R<251 554<Z<largeValue
1570  else if (lZ < 564 + ecShift) { // 129<R<287 554<Z<564
1571  if (lR < 251) {
1572  rzLims = MatBounds(129, 251, 554 + ecShift, 564 + ecShift); // HE support 129<R<251 554<Z<564
1573  dEdx = dEdX_Fe;
1574  radX0 = radX0_Fe;
1575  } //HE support
1576  else {
1577  rzLims = MatBounds(251, 287, 554 + ecShift, 564 + ecShift); //HE/ME gap 251<R<287 554<Z<564
1578  dEdx = dEdX_MCh;
1579  radX0 = radX0_MCh;
1580  }
1581  } else if (lZ < 625 + ecShift) { // ME/1/1 129<R<287 564<Z<625
1582  rzLims = MatBounds(129, 287, 564 + ecShift, 625 + ecShift);
1583  dEdx = dEdX_MCh;
1584  radX0 = radX0_MCh;
1585  } else if (lZ < 785 + ecShift) { //129<R<287 625<Z<785
1586  if (lR < 275)
1587  rzLims = MatBounds(129, 275, 625 + ecShift, 785 + ecShift); //YE/1 129<R<275 625<Z<785
1588  else { // 275<R<287 625<Z<785
1589  if (lZ < 720 + ecShift)
1590  rzLims = MatBounds(275, 287, 625 + ecShift, 720 + ecShift); // ME/1/2 275<R<287 625<Z<720
1591  else
1592  rzLims = MatBounds(275, 287, 720 + ecShift, 785 + ecShift); // YE/1 275<R<287 720<Z<785
1593  }
1594  if (!(lR > 275 && lZ < 720 + ecShift)) {
1595  dEdx = dEdX_Fe;
1596  radX0 = radX0_Fe;
1597  } //iron
1598  else {
1599  dEdx = dEdX_MCh;
1600  radX0 = radX0_MCh;
1601  }
1602  } else if (lZ < 850 + ecShift) {
1603  rzLims = MatBounds(129, 287, 785 + ecShift, 850 + ecShift);
1604  dEdx = dEdX_MCh;
1605  radX0 = radX0_MCh;
1606  } else if (lZ < 910 + ecShift) {
1607  rzLims = MatBounds(129, 287, 850 + ecShift, 910 + ecShift);
1608  dEdx = dEdX_Fe;
1609  radX0 = radX0_Fe;
1610  } //iron
1611  else if (lZ < 975 + ecShift) {
1612  rzLims = MatBounds(129, 287, 910 + ecShift, 975 + ecShift);
1613  dEdx = dEdX_MCh;
1614  radX0 = radX0_MCh;
1615  } else if (lZ < 1000 + ecShift) {
1616  rzLims = MatBounds(129, 287, 975 + ecShift, 1000 + ecShift);
1617  dEdx = dEdX_Fe;
1618  radX0 = radX0_Fe;
1619  } //iron
1620  else if (lZ < 1063 + ecShift) {
1621  rzLims = MatBounds(129, 287, 1000 + ecShift, 1063 + ecShift);
1622  dEdx = dEdX_MCh;
1623  radX0 = radX0_MCh;
1624  } else if (lZ < 1073 + ecShift) {
1625  rzLims = MatBounds(129, 287, 1063 + ecShift, 1073 + ecShift);
1626  dEdx = dEdX_Fe;
1627  radX0 = radX0_Fe;
1628  } else if (lZ < 1E4) {
1629  rzLims = MatBounds(129, 287, 1073 + ecShift, 1E4);
1630  dEdx = dEdX_Air;
1631  radX0 = radX0_Air;
1632  } else {
1633  dEdx = dEdX_Air;
1634  radX0 = radX0_Air;
1635  }
1636  } else if (lR < 380 && lZ < 667) {
1637  if (lZ < 630) {
1638  if (lR < 315)
1639  rzLims = MatBounds(287, 315, 0, 630);
1640  else if (lR < 341)
1641  rzLims = MatBounds(315, 341, 0, 630); //b-field ~linear rapid fall here
1642  else
1643  rzLims = MatBounds(341, 380, 0, 630);
1644  } else
1645  rzLims = MatBounds(287, 380, 630, 667);
1646 
1647  if (lZ < 630) {
1648  dEdx = dEdX_coil;
1649  radX0 = radX0_coil;
1650  } //a guess for the solenoid average
1651  else {
1652  dEdx = dEdX_Air;
1653  radX0 = radX0_Air;
1654  } //endcap gap
1655  } else {
1656  if (lZ < 667) {
1657  if (lR < 850) {
1658  bool isIron = false;
1659  //sanity check in addition to flags
1660  if (useIsYokeFlag_ && useMagVolumes_ && sv.magVol != nullptr) {
1661  isIron = sv.isYokeVol;
1662  } else {
1663  double bMag = sv.bf.mag();
1664  isIron = (bMag > 0.75 && !(lZ > 500 && lR < 500 && bMag < 1.15) && !(lZ < 450 && lR > 420 && bMag < 1.15));
1665  }
1666  //tell the material stepper where mat bounds are
1667  rzLims = MatBounds(380, 850, 0, 667);
1668  if (isIron) {
1669  dEdx = dEdX_Fe;
1670  radX0 = radX0_Fe;
1671  } //iron
1672  else {
1673  dEdx = dEdX_MCh;
1674  radX0 = radX0_MCh;
1675  }
1676  } else {
1677  rzLims = MatBounds(850, 1E4, 0, 667);
1678  dEdx = dEdX_Air;
1679  radX0 = radX0_Air;
1680  }
1681  } else if (lR > 750) {
1682  rzLims = MatBounds(750, 1E4, 667, 1E4);
1683  dEdx = dEdX_Air;
1684  radX0 = radX0_Air;
1685  } else if (lZ < 667 + ecShift) {
1686  rzLims = MatBounds(287, 750, 667, 667 + ecShift);
1687  dEdx = dEdX_Air;
1688  radX0 = radX0_Air;
1689  }
1690  //the rest is endcap piece with 287<R<750 Z>667
1691  else if (lZ < 724 + ecShift) {
1692  if (lR < 380)
1693  rzLims = MatBounds(287, 380, 667 + ecShift, 724 + ecShift);
1694  else
1695  rzLims = MatBounds(380, 750, 667 + ecShift, 724 + ecShift);
1696  dEdx = dEdX_MCh;
1697  radX0 = radX0_MCh;
1698  } else if (lZ < 785 + ecShift) {
1699  if (lR < 380)
1700  rzLims = MatBounds(287, 380, 724 + ecShift, 785 + ecShift);
1701  else
1702  rzLims = MatBounds(380, 750, 724 + ecShift, 785 + ecShift);
1703  dEdx = dEdX_Fe;
1704  radX0 = radX0_Fe;
1705  } //iron
1706  else if (lZ < 850 + ecShift) {
1707  rzLims = MatBounds(287, 750, 785 + ecShift, 850 + ecShift);
1708  dEdx = dEdX_MCh;
1709  radX0 = radX0_MCh;
1710  } else if (lZ < 910 + ecShift) {
1711  rzLims = MatBounds(287, 750, 850 + ecShift, 910 + ecShift);
1712  dEdx = dEdX_Fe;
1713  radX0 = radX0_Fe;
1714  } //iron
1715  else if (lZ < 975 + ecShift) {
1716  rzLims = MatBounds(287, 750, 910 + ecShift, 975 + ecShift);
1717  dEdx = dEdX_MCh;
1718  radX0 = radX0_MCh;
1719  } else if (lZ < 1000 + ecShift) {
1720  rzLims = MatBounds(287, 750, 975 + ecShift, 1000 + ecShift);
1721  dEdx = dEdX_Fe;
1722  radX0 = radX0_Fe;
1723  } //iron
1724  else if (lZ < 1063 + ecShift) {
1725  if (lR < 360) {
1726  rzLims = MatBounds(287, 360, 1000 + ecShift, 1063 + ecShift);
1727  dEdx = dEdX_MCh;
1728  radX0 = radX0_MCh;
1729  }
1730  //put empty air where me4/2 should be (future)
1731  else {
1732  rzLims = MatBounds(360, 750, 1000 + ecShift, 1063 + ecShift);
1733  dEdx = dEdX_Air;
1734  radX0 = radX0_Air;
1735  }
1736  } else if (lZ < 1073 + ecShift) {
1737  rzLims = MatBounds(287, 750, 1063 + ecShift, 1073 + ecShift);
1738  //this plate does not exist: air
1739  dEdx = dEdX_Air;
1740  radX0 = radX0_Air;
1741  } else if (lZ < 1E4) {
1742  rzLims = MatBounds(287, 750, 1073 + ecShift, 1E4);
1743  dEdx = dEdX_Air;
1744  radX0 = radX0_Air;
1745  } else {
1746  dEdx = dEdX_Air;
1747  radX0 = radX0_Air;
1748  } //air
1749  }
1750 
1751  //dEdx so far is a relative number (relative to iron)
1752  //scale by what's expected for iron (the function was fit from pdg table)
1753  //0.065 (PDG) --> 0.044 to better match with MPV
1754  double p0 = sv.p3.mag();
1755  double logp0 = log(p0);
1756  double p0powN33 = 0;
1757  if (p0 > 3.) {
1758  // p0powN33 = exp(-0.33*logp0); //calculate for p>3GeV
1759  double xx = 1. / p0;
1760  xx = sqrt(sqrt(sqrt(sqrt(xx))));
1761  xx = xx * xx * xx * xx * xx; // this is (p0)**(-5/16), close enough to -0.33
1762  p0powN33 = xx;
1763  }
1764  double dEdX_mat = -(11.4 + 0.96 * fabs(logp0 + log(2.8)) + 0.033 * p0 * (1.0 - p0powN33)) * 1e-3;
1765  //in GeV/cm .. 0.8 to get closer to the median or MPV
1766 
1767  dEdXPrime = dEdx == 0 ? 0 : -dEdx * (0.96 / p0 + 0.033 - 0.022 * p0powN33) * 1e-3; //== d(dEdX)/dp
1768  dEdx *= dEdX_mat;
1769 
1770  return dEdx;
1771 }
log
static std::vector< std::string > checklist log
Definition: LHERunInfoProduct.cc:179
SteppingHelixPropagator::MatBounds
SteppingHelixStateInfo::VolumeBounds MatBounds
Definition: SteppingHelixPropagator.h:170
plot_hgcal_utils.dEdx
tuple dEdx
Definition: plot_hgcal_utils.py:216
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
void SteppingHelixPropagator::getNextState ( const SteppingHelixPropagator::StateInfo &  svPrevious,
SteppingHelixPropagator::StateInfo &  svNext,
double  dP,
const SteppingHelixPropagator::Vector &  tau,
const SteppingHelixPropagator::Vector &  drVec,
double  dS,
double  dX0,
const AlgebraicMatrix55 &  dCovCurv 
) const
protected

(Internals) compute transients for current point (increments step counter). Called by makeAtomStep

Definition at line 775 of file SteppingHelixPropagator.cc.

References SteppingHelixStateInfo::bf, SteppingHelixStateInfo::covCurv, debug_, SteppingHelixStateInfo::dEdx, plot_hgcal_utils::dEdx, SteppingHelixStateInfo::dEdXPrime, SteppingHelixStateInfo::dir, alignCSCRings::e, field_, VolumeBasedMagneticField::findVolume(), getDeDx(), SteppingHelixStateInfo::hasErrorPropagated_, MagneticField::inTesla(), MagVolume::inTesla(), SteppingHelixStateInfo::isYokeVol, isYokeVolume(), LogTrace, SteppingHelixStateInfo::magVol, SteppingHelixStateInfo::matDCovCurv, metname, SteppingHelixStateInfo::p3, SteppingHelixStateInfo::path(), SteppingHelixStateInfo::path_, SteppingHelixStateInfo::q, SteppingHelixStateInfo::r3, SteppingHelixStateInfo::radPath(), SteppingHelixStateInfo::radPath_, SteppingHelixStateInfo::radX0, SteppingHelixStateInfo::rzLims, AlCaHLTBitMon_QueryRunRegistry::string, metsig::tau, createJobs::tmp, useInTeslaFromMagField_, useIsYokeFlag_, useMagVolumes_, and vbField_.

Referenced by makeAtomStep().

782  {
783  static const std::string metname = "SteppingHelixPropagator";
784  svNext.q = svPrevious.q;
785  svNext.dir = dS > 0.0 ? 1. : -1.;
786  svNext.p3 = tau;
787  svNext.p3 *= (svPrevious.p3.mag() - svNext.dir * fabs(dP));
788 
789  svNext.r3 = svPrevious.r3;
790  svNext.r3 += drVec;
791 
792  svNext.path_ = svPrevious.path() + dS;
793  svNext.radPath_ = svPrevious.radPath() + dX0;
794 
795  GlobalPoint gPointNorZ(svNext.r3.x(), svNext.r3.y(), svNext.r3.z());
796 
797  GlobalVector bf(0, 0, 0);
798 
799  if (useMagVolumes_) {
800  if (vbField_ != nullptr) {
801  svNext.magVol = vbField_->findVolume(gPointNorZ);
802  if (useIsYokeFlag_) {
803  double curRad = svNext.r3.perp();
804  if (curRad > 380 && curRad < 850 && fabs(svNext.r3.z()) < 667) {
805  svNext.isYokeVol = isYokeVolume(svNext.magVol);
806  } else {
807  svNext.isYokeVol = false;
808  }
809  }
810  } else {
811  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
812  << "Failed to cast into VolumeBasedMagneticField" << std::endl;
813  svNext.magVol = nullptr;
814  }
815  if (debug_) {
816  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Got volume at "
817  << svNext.magVol << std::endl;
818  }
819  }
820 
821  if (useMagVolumes_ && svNext.magVol != nullptr && !useInTeslaFromMagField_) {
822  bf = svNext.magVol->inTesla(gPointNorZ);
823  svNext.bf.set(bf.x(), bf.y(), bf.z());
824  } else {
825  GlobalPoint gPoint(svNext.r3.x(), svNext.r3.y(), svNext.r3.z());
826  bf = field_->inTesla(gPoint);
827  svNext.bf.set(bf.x(), bf.y(), bf.z());
828  }
829  if (svNext.bf.mag2() < 1e-32)
830  svNext.bf.set(0., 0., 1e-16);
831 
832  double dEdXPrime = 0;
833  double dEdx = 0;
834  double radX0 = 0;
835  MatBounds rzLims;
836  dEdx = getDeDx(svNext, dEdXPrime, radX0, rzLims);
837  svNext.dEdx = dEdx;
838  svNext.dEdXPrime = dEdXPrime;
839  svNext.radX0 = radX0;
840  svNext.rzLims = rzLims;
841 
842  //update Emat only if it's valid
843  svNext.hasErrorPropagated_ = svPrevious.hasErrorPropagated_;
844  if (svPrevious.hasErrorPropagated_) {
845  {
846  AlgebraicMatrix55 tmp = dCovCurvTransform * svPrevious.covCurv;
847  ROOT::Math::AssignSym::Evaluate(svNext.covCurv, tmp * ROOT::Math::Transpose(dCovCurvTransform));
848 
849  svNext.covCurv += svPrevious.matDCovCurv;
850  }
851  } else {
852  //could skip dragging along the unprop. cov later .. now
853  // svNext.cov = svPrevious.cov;
854  }
855 
856  if (debug_) {
857  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Now at path: " << svNext.path()
858  << " radPath: " << svNext.radPath() << " p3 "
859  << " pt: " << svNext.p3.perp() << " phi: " << svNext.p3.phi() << " eta: " << svNext.p3.eta()
860  << " " << svNext.p3 << " r3: " << svNext.r3
861  << " dPhi: " << acos(svNext.p3.unit().dot(svPrevious.p3.unit())) << " bField: " << svNext.bf.mag()
862  << " dedx: " << svNext.dEdx << std::endl;
863  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "New Covariance "
864  << svNext.covCurv << std::endl;
865  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Transf by dCovTransform "
866  << dCovCurvTransform << std::endl;
867  }
868 }
MagVolume::inTesla
::GlobalVector inTesla(const ::GlobalPoint &gp) const override
Definition: MagVolume.h:42
metname
const std::string metname
Definition: MuonSeedOrcaPatternRecognition.cc:40
MagneticField::inTesla
virtual GlobalVector inTesla(const GlobalPoint &gp) const =0
Field value ad specified global point, in Tesla.
SteppingHelixPropagator::isYokeVolume
bool isYokeVolume(const MagVolume *vol) const
check if it&#39;s a yoke/iron based on this MagVol internals
Definition: SteppingHelixPropagator.cc:2432
SteppingHelixPropagator::MatBounds
SteppingHelixStateInfo::VolumeBounds MatBounds
Definition: SteppingHelixPropagator.h:170
SteppingHelixStateInfo::covCurv
AlgebraicSymMatrix55 covCurv
Definition: SteppingHelixStateInfo.h:87
AlgebraicMatrix55
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepStd< double, 5, 5 > > AlgebraicMatrix55
Definition: AlgebraicROOTObjects.h:55
VolumeBasedMagneticField::findVolume
const MagVolume * findVolume(const GlobalPoint &gp) const
Definition: VolumeBasedMagneticField.cc:62
SteppingHelixPropagator::getDeDx
double getDeDx(const SteppingHelixPropagator::StateInfo &sv, double &dEdXPrime, double &radX0, MatBounds &rzLims) const
estimate average (in fact smth. close to MPV and median) energy loss per unit path length ...
Definition: SteppingHelixPropagator.cc:1317
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
LogTrace
#define LogTrace(id)
Definition: MessageLogger.h:234
plot_hgcal_utils.dEdx
tuple dEdx
Definition: plot_hgcal_utils.py:216
SteppingHelixPropagator::vbField_
const VolumeBasedMagneticField * vbField_
Definition: SteppingHelixPropagator.h:266
SteppingHelixStateInfo::matDCovCurv
AlgebraicSymMatrix55 matDCovCurv
Definition: SteppingHelixStateInfo.h:88
createJobs.tmp
tmp
align.sh
Definition: createJobs.py:716
SteppingHelixPropagator::field_
const MagneticField * field_
Definition: SteppingHelixPropagator.h:265
void SteppingHelixPropagator::initStateArraySHPSpecific ( StateArray &  svBuf,
bool  flagsOnly 
) const
protected

(Internals) set defaults needed for states used inside propagate methods

Definition at line 38 of file SteppingHelixPropagator.cc.

References mps_fire::i, MAX_POINTS, and noErrorPropagation_.

Referenced by propagate(), and propagateWithPath().

38  {
39  for (int i = 0; i <= MAX_POINTS; i++) {
40  svBuf[i].isComplete = true;
41  svBuf[i].isValid_ = true;
42  svBuf[i].hasErrorPropagated_ = !noErrorPropagation_;
43  if (!flagsOnly) {
44  svBuf[i].p3 = Vector(0, 0, 0);
45  svBuf[i].r3 = Point(0, 0, 0);
46  svBuf[i].bf = Vector(0, 0, 0);
47  svBuf[i].bfGradLoc = Vector(0, 0, 0);
48  svBuf[i].covCurv = AlgebraicSymMatrix55();
49  svBuf[i].matDCovCurv = AlgebraicSymMatrix55();
50  }
51  }
52 }
mps_fire.i
i
Definition: mps_fire.py:428
AlgebraicSymMatrix55
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepSym< double, 5 > > AlgebraicSymMatrix55
Definition: AlgebraicROOTObjects.h:23
bool SteppingHelixPropagator::isYokeVolume ( const MagVolume *  vol) const
protected

check if it's a yoke/iron based on this MagVol internals

Definition at line 2432 of file SteppingHelixPropagator.cc.

References debug_, MagVolume::isIron(), LogTrace, metname, GloballyPositioned< T >::position(), mps_fire::result, and AlCaHLTBitMon_QueryRunRegistry::string.

Referenced by getNextState(), and loadState().

2432  {
2433  static const std::string metname = "SteppingHelixPropagator";
2434  if (vol == nullptr)
2435  return false;
2436  /*
2437  const MFGrid* mGrid = reinterpret_cast<const MFGrid*>(vol->provider());
2438  std::vector<int> dims(mGrid->dimensions());
2439 
2440  LocalVector lVCen(mGrid->nodeValue(dims[0]/2, dims[1]/2, dims[2]/2));
2441  LocalVector lVZLeft(mGrid->nodeValue(dims[0]/2, dims[1]/2, dims[2]/5));
2442  LocalVector lVZRight(mGrid->nodeValue(dims[0]/2, dims[1]/2, (dims[2]*4)/5));
2443 
2444  double mag2VCen = lVCen.mag2();
2445  double mag2VZLeft = lVZLeft.mag2();
2446  double mag2VZRight = lVZRight.mag2();
2447 
2448  bool result = false;
2449  if (mag2VCen > 0.6 && mag2VZLeft > 0.6 && mag2VZRight > 0.6){
2450  if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Volume is magnetic, located at "<<vol->position()<<std::endl;
2451  result = true;
2452  } else {
2453  if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Volume is not magnetic, located at "<<vol->position()<<std::endl;
2454  }
2455 
2456  */
2457  bool result = vol->isIron();
2458  if (result) {
2459  if (debug_)
2460  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
2461  << "Volume is magnetic, located at " << vol->position() << std::endl;
2462  } else {
2463  if (debug_)
2464  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
2465  << "Volume is not magnetic, located at " << vol->position() << std::endl;
2466  }
2467 
2468  return result;
2469 }
metname
const std::string metname
Definition: MuonSeedOrcaPatternRecognition.cc:40
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
LogTrace
#define LogTrace(id)
Definition: MessageLogger.h:234
mps_fire.result
tuple result
Definition: mps_fire.py:311
MagVolume::isIron
bool isIron() const
Temporary hack to pass information on material. Will eventually be replaced!
Definition: MagVolume.h:45
GloballyPositioned::position
const PositionType & position() const
Definition: GloballyPositioned.h:36
void SteppingHelixPropagator::loadState ( SteppingHelixPropagator::StateInfo &  svCurrent,
const SteppingHelixPropagator::Vector &  p3,
const SteppingHelixPropagator::Point &  r3,
int  charge,
PropagationDirection  dir,
const AlgebraicSymMatrix55 &  covCurv 
) const
protected

(Internals) compute transient values for initial point (resets step counter). Called by setIState

Definition at line 657 of file SteppingHelixPropagator.cc.

References alongMomentum, SteppingHelixStateInfo::bf, RecoTauCleanerPlugins::charge, SteppingHelixStateInfo::covCurv, debug_, SteppingHelixStateInfo::dEdx, plot_hgcal_utils::dEdx, SteppingHelixStateInfo::dEdXPrime, SteppingHelixStateInfo::dir, alignCSCRings::e, validate-o2o-wbm::f, field_, MagVolume::fieldInTesla(), VolumeBasedMagneticField::findVolume(), getDeDx(), MagneticField::inTesla(), MagVolume::inTesla(), SteppingHelixStateInfo::isComplete, SteppingHelixStateInfo::isValid_, SteppingHelixStateInfo::isYokeVol, isYokeVolume(), LogTrace, PV3DBase< T, PVType, FrameType >::mag2(), SteppingHelixStateInfo::magVol, metname, SteppingHelixStateInfo::p3, SteppingHelixStateInfo::path(), SteppingHelixStateInfo::path_, SteppingHelixStateInfo::q, SteppingHelixStateInfo::r3, SteppingHelixStateInfo::radPath(), SteppingHelixStateInfo::radPath_, SteppingHelixStateInfo::radX0, SteppingHelixStateInfo::rzLims, AlCaHLTBitMon_QueryRunRegistry::string, GloballyPositioned< T >::toLocal(), useInTeslaFromMagField_, useIsYokeFlag_, useMagVolumes_, vbField_, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by setIState().

662  {
663  static const std::string metname = "SteppingHelixPropagator";
664 
665  svCurrent.q = charge;
666  svCurrent.p3 = p3;
667  svCurrent.r3 = r3;
668  svCurrent.dir = dir == alongMomentum ? 1. : -1.;
669 
670  svCurrent.path_ = 0; // this could've held the initial path
671  svCurrent.radPath_ = 0;
672 
673  GlobalPoint gPointNorZ(svCurrent.r3.x(), svCurrent.r3.y(), svCurrent.r3.z());
674 
675  float gpmag = gPointNorZ.mag2();
676  float pmag2 = p3.mag2();
677  if (gpmag > 1e20f) {
678  LogTrace(metname) << "Initial point is too far";
679  svCurrent.isValid_ = false;
680  return;
681  }
682  if (pmag2 < 1e-18f) {
683  LogTrace(metname) << "Initial momentum is too low";
684  svCurrent.isValid_ = false;
685  return;
686  }
687  if (!(gpmag == gpmag)) {
688  LogTrace(metname) << "Initial point is a nan";
689  edm::LogWarning("SteppingHelixPropagatorNAN") << "Initial point is a nan";
690  svCurrent.isValid_ = false;
691  return;
692  }
693  if (!(pmag2 == pmag2)) {
694  LogTrace(metname) << "Initial momentum is a nan";
695  edm::LogWarning("SteppingHelixPropagatorNAN") << "Initial momentum is a nan";
696  svCurrent.isValid_ = false;
697  return;
698  }
699 
700  GlobalVector bf(0, 0, 0);
701  // = field_->inTesla(gPoint);
702  if (useMagVolumes_) {
703  if (vbField_) {
704  svCurrent.magVol = vbField_->findVolume(gPointNorZ);
705  if (useIsYokeFlag_) {
706  double curRad = svCurrent.r3.perp();
707  if (curRad > 380 && curRad < 850 && fabs(svCurrent.r3.z()) < 667) {
708  svCurrent.isYokeVol = isYokeVolume(svCurrent.magVol);
709  } else {
710  svCurrent.isYokeVol = false;
711  }
712  }
713  } else {
714  edm::LogWarning(metname) << std::setprecision(17) << std::setw(20) << std::scientific
715  << "Failed to cast into VolumeBasedMagneticField: fall back to the default behavior"
716  << std::endl;
717  svCurrent.magVol = nullptr;
718  }
719  if (debug_) {
720  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Got volume at "
721  << svCurrent.magVol << std::endl;
722  }
723  }
724 
725  if (useMagVolumes_ && svCurrent.magVol != nullptr && !useInTeslaFromMagField_) {
726  bf = svCurrent.magVol->inTesla(gPointNorZ);
727  if (debug_) {
728  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Loaded bfield float: " << bf
729  << " at global float " << gPointNorZ << " double " << svCurrent.r3 << std::endl;
730  LocalPoint lPoint(svCurrent.magVol->toLocal(gPointNorZ));
731  LocalVector lbf = svCurrent.magVol->fieldInTesla(lPoint);
732  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "\t cf in local locF: " << lbf
733  << " at " << lPoint << std::endl;
734  }
735  svCurrent.bf.set(bf.x(), bf.y(), bf.z());
736  } else {
737  GlobalPoint gPoint(r3.x(), r3.y(), r3.z());
738  bf = field_->inTesla(gPoint);
739  svCurrent.bf.set(bf.x(), bf.y(), bf.z());
740  }
741  if (svCurrent.bf.mag2() < 1e-32)
742  svCurrent.bf.set(0., 0., 1e-16);
743  if (debug_) {
744  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
745  << "Loaded bfield double: " << svCurrent.bf << " from float: " << bf << " at float "
746  << gPointNorZ << " double " << svCurrent.r3 << std::endl;
747  }
748 
749  double dEdXPrime = 0;
750  double dEdx = 0;
751  double radX0 = 0;
752  MatBounds rzLims;
753  dEdx = getDeDx(svCurrent, dEdXPrime, radX0, rzLims);
754  svCurrent.dEdx = dEdx;
755  svCurrent.dEdXPrime = dEdXPrime;
756  svCurrent.radX0 = radX0;
757  svCurrent.rzLims = rzLims;
758 
759  svCurrent.covCurv = covCurv;
760 
761  svCurrent.isComplete = true;
762  svCurrent.isValid_ = true;
763 
764  if (debug_) {
765  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
766  << "Loaded at path: " << svCurrent.path() << " radPath: " << svCurrent.radPath() << " p3 "
767  << " pt: " << svCurrent.p3.perp() << " phi: " << svCurrent.p3.phi()
768  << " eta: " << svCurrent.p3.eta() << " " << svCurrent.p3 << " r3: " << svCurrent.r3
769  << " bField: " << svCurrent.bf.mag() << std::endl;
770  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
771  << "Input Covariance in Curvilinear RF " << covCurv << std::endl;
772  }
773 }
PV3DBase::mag2
T mag2() const
Definition: PV3DBase.h:63
MagVolume::inTesla
::GlobalVector inTesla(const ::GlobalPoint &gp) const override
Definition: MagVolume.h:42
metname
const std::string metname
Definition: MuonSeedOrcaPatternRecognition.cc:40
MagneticField::inTesla
virtual GlobalVector inTesla(const GlobalPoint &gp) const =0
Field value ad specified global point, in Tesla.
SteppingHelixPropagator::isYokeVolume
bool isYokeVolume(const MagVolume *vol) const
check if it&#39;s a yoke/iron based on this MagVol internals
Definition: SteppingHelixPropagator.cc:2432
SteppingHelixPropagator::MatBounds
SteppingHelixStateInfo::VolumeBounds MatBounds
Definition: SteppingHelixPropagator.h:170
SteppingHelixStateInfo::covCurv
AlgebraicSymMatrix55 covCurv
Definition: SteppingHelixStateInfo.h:87
VolumeBasedMagneticField::findVolume
const MagVolume * findVolume(const GlobalPoint &gp) const
Definition: VolumeBasedMagneticField.cc:62
SteppingHelixPropagator::getDeDx
double getDeDx(const SteppingHelixPropagator::StateInfo &sv, double &dEdXPrime, double &radX0, MatBounds &rzLims) const
estimate average (in fact smth. close to MPV and median) energy loss per unit path length ...
Definition: SteppingHelixPropagator.cc:1317
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
LogTrace
#define LogTrace(id)
Definition: MessageLogger.h:234
plot_hgcal_utils.dEdx
tuple dEdx
Definition: plot_hgcal_utils.py:216
DeadROC_duringRun.dir
string dir
Definition: DeadROC_duringRun.py:23
GloballyPositioned::toLocal
LocalPoint toLocal(const GlobalPoint &gp) const
Definition: GloballyPositioned.h:98
SteppingHelixPropagator::vbField_
const VolumeBasedMagneticField * vbField_
Definition: SteppingHelixPropagator.h:266
MagVolume::fieldInTesla
LocalVector fieldInTesla(const LocalPoint &lp) const
Definition: MagVolume.cc:11
edm::LogWarning
Log< level::Warning, false > LogWarning
Definition: MessageLogger.h:122
SteppingHelixPropagator::field_
const MagneticField * field_
Definition: SteppingHelixPropagator.h:265
const MagneticField* SteppingHelixPropagator::magneticField ( ) const
inlineoverridevirtual

Implements Propagator.

Definition at line 81 of file SteppingHelixPropagator.h.

References field_.

81 { return field_; }
SteppingHelixPropagator::field_
const MagneticField * field_
Definition: SteppingHelixPropagator.h:265
bool SteppingHelixPropagator::makeAtomStep ( SteppingHelixPropagator::StateInfo &  svCurrent,
SteppingHelixPropagator::StateInfo &  svNext,
double  dS,
PropagationDirection  dir,
SteppingHelixPropagator::Fancy  fancy 
) const
protected

main stepping function: compute next state vector after a step of length dS

Definition at line 879 of file SteppingHelixPropagator.cc.

References alongMomentum, applyRadX0Correction_, b0, SteppingHelixStateInfo::bf, funct::cos(), Phase1L1TJetProducer_cfi::cosPhi, debug_, SteppingHelixStateInfo::dEdx, plot_hgcal_utils::dEdx, SteppingHelixStateInfo::dEdXPrime, PVValHelper::dx, alignCSCRings::e, field_, getDeDx(), getNextState(), SteppingHelixStateInfo::hasErrorPropagated_, HEL_ALL_F, HEL_AS_F, MagneticField::inTesla(), MagVolume::inTesla(), SteppingHelixStateInfo::isYokeVol, log, LogTrace, SteppingHelixPropagator::Basis::lX, SteppingHelixPropagator::Basis::lY, SteppingHelixPropagator::Basis::lZ, SteppingHelixStateInfo::magVol, SteppingHelixStateInfo::matDCovCurv, metname, oppositeToMomentum, SteppingHelixStateInfo::p3, SteppingHelixStateInfo::path(), phi, createTree::pp, submitPVResolutionJobs::q, SteppingHelixStateInfo::q, SteppingHelixStateInfo::r3, SteppingHelixStateInfo::radPath(), SteppingHelixStateInfo::radX0, cuy::rep, setRep(), funct::sin(), Phase1L1TJetProducer_cfi::sinPhi, mathSSE::sqrt(), AlCaHLTBitMon_QueryRunRegistry::string, metsig::tau, unit55_, useInTeslaFromMagField_, useMagVolumes_, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by propagate().

883  {
884  static const std::string metname = "SteppingHelixPropagator";
885  if (debug_) {
886  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Make atom step "
887  << svCurrent.path() << " with step " << dS << " in direction " << dir << std::endl;
888  }
889 
890  AlgebraicMatrix55 dCCurvTransform(unit55_);
891 
892  double dP = 0;
893  double curP = svCurrent.p3.mag();
894  Vector tau = svCurrent.p3;
895  tau *= 1. / curP;
896  Vector tauNext(tau);
897  Vector drVec(0, 0, 0);
898 
899  dS = dir == alongMomentum ? fabs(dS) : -fabs(dS);
900 
901  double radX0 = 1e24;
902 
903  switch (fancy) {
904  case HEL_AS_F:
905  case HEL_ALL_F: {
906  double p0 = curP; // see above = svCurrent.p3.mag();
907  double p0Inv = 1. / p0;
908  double b0 = svCurrent.bf.mag();
909 
910  //get to the mid-point first
911  double phi = (2.99792458e-3 * svCurrent.q * b0 * p0Inv) * dS / 2.;
912  bool phiSmall = fabs(phi) < 1e-4;
913 
914  double cosPhi = 0;
915  double sinPhi = 0;
916 
917  double oneLessCosPhi = 0;
918  double oneLessCosPhiOPhi = 0;
919  double phiLessSinPhiOPhi = 0;
920 
921  if (phiSmall) {
922  double phi2 = phi * phi;
923  double phi3 = phi2 * phi;
924  double phi4 = phi3 * phi;
925  oneLessCosPhi = phi2 / 2. - phi4 / 24. + phi2 * phi4 / 720.; // 0.5*phi*phi;//*(1.- phi*phi/12.);
926  oneLessCosPhiOPhi = 0.5 * phi - phi3 / 24. + phi2 * phi3 / 720.; //*(1.- phi*phi/12.);
927  phiLessSinPhiOPhi = phi * phi / 6. - phi4 / 120. + phi4 * phi2 / 5040.; //*(1. - phi*phi/20.);
928  } else {
929  cosPhi = cos(phi);
930  sinPhi = sin(phi);
931  oneLessCosPhi = 1. - cosPhi;
932  oneLessCosPhiOPhi = oneLessCosPhi / phi;
933  double sinPhiOPhi = sinPhi / phi;
934  phiLessSinPhiOPhi = 1 - sinPhiOPhi;
935  }
936 
937  Vector bHat = svCurrent.bf;
938  bHat *= 1. / b0; //bHat.mag();
939  // bool bAlongZ = fabs(bHat.z()) > 0.9999;
940 
941  Vector btVec(bHat.cross(tau)); // for balong z btVec.z()==0
942  double tauB = tau.dot(bHat);
943  Vector bbtVec(bHat * tauB - tau); // (-tau.x(), -tau.y(), 0)
944 
945  //don't need it here tauNext = tau + bbtVec*oneLessCosPhi - btVec*sinPhi;
946  drVec = bbtVec * phiLessSinPhiOPhi;
947  drVec -= btVec * oneLessCosPhiOPhi;
948  drVec += tau;
949  drVec *= dS / 2.;
950 
951  double dEdx = svCurrent.dEdx;
952  double dEdXPrime = svCurrent.dEdXPrime;
953  radX0 = svCurrent.radX0;
954  dP = dEdx * dS;
955 
956  //improve with above values:
957  drVec += svCurrent.r3;
958  GlobalVector bfGV(0, 0, 0);
959  Vector bf(0, 0, 0);
960  if (useMagVolumes_ && svCurrent.magVol != nullptr && !useInTeslaFromMagField_) {
961  bfGV = svCurrent.magVol->inTesla(GlobalPoint(drVec.x(), drVec.y(), drVec.z()));
962  bf.set(bfGV.x(), bfGV.y(), bfGV.z());
963  } else {
964  bfGV = field_->inTesla(GlobalPoint(drVec.x(), drVec.y(), drVec.z()));
965  bf.set(bfGV.x(), bfGV.y(), bfGV.z());
966  }
967  double b0Init = b0;
968  b0 = bf.mag();
969  if (b0 < 1e-16) {
970  b0 = 1e-16;
971  bf.set(0., 0., 1e-16);
972  }
973  if (debug_) {
974  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Improved b " << b0
975  << " at r3 " << drVec << std::endl;
976  }
977 
978  if (fabs((b0 - b0Init) * dS) > 1) {
979  //missed the mag volume boundary?
980  if (debug_) {
981  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Large bf*dS change "
982  << fabs((b0 - svCurrent.bf.mag()) * dS) << " --> recalc dedx" << std::endl;
983  }
984  svNext.r3 = drVec;
985  svNext.bf = bf;
986  svNext.p3 = svCurrent.p3;
987  svNext.isYokeVol = svCurrent.isYokeVol;
988  svNext.magVol = svCurrent.magVol;
989  MatBounds rzTmp;
990  dEdx = getDeDx(svNext, dEdXPrime, radX0, rzTmp);
991  dP = dEdx * dS;
992  if (debug_) {
993  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "New dEdX= " << dEdx
994  << " dP= " << dP << " for p0 " << p0 << std::endl;
995  }
996  }
997  //p0 is mid-way and b0 from mid-point
998  p0 += dP / 2.;
999  if (p0 < 1e-2)
1000  p0 = 1e-2;
1001  p0Inv = 1. / p0;
1002 
1003  phi = (2.99792458e-3 * svCurrent.q * b0 * p0Inv) * dS;
1004  phiSmall = fabs(phi) < 1e-4;
1005 
1006  if (phiSmall) {
1007  double phi2 = phi * phi;
1008  double phi3 = phi2 * phi;
1009  double phi4 = phi3 * phi;
1010  sinPhi = phi - phi3 / 6. + phi4 * phi / 120.;
1011  cosPhi = 1. - phi2 / 2. + phi4 / 24.;
1012  oneLessCosPhi = phi2 / 2. - phi4 / 24. + phi2 * phi4 / 720.; // 0.5*phi*phi;//*(1.- phi*phi/12.);
1013  oneLessCosPhiOPhi = 0.5 * phi - phi3 / 24. + phi2 * phi3 / 720.; //*(1.- phi*phi/12.);
1014  phiLessSinPhiOPhi = phi * phi / 6. - phi4 / 120. + phi4 * phi2 / 5040.; //*(1. - phi*phi/20.);
1015  } else {
1016  cosPhi = cos(phi);
1017  sinPhi = sin(phi);
1018  oneLessCosPhi = 1. - cosPhi;
1019  oneLessCosPhiOPhi = oneLessCosPhi / phi;
1020  double sinPhiOPhi = sinPhi / phi;
1021  phiLessSinPhiOPhi = 1. - sinPhiOPhi;
1022  }
1023 
1024  bHat = bf;
1025  bHat *= 1. / b0; // as above =1./bHat.mag();
1026  // bAlongZ = fabs(bHat.z()) > 0.9999;
1027  btVec = bHat.cross(tau); // for b||z (-tau.y(), tau.x() ,0)
1028  tauB = tau.dot(bHat);
1029  bbtVec = bHat * tauB - tau; //bHat.cross(btVec); for b||z: (-tau.x(), -tau.y(), 0)
1030 
1031  tauNext = bbtVec * oneLessCosPhi;
1032  tauNext -= btVec * sinPhi;
1033  tauNext += tau; //for b||z tauNext.z() == tau.z()
1034  double tauNextPerpInv = 1. / tauNext.perp();
1035  drVec = bbtVec * phiLessSinPhiOPhi;
1036  drVec -= btVec * oneLessCosPhiOPhi;
1037  drVec += tau;
1038  drVec *= dS;
1039 
1040  if (svCurrent.hasErrorPropagated_) {
1041  double theta02 = 0;
1042  double dX0 = fabs(dS) / radX0;
1043 
1044  if (applyRadX0Correction_) {
1045  // this provides the integrand for theta^2
1046  // if summed up along the path, should result in
1047  // theta_total^2 = Int_0^x0{ f(x)dX} = (13.6/p0)^2*x0*(1+0.036*ln(x0+1))
1048  // x0+1 above is to make the result infrared safe.
1049  double x0 = fabs(svCurrent.radPath());
1050  double alphaX0 = 13.6e-3 * p0Inv;
1051  alphaX0 *= alphaX0;
1052  double betaX0 = 0.038;
1053  double logx0p1 = log(x0 + 1);
1054  theta02 = dX0 * alphaX0 * (1 + betaX0 * logx0p1) * (1 + betaX0 * logx0p1 + 2. * betaX0 * x0 / (x0 + 1));
1055  } else {
1056  theta02 = 196e-6 * p0Inv * p0Inv *
1057  dX0; //14.e-3/p0*sqrt(fabs(dS)/radX0); // .. drop log term (this is non-additive)
1058  }
1059 
1060  double epsilonP0 = 1. + dP / (p0 - 0.5 * dP);
1061  // double omegaP0 = -dP/(p0-0.5*dP) + dS*dEdXPrime;
1062  // double dsp = dS/(p0-0.5*dP); //use the initial p0 (not the mid-point) to keep the transport properly additive
1063 
1064  Vector tbtVec(bHat - tauB * tau); // for b||z tau.z()*(-tau.x(), -tau.y(), 1.-tau.z())
1065 
1066  {
1067  //Slightly modified copy of the curvilinear jacobian (don't use the original just because it's in float precision
1068  // and seems to have some assumptions about the field values
1069  // notation changes: p1--> tau, p2-->tauNext
1070  // theta --> phi
1071  // Vector p1 = tau;
1072  // Vector p2 = tauNext;
1073  Point xStart = svCurrent.r3;
1074  const Vector& dx = drVec;
1075  //GlobalVector h = MagneticField::inInverseGeV(xStart);
1076  // Martijn: field is now given as parameter.. GlobalVector h = globalParameters.magneticFieldInInverseGeV(xStart);
1077 
1078  //double qbp = fts.signedInverseMomentum();
1079  double qbp = svCurrent.q * p0Inv;
1080  // double absS = dS;
1081 
1082  // calculate transport matrix
1083  // Origin: TRPRFN
1084  double t11 = tau.x();
1085  double t12 = tau.y();
1086  double t13 = tau.z();
1087  double t21 = tauNext.x();
1088  double t22 = tauNext.y();
1089  double t23 = tauNext.z();
1090  double cosl0 = tau.perp();
1091  // double cosl1 = 1./tauNext.perp(); //not quite a cos .. it's a cosec--> change to cosecl1 below
1092  double cosecl1 = tauNextPerpInv;
1093  //AlgebraicMatrix a(5,5,1);
1094  // define average magnetic field and gradient
1095  // at initial point - inlike TRPRFN
1096  const Vector& hn = bHat;
1097  // double qp = -2.99792458e-3*b0;
1098  // double q = -h.mag()*qbp;
1099 
1100  double q = -phi / dS; //qp*qbp; // -phi/dS
1101  // double theta = -phi;
1102  double sint = -sinPhi;
1103  double cost = cosPhi;
1104  double hn1 = hn.x();
1105  double hn2 = hn.y();
1106  double hn3 = hn.z();
1107  double dx1 = dx.x();
1108  double dx2 = dx.y();
1109  double dx3 = dx.z();
1110  // double hnDt1 = hn1*t11 + hn2*t12 + hn3*t13;
1111 
1112  double gamma = hn1 * t21 + hn2 * t22 + hn3 * t23;
1113  double an1 = hn2 * t23 - hn3 * t22;
1114  double an2 = hn3 * t21 - hn1 * t23;
1115  double an3 = hn1 * t22 - hn2 * t21;
1116  // double auInv = sqrt(1.- t13*t13); double au = auInv>0 ? 1./auInv : 1e24;
1117  double auInv = cosl0;
1118  double au = auInv > 0 ? 1. / auInv : 1e24;
1119  // double auInv = sqrt(t11*t11 + t12*t12); double au = auInv>0 ? 1./auInv : 1e24;
1120  double u11 = -au * t12;
1121  double u12 = au * t11;
1122  double v11 = -t13 * u12;
1123  double v12 = t13 * u11;
1124  double v13 = auInv; //t11*u12 - t12*u11;
1125  auInv = sqrt(1. - t23 * t23);
1126  au = auInv > 0 ? 1. / auInv : 1e24;
1127  // auInv = sqrt(t21*t21 + t22*t22); au = auInv>0 ? 1./auInv : 1e24;
1128  double u21 = -au * t22;
1129  double u22 = au * t21;
1130  double v21 = -t23 * u22;
1131  double v22 = t23 * u21;
1132  double v23 = auInv; //t21*u22 - t22*u21;
1133  // now prepare the transport matrix
1134  // pp only needed in high-p case (WA)
1135  // double pp = 1./qbp;
1137  // moved up (where -h.mag() is needed()
1138  // double qp = q*pp;
1139  double anv = -(hn1 * u21 + hn2 * u22);
1140  double anu = (hn1 * v21 + hn2 * v22 + hn3 * v23);
1141  double omcost = oneLessCosPhi;
1142  double tmsint = -phi * phiLessSinPhiOPhi;
1143 
1144  double hu1 = -hn3 * u12;
1145  double hu2 = hn3 * u11;
1146  double hu3 = hn1 * u12 - hn2 * u11;
1147 
1148  double hv1 = hn2 * v13 - hn3 * v12;
1149  double hv2 = hn3 * v11 - hn1 * v13;
1150  double hv3 = hn1 * v12 - hn2 * v11;
1151 
1152  // 1/p - doesn't change since |tau| = |tauNext| ... not. It changes now
1153  dCCurvTransform(0, 0) = 1. / (epsilonP0 * epsilonP0) * (1. + dS * dEdXPrime);
1154 
1155  // lambda
1156 
1157  dCCurvTransform(1, 0) = phi * p0 / svCurrent.q * cosecl1 * (sinPhi * bbtVec.z() - cosPhi * btVec.z());
1158  //was dCCurvTransform(1,0) = -qp*anv*(t21*dx1 + t22*dx2 + t23*dx3); //NOTE (SK) this was found to have an opposite sign
1159  //from independent re-calculation ... in fact the tauNext.dot.dR piece isnt reproduced
1160 
1161  dCCurvTransform(1, 1) =
1162  cost * (v11 * v21 + v12 * v22 + v13 * v23) + sint * (hv1 * v21 + hv2 * v22 + hv3 * v23) +
1163  omcost * (hn1 * v11 + hn2 * v12 + hn3 * v13) * (hn1 * v21 + hn2 * v22 + hn3 * v23) +
1164  anv * (-sint * (v11 * t21 + v12 * t22 + v13 * t23) + omcost * (v11 * an1 + v12 * an2 + v13 * an3) -
1165  tmsint * gamma * (hn1 * v11 + hn2 * v12 + hn3 * v13));
1166 
1167  dCCurvTransform(1, 2) = cost * (u11 * v21 + u12 * v22) + sint * (hu1 * v21 + hu2 * v22 + hu3 * v23) +
1168  omcost * (hn1 * u11 + hn2 * u12) * (hn1 * v21 + hn2 * v22 + hn3 * v23) +
1169  anv * (-sint * (u11 * t21 + u12 * t22) + omcost * (u11 * an1 + u12 * an2) -
1170  tmsint * gamma * (hn1 * u11 + hn2 * u12));
1171  dCCurvTransform(1, 2) *= cosl0;
1172 
1173  // Commented out in part for reproducibility purposes: these terms are zero in cart->curv
1174  // dCCurvTransform(1,3) = -q*anv*(u11*t21 + u12*t22 ); //don't show up in cartesian setup-->curv
1175  //why would lambdaNext depend explicitely on initial position ? any arbitrary init point can be chosen not
1176  // affecting the final state's momentum direction ... is this the field gradient in curvilinear coord?
1177  // dCCurvTransform(1,4) = -q*anv*(v11*t21 + v12*t22 + v13*t23); //don't show up in cartesian setup-->curv
1178 
1179  // phi
1180 
1181  dCCurvTransform(2, 0) = -phi * p0 / svCurrent.q * cosecl1 * cosecl1 *
1182  (oneLessCosPhi * bHat.z() * btVec.mag2() + sinPhi * btVec.z() + cosPhi * tbtVec.z());
1183  //was dCCurvTransform(2,0) = -qp*anu*(t21*dx1 + t22*dx2 + t23*dx3)*cosecl1;
1184 
1185  dCCurvTransform(2, 1) =
1186  cost * (v11 * u21 + v12 * u22) + sint * (hv1 * u21 + hv2 * u22) +
1187  omcost * (hn1 * v11 + hn2 * v12 + hn3 * v13) * (hn1 * u21 + hn2 * u22) +
1188  anu * (-sint * (v11 * t21 + v12 * t22 + v13 * t23) + omcost * (v11 * an1 + v12 * an2 + v13 * an3) -
1189  tmsint * gamma * (hn1 * v11 + hn2 * v12 + hn3 * v13));
1190  dCCurvTransform(2, 1) *= cosecl1;
1191 
1192  dCCurvTransform(2, 2) = cost * (u11 * u21 + u12 * u22) + sint * (hu1 * u21 + hu2 * u22) +
1193  omcost * (hn1 * u11 + hn2 * u12) * (hn1 * u21 + hn2 * u22) +
1194  anu * (-sint * (u11 * t21 + u12 * t22) + omcost * (u11 * an1 + u12 * an2) -
1195  tmsint * gamma * (hn1 * u11 + hn2 * u12));
1196  dCCurvTransform(2, 2) *= cosecl1 * cosl0;
1197 
1198  // Commented out in part for reproducibility purposes: these terms are zero in cart->curv
1199  // dCCurvTransform(2,3) = -q*anu*(u11*t21 + u12*t22 )*cosecl1;
1200  //why would lambdaNext depend explicitely on initial position ? any arbitrary init point can be chosen not
1201  // affecting the final state's momentum direction ... is this the field gradient in curvilinear coord?
1202  // dCCurvTransform(2,4) = -q*anu*(v11*t21 + v12*t22 + v13*t23)*cosecl1;
1203 
1204  // yt
1205 
1206  double pp = 1. / qbp;
1207  // (SK) these terms seem to consistently have a sign opp from private derivation
1208  dCCurvTransform(3, 0) = -pp * (u21 * dx1 + u22 * dx2); //NB: modified from the original: changed the sign
1209  dCCurvTransform(4, 0) = -pp * (v21 * dx1 + v22 * dx2 + v23 * dx3);
1210 
1211  dCCurvTransform(3, 1) = (sint * (v11 * u21 + v12 * u22) + omcost * (hv1 * u21 + hv2 * u22) +
1212  tmsint * (hn1 * u21 + hn2 * u22) * (hn1 * v11 + hn2 * v12 + hn3 * v13)) /
1213  q;
1214 
1215  dCCurvTransform(3, 2) = (sint * (u11 * u21 + u12 * u22) + omcost * (hu1 * u21 + hu2 * u22) +
1216  tmsint * (hn1 * u21 + hn2 * u22) * (hn1 * u11 + hn2 * u12)) *
1217  cosl0 / q;
1218 
1219  dCCurvTransform(3, 3) = (u11 * u21 + u12 * u22);
1220 
1221  dCCurvTransform(3, 4) = (v11 * u21 + v12 * u22);
1222 
1223  // zt
1224 
1225  dCCurvTransform(4, 1) =
1226  (sint * (v11 * v21 + v12 * v22 + v13 * v23) + omcost * (hv1 * v21 + hv2 * v22 + hv3 * v23) +
1227  tmsint * (hn1 * v21 + hn2 * v22 + hn3 * v23) * (hn1 * v11 + hn2 * v12 + hn3 * v13)) /
1228  q;
1229 
1230  dCCurvTransform(4, 2) = (sint * (u11 * v21 + u12 * v22) + omcost * (hu1 * v21 + hu2 * v22 + hu3 * v23) +
1231  tmsint * (hn1 * v21 + hn2 * v22 + hn3 * v23) * (hn1 * u11 + hn2 * u12)) *
1232  cosl0 / q;
1233 
1234  dCCurvTransform(4, 3) = (u11 * v21 + u12 * v22);
1235 
1236  dCCurvTransform(4, 4) = (v11 * v21 + v12 * v22 + v13 * v23);
1237  // end of TRPRFN
1238  }
1239 
1240  if (debug_) {
1241  Basis rep;
1242  setRep(rep, tauNext);
1243  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "rep X: " << rep.lX << " "
1244  << rep.lX.mag() << "\t Y: " << rep.lY << " " << rep.lY.mag() << "\t Z: " << rep.lZ << " "
1245  << rep.lZ.mag();
1246  }
1247  //mind the sign of dS and dP (dS*dP < 0 allways)
1248  //covariance should grow no matter which direction you propagate
1249  //==> take abs values.
1250  //reset not needed: fill all below svCurrent.matDCov *= 0.;
1251  double mulRR = theta02 * dS * dS / 3.;
1252  double mulRP = theta02 * fabs(dS) * p0 / 2.;
1253  double mulPP = theta02 * p0 * p0;
1254  double losPP = dP * dP * 1.6 / fabs(dS) * (1.0 + p0 * 1e-3);
1255  //another guess .. makes sense for 1 cm steps 2./dS == 2 [cm] / dS [cm] at low pt
1256  //double it by 1TeV
1257  //not gaussian anyways
1258  // derived from the fact that sigma_p/eLoss ~ 0.08 after ~ 200 steps
1259 
1260  //curvilinear
1261  double sinThetaInv = tauNextPerpInv;
1262  double p0Mat =
1263  p0 + 0.5 * dP; // FIXME change this to p0 after it's clear that there's agreement in everything else
1264  double p0Mat2 = p0Mat * p0Mat;
1265  // with 6x6 formulation
1266  svCurrent.matDCovCurv *= 0;
1267 
1268  svCurrent.matDCovCurv(0, 0) = losPP / (p0Mat2 * p0Mat2);
1269  // svCurrent.matDCovCurv(0,1) = 0;
1270  // svCurrent.matDCovCurv(0,2) = 0;
1271  // svCurrent.matDCovCurv(0,3) = 0;
1272  // svCurrent.matDCovCurv(0,4) = 0;
1273 
1274  // svCurrent.matDCovCurv(1,0) = 0;
1275  svCurrent.matDCovCurv(1, 1) = mulPP / p0Mat2;
1276  // svCurrent.matDCovCurv(1,2) = 0;
1277  // svCurrent.matDCovCurv(1,3) = 0;
1278  svCurrent.matDCovCurv(1, 4) = mulRP / p0Mat;
1279 
1280  // svCurrent.matDCovCurv(2,0) = 0;
1281  // svCurrent.matDCovCurv(2,1) = 0;
1282  svCurrent.matDCovCurv(2, 2) = mulPP / p0Mat2 * (sinThetaInv * sinThetaInv);
1283  svCurrent.matDCovCurv(2, 3) = mulRP / p0Mat * sinThetaInv;
1284  // svCurrent.matDCovCurv(2,4) = 0;
1285 
1286  // svCurrent.matDCovCurv(3,0) = 0;
1287  // svCurrent.matDCovCurv(3,1) = 0;
1288  svCurrent.matDCovCurv(3, 2) = mulRP / p0Mat * sinThetaInv;
1289  // svCurrent.matDCovCurv(3,0) = 0;
1290  svCurrent.matDCovCurv(3, 3) = mulRR;
1291  // svCurrent.matDCovCurv(3,4) = 0;
1292 
1293  // svCurrent.matDCovCurv(4,0) = 0;
1294  svCurrent.matDCovCurv(4, 1) = mulRP / p0Mat;
1295  // svCurrent.matDCovCurv(4,2) = 0;
1296  // svCurrent.matDCovCurv(4,3) = 0;
1297  svCurrent.matDCovCurv(4, 4) = mulRR;
1298  }
1299  break;
1300  }
1301  default:
1302  break;
1303  }
1304 
1305  double pMag = curP; //svCurrent.p3.mag();
1306 
1307  if (dir == oppositeToMomentum)
1308  dP = fabs(dP);
1309  else if (dP < 0) { //the case of negative dP
1310  dP = -dP > pMag ? -pMag + 1e-5 : dP;
1311  }
1312 
1313  getNextState(svCurrent, svNext, dP, tauNext, drVec, dS, dS / radX0, dCCurvTransform);
1314  return true;
1315 }
cuy.rep
string rep
Definition: cuy.py:1189
log
static std::vector< std::string > checklist log
Definition: LHERunInfoProduct.cc:179
createTree.pp
tuple pp
Definition: createTree.py:17
MagVolume::inTesla
::GlobalVector inTesla(const ::GlobalPoint &gp) const override
Definition: MagVolume.h:42
metname
const std::string metname
Definition: MuonSeedOrcaPatternRecognition.cc:40
Vector
ROOT::Math::Plane3D::Vector Vector
Definition: EcalHitMaker.cc:29
MagneticField::inTesla
virtual GlobalVector inTesla(const GlobalPoint &gp) const =0
Field value ad specified global point, in Tesla.
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
GlobalPoint
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
SteppingHelixPropagator::MatBounds
SteppingHelixStateInfo::VolumeBounds MatBounds
Definition: SteppingHelixPropagator.h:170
AlgebraicMatrix55
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepStd< double, 5, 5 > > AlgebraicMatrix55
Definition: AlgebraicROOTObjects.h:55
SteppingHelixPropagator::getDeDx
double getDeDx(const SteppingHelixPropagator::StateInfo &sv, double &dEdXPrime, double &radX0, MatBounds &rzLims) const
estimate average (in fact smth. close to MPV and median) energy loss per unit path length ...
Definition: SteppingHelixPropagator.cc:1317
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
LogTrace
#define LogTrace(id)
Definition: MessageLogger.h:234
plot_hgcal_utils.dEdx
tuple dEdx
Definition: plot_hgcal_utils.py:216
DeadROC_duringRun.dir
string dir
Definition: DeadROC_duringRun.py:23
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
Phase1L1TJetProducer_cfi.cosPhi
cosPhi
Definition: Phase1L1TJetProducer_cfi.py:40
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
SteppingHelixStateInfo::matDCovCurv
AlgebraicSymMatrix55 matDCovCurv
Definition: SteppingHelixStateInfo.h:88
Phase1L1TJetProducer_cfi.sinPhi
sinPhi
Definition: Phase1L1TJetProducer_cfi.py:39
SteppingHelixPropagator::unit55_
const AlgebraicSymMatrix55 unit55_
Definition: SteppingHelixPropagator.h:267
SteppingHelixPropagator::setRep
void setRep(SteppingHelixPropagator::Basis &rep, const SteppingHelixPropagator::Vector &tau) const
Set/compute basis vectors for local coordinates at current step (called by incrementState) ...
Definition: SteppingHelixPropagator.cc:870
Point
Structure Point Contains parameters of Gaussian fits to DMRs.
Definition: DMRtrends.cc:57
SteppingHelixPropagator::getNextState
void getNextState(const SteppingHelixPropagator::StateInfo &svPrevious, SteppingHelixPropagator::StateInfo &svNext, double dP, const SteppingHelixPropagator::Vector &tau, const SteppingHelixPropagator::Vector &drVec, double dS, double dX0, const AlgebraicMatrix55 &dCovCurv) const
Definition: SteppingHelixPropagator.cc:775
b0
static constexpr float b0
Definition: L1EGammaCrystalsEmulatorProducer.cc:83
SteppingHelixPropagator::field_
const MagneticField * field_
Definition: SteppingHelixPropagator.h:265
void SteppingHelixPropagator::propagate ( const SteppingHelixStateInfo &  ftsStart,
const Surface &  sDest,
SteppingHelixStateInfo &  out 
) const

Propagate to Plane given a starting point.

Definition at line 152 of file SteppingHelixPropagator.cc.

References invalidState_, SteppingHelixStateInfo::isValid(), and sendLogWarning_.

Referenced by AlCaHOCalibProducer::fillHOStore(), propagate(), and propagateWithPath().

154  {
155  if (!sStart.isValid()) {
156  if (sendLogWarning_) {
157  edm::LogWarning("SteppingHelixPropagator") << "Can't propagate: invalid input state" << std::endl;
158  }
159  out = invalidState_;
160  return;
161  }
162 
163  const Plane* pDest = dynamic_cast<const Plane*>(&sDest);
164  if (pDest != nullptr) {
165  propagate(sStart, *pDest, out);
166  return;
167  }
168 
169  const Cylinder* cDest = dynamic_cast<const Cylinder*>(&sDest);
170  if (cDest != nullptr) {
171  propagate(sStart, *cDest, out);
172  return;
173  }
174 
175  throw PropagationException("The surface is neither Cylinder nor Plane");
176 }
Plane
Definition: Plane.h:16
SteppingHelixPropagator::propagate
void propagate(const SteppingHelixStateInfo &ftsStart, const Surface &sDest, SteppingHelixStateInfo &out) const
Propagate to Plane given a starting point.
Definition: SteppingHelixPropagator.cc:152
PropagationException
Common base class.
Definition: PropagationExceptions.h:14
edm::LogWarning
Log< level::Warning, false > LogWarning
Definition: MessageLogger.h:122
void SteppingHelixPropagator::propagate ( const SteppingHelixStateInfo &  ftsStart,
const Plane &  pDest,
SteppingHelixStateInfo &  out 
) const

Definition at line 178 of file SteppingHelixPropagator.cc.

References cIndex_(), initStateArraySHPSpecific(), invalidState_, SteppingHelixStateInfo::isValid(), SteppingHelixStateInfo::path(), PLANE_DT, GloballyPositioned< T >::position(), propagate(), GloballyPositioned< T >::rotation(), sendLogWarning_, setIState(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), PV3DBase< T, PVType, FrameType >::z(), TkRotation< T >::zx(), TkRotation< T >::zy(), and TkRotation< T >::zz().

180  {
181  if (!sStart.isValid()) {
182  if (sendLogWarning_) {
183  edm::LogWarning("SteppingHelixPropagator") << "Can't propagate: invalid input state" << std::endl;
184  }
185  out = invalidState_;
186  return;
187  }
188  StateArray svBuf;
189  initStateArraySHPSpecific(svBuf, true);
190  int nPoints = 0;
191  setIState(sStart, svBuf, nPoints);
192 
193  Point rPlane(pDest.position().x(), pDest.position().y(), pDest.position().z());
194  Vector nPlane(pDest.rotation().zx(), pDest.rotation().zy(), pDest.rotation().zz());
195  nPlane /= nPlane.mag();
196 
197  double pars[6] = {rPlane.x(), rPlane.y(), rPlane.z(), nPlane.x(), nPlane.y(), nPlane.z()};
198 
199  propagate(svBuf, nPoints, PLANE_DT, pars);
200 
201  //(re)set it before leaving: dir =1 (-1) if path increased (decreased) and 0 if it didn't change
202  //need to implement this somewhere else as a separate function
203  double lDir = 0;
204  if (sStart.path() < svBuf[cIndex_(nPoints - 1)].path())
205  lDir = 1.;
206  if (sStart.path() > svBuf[cIndex_(nPoints - 1)].path())
207  lDir = -1.;
208  svBuf[cIndex_(nPoints - 1)].dir = lDir;
209 
210  out = svBuf[cIndex_(nPoints - 1)];
211  return;
212 }
SteppingHelixPropagator::StateArray
StateInfo StateArray[MAX_POINTS+1]
Definition: SteppingHelixPropagator.h:172
Vector
ROOT::Math::Plane3D::Vector Vector
Definition: EcalHitMaker.cc:29
PV3DBase::y
T y() const
Definition: PV3DBase.h:60
SteppingHelixPropagator::cIndex_
int cIndex_(int ind) const
(Internals) circular index for array of transients
Definition: SteppingHelixPropagator.cc:1773
TkRotation::zx
T zx() const
Definition: extTkRotation.h:199
TkRotation::zz
T zz() const
Definition: extTkRotation.h:201
SteppingHelixPropagator::initStateArraySHPSpecific
void initStateArraySHPSpecific(StateArray &svBuf, bool flagsOnly) const
(Internals) set defaults needed for states used inside propagate methods
Definition: SteppingHelixPropagator.cc:38
PV3DBase::z
T z() const
Definition: PV3DBase.h:61
TkRotation::zy
T zy() const
Definition: extTkRotation.h:200
SteppingHelixPropagator::propagate
void propagate(const SteppingHelixStateInfo &ftsStart, const Surface &sDest, SteppingHelixStateInfo &out) const
Propagate to Plane given a starting point.
Definition: SteppingHelixPropagator.cc:152
Point
Structure Point Contains parameters of Gaussian fits to DMRs.
Definition: DMRtrends.cc:57
SteppingHelixPropagator::setIState
void setIState(const SteppingHelixStateInfo &sStart, StateArray &svBuff, int &nPoints) const
(Internals) Init starting point
Definition: SteppingHelixPropagator.cc:296
GloballyPositioned::rotation
const RotationType & rotation() const
Definition: GloballyPositioned.h:38
edm::LogWarning
Log< level::Warning, false > LogWarning
Definition: MessageLogger.h:122
PV3DBase::x
T x() const
Definition: PV3DBase.h:59
GloballyPositioned::position
const PositionType & position() const
Definition: GloballyPositioned.h:36
void SteppingHelixPropagator::propagate ( const SteppingHelixStateInfo &  ftsStart,
const Cylinder &  cDest,
SteppingHelixStateInfo &  out 
) const

Propagate to Cylinder given a starting point (a Cylinder is assumed to be positioned at 0,0,0)

Definition at line 214 of file SteppingHelixPropagator.cc.

References cIndex_(), initStateArraySHPSpecific(), invalidState_, SteppingHelixStateInfo::isValid(), SteppingHelixStateInfo::path(), propagate(), Cylinder::radius(), RADIUS_DT, RADIUS_P, sendLogWarning_, and setIState().

216  {
217  if (!sStart.isValid()) {
218  if (sendLogWarning_) {
219  edm::LogWarning("SteppingHelixPropagator") << "Can't propagate: invalid input state" << std::endl;
220  }
221  out = invalidState_;
222  return;
223  }
224  StateArray svBuf;
225  initStateArraySHPSpecific(svBuf, true);
226  int nPoints = 0;
227  setIState(sStart, svBuf, nPoints);
228 
229  double pars[6] = {0, 0, 0, 0, 0, 0};
230  pars[RADIUS_P] = cDest.radius();
231 
232  propagate(svBuf, nPoints, RADIUS_DT, pars);
233 
234  //(re)set it before leaving: dir =1 (-1) if path increased (decreased) and 0 if it didn't change
235  //need to implement this somewhere else as a separate function
236  double lDir = 0;
237  if (sStart.path() < svBuf[cIndex_(nPoints - 1)].path())
238  lDir = 1.;
239  if (sStart.path() > svBuf[cIndex_(nPoints - 1)].path())
240  lDir = -1.;
241  svBuf[cIndex_(nPoints - 1)].dir = lDir;
242  out = svBuf[cIndex_(nPoints - 1)];
243  return;
244 }
SteppingHelixPropagator::StateArray
StateInfo StateArray[MAX_POINTS+1]
Definition: SteppingHelixPropagator.h:172
SteppingHelixPropagator::cIndex_
int cIndex_(int ind) const
(Internals) circular index for array of transients
Definition: SteppingHelixPropagator.cc:1773
SteppingHelixPropagator::initStateArraySHPSpecific
void initStateArraySHPSpecific(StateArray &svBuf, bool flagsOnly) const
(Internals) set defaults needed for states used inside propagate methods
Definition: SteppingHelixPropagator.cc:38
Cylinder::radius
Scalar radius() const
Radius of the cylinder.
Definition: Cylinder.h:64
SteppingHelixPropagator::propagate
void propagate(const SteppingHelixStateInfo &ftsStart, const Surface &sDest, SteppingHelixStateInfo &out) const
Propagate to Plane given a starting point.
Definition: SteppingHelixPropagator.cc:152
SteppingHelixPropagator::setIState
void setIState(const SteppingHelixStateInfo &sStart, StateArray &svBuff, int &nPoints) const
(Internals) Init starting point
Definition: SteppingHelixPropagator.cc:296
edm::LogWarning
Log< level::Warning, false > LogWarning
Definition: MessageLogger.h:122
void SteppingHelixPropagator::propagate ( const SteppingHelixStateInfo &  ftsStart,
const GlobalPoint &  pDest,
SteppingHelixStateInfo &  out 
) const

Propagate to PCA to point given a starting point.

Definition at line 246 of file SteppingHelixPropagator.cc.

References cIndex_(), initStateArraySHPSpecific(), invalidState_, SteppingHelixStateInfo::isValid(), POINT_PCA_DT, propagate(), sendLogWarning_, setIState(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

248  {
249  if (!sStart.isValid()) {
250  if (sendLogWarning_) {
251  edm::LogWarning("SteppingHelixPropagator") << "Can't propagate: invalid input state" << std::endl;
252  }
253  out = invalidState_;
254  return;
255  }
256  StateArray svBuf;
257  initStateArraySHPSpecific(svBuf, true);
258  int nPoints = 0;
259  setIState(sStart, svBuf, nPoints);
260 
261  double pars[6] = {pDest.x(), pDest.y(), pDest.z(), 0, 0, 0};
262 
263  propagate(svBuf, nPoints, POINT_PCA_DT, pars);
264 
265  out = svBuf[cIndex_(nPoints - 1)];
266  return;
267 }
SteppingHelixPropagator::StateArray
StateInfo StateArray[MAX_POINTS+1]
Definition: SteppingHelixPropagator.h:172
PV3DBase::y
T y() const
Definition: PV3DBase.h:60
SteppingHelixPropagator::cIndex_
int cIndex_(int ind) const
(Internals) circular index for array of transients
Definition: SteppingHelixPropagator.cc:1773
SteppingHelixPropagator::initStateArraySHPSpecific
void initStateArraySHPSpecific(StateArray &svBuf, bool flagsOnly) const
(Internals) set defaults needed for states used inside propagate methods
Definition: SteppingHelixPropagator.cc:38
PV3DBase::z
T z() const
Definition: PV3DBase.h:61
SteppingHelixPropagator::propagate
void propagate(const SteppingHelixStateInfo &ftsStart, const Surface &sDest, SteppingHelixStateInfo &out) const
Propagate to Plane given a starting point.
Definition: SteppingHelixPropagator.cc:152
SteppingHelixPropagator::setIState
void setIState(const SteppingHelixStateInfo &sStart, StateArray &svBuff, int &nPoints) const
(Internals) Init starting point
Definition: SteppingHelixPropagator.cc:296
edm::LogWarning
Log< level::Warning, false > LogWarning
Definition: MessageLogger.h:122
PV3DBase::x
T x() const
Definition: PV3DBase.h:59
void SteppingHelixPropagator::propagate ( const SteppingHelixStateInfo &  ftsStart,
const GlobalPoint &  pDest1,
const GlobalPoint &  pDest2,
SteppingHelixStateInfo &  out 
) const

Propagate to PCA to a line (given by 2 points) given a starting point.

Definition at line 269 of file SteppingHelixPropagator.cc.

References cIndex_(), alignCSCRings::e, initStateArraySHPSpecific(), invalidState_, SteppingHelixStateInfo::isValid(), LINE_PCA_DT, mag(), propagate(), sendLogWarning_, setIState(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

272  {
273  if ((pDest1 - pDest2).mag() < 1e-10 || !sStart.isValid()) {
274  if (sendLogWarning_) {
275  if ((pDest1 - pDest2).mag() < 1e-10)
276  edm::LogWarning("SteppingHelixPropagator") << "Can't propagate: points are too close" << std::endl;
277  if (!sStart.isValid())
278  edm::LogWarning("SteppingHelixPropagator") << "Can't propagate: invalid input state" << std::endl;
279  }
280  out = invalidState_;
281  return;
282  }
283  StateArray svBuf;
284  initStateArraySHPSpecific(svBuf, true);
285  int nPoints = 0;
286  setIState(sStart, svBuf, nPoints);
287 
288  double pars[6] = {pDest1.x(), pDest1.y(), pDest1.z(), pDest2.x(), pDest2.y(), pDest2.z()};
289 
290  propagate(svBuf, nPoints, LINE_PCA_DT, pars);
291 
292  out = svBuf[cIndex_(nPoints - 1)];
293  return;
294 }
SteppingHelixPropagator::StateArray
StateInfo StateArray[MAX_POINTS+1]
Definition: SteppingHelixPropagator.h:172
mag
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
Definition: Basic3DVectorLD.h:127
PV3DBase::y
T y() const
Definition: PV3DBase.h:60
SteppingHelixPropagator::cIndex_
int cIndex_(int ind) const
(Internals) circular index for array of transients
Definition: SteppingHelixPropagator.cc:1773
SteppingHelixPropagator::initStateArraySHPSpecific
void initStateArraySHPSpecific(StateArray &svBuf, bool flagsOnly) const
(Internals) set defaults needed for states used inside propagate methods
Definition: SteppingHelixPropagator.cc:38
PV3DBase::z
T z() const
Definition: PV3DBase.h:61
SteppingHelixPropagator::propagate
void propagate(const SteppingHelixStateInfo &ftsStart, const Surface &sDest, SteppingHelixStateInfo &out) const
Propagate to Plane given a starting point.
Definition: SteppingHelixPropagator.cc:152
SteppingHelixPropagator::setIState
void setIState(const SteppingHelixStateInfo &sStart, StateArray &svBuff, int &nPoints) const
(Internals) Init starting point
Definition: SteppingHelixPropagator.cc:296
PV3DBase::x
T x() const
Definition: PV3DBase.h:59
SteppingHelixPropagator::Result SteppingHelixPropagator::propagate ( StateArray &  svBuff,
int &  nPoints,
SteppingHelixPropagator::DestType  type,
const double  pars[6],
double  epsilon = 1e-3 
) const
protected

propagate: chose stop point by type argument propagate to fixed radius [ r = sqrt(x**2+y**2) ] with precision epsilon propagate to plane by [x0,y0,z0, n_x, n_y, n_z] parameters

define a fast-skip distance: should be the shortest of a possible step or distance

Definition at line 309 of file SteppingHelixPropagator.cc.

References alongMomentum, anyDirection, cIndex_(), debug_, SteppingHelixStateInfo::dEdx, defaultStep_, DeadROC_duringRun::dir, alignCSCRings::e, SteppingHelixStateInfo::FAULT, SteppingHelixStateInfo::field, field_, HEL_AS_F, SteppingHelixStateInfo::INACC, SteppingHelixStateInfo::isValid_, LINE_PCA_DT, LogTrace, makeAtomStep(), MAX_STEPS, metname, SteppingHelixStateInfo::OK, oppositeToMomentum, SteppingHelixStateInfo::p3, PATHL_DT, PATHL_P, PLANE_DT, POINT_PCA_DT, Propagator::propagationDirection(), SteppingHelixStateInfo::r3, RADIUS_DT, RADIUS_P, SteppingHelixStateInfo::RANGEOUT, refToDest(), refToMagVolume(), refToMatVolume(), mps_fire::result, SteppingHelixStateInfo::ResultName, sendLogWarning_, SteppingHelixStateInfo::status_, AlCaHLTBitMon_QueryRunRegistry::string, SteppingHelixStateInfo::UNDEFINED, useIsYokeFlag_, useMagVolumes_, useMatVolumes_, useTuningForL2Speed_, SteppingHelixStateInfo::WRONG_VOLUME, Z_DT, and Z_P.

313  {
314  static const std::string metname = "SteppingHelixPropagator";
315  StateInfo* svCurrent = &svBuf[cIndex_(nPoints - 1)];
316 
317  //check if it's going to work at all
318  double tanDist = 0;
319  double dist = 0;
320  PropagationDirection refDirection = anyDirection;
321  Result result = refToDest(type, (*svCurrent), pars, dist, tanDist, refDirection);
322 
323  if (result != SteppingHelixStateInfo::OK || fabs(dist) > 1e6) {
324  svCurrent->status_ = result;
325  if (fabs(dist) > 1e6)
326  svCurrent->status_ = SteppingHelixStateInfo::INACC;
327  svCurrent->isValid_ = false;
328  svCurrent->field = field_;
329  if (sendLogWarning_) {
330  edm::LogWarning(metname) << std::setprecision(17) << std::setw(20) << std::scientific
331  << " Failed after first refToDest check with status "
332  << SteppingHelixStateInfo::ResultName[result] << std::endl;
333  }
334  return result;
335  }
336 
337  result = SteppingHelixStateInfo::UNDEFINED;
338  bool makeNextStep = true;
339  double dStep = defaultStep_;
340  PropagationDirection dir, oldDir;
341  dir = propagationDirection();
342  oldDir = dir;
343  int nOsc = 0;
344 
345  double distMag = 1e12;
346  double tanDistMag = 1e12;
347 
348  double distMat = 1e12;
349  double tanDistMat = 1e12;
350 
351  double tanDistNextCheck = -0.1; //just need a negative start val
352  double tanDistMagNextCheck = -0.1;
353  double tanDistMatNextCheck = -0.1;
354  double oldDStep = 0;
355  PropagationDirection oldRefDirection = propagationDirection();
356 
357  Result resultToMat = SteppingHelixStateInfo::UNDEFINED;
358  Result resultToMag = SteppingHelixStateInfo::UNDEFINED;
359 
360  bool isFirstStep = true;
361  bool expectNewMagVolume = false;
362 
363  int loopCount = 0;
364  while (makeNextStep) {
365  dStep = defaultStep_;
366  svCurrent = &svBuf[cIndex_(nPoints - 1)];
367  double curZ = svCurrent->r3.z();
368  double curR = svCurrent->r3.perp();
369  if (fabs(curZ) < 440 && curR < 260)
370  dStep = defaultStep_ * 2;
371 
372  //more such ifs might be scattered around
373  //even though dStep is large, it will still make stops at each volume boundary
374  if (useTuningForL2Speed_) {
375  dStep = 100.;
376  if (!useIsYokeFlag_ && fabs(curZ) < 667 && curR > 380 && curR < 850) {
377  dStep = 5 * (1 + 0.2 * svCurrent->p3.mag());
378  }
379  }
380 
381  // refDirection = propagationDirection();
382 
383  tanDistNextCheck -= oldDStep;
384  tanDistMagNextCheck -= oldDStep;
385  tanDistMatNextCheck -= oldDStep;
386 
387  if (tanDistNextCheck < 0) {
388  //use pre-computed values if it's the first step
389  if (!isFirstStep)
390  refToDest(type, (*svCurrent), pars, dist, tanDist, refDirection);
391  // constrain allowed path for a tangential approach
392  if (fabs(tanDist) > 4. * (fabs(dist)))
393  tanDist *= tanDist == 0 ? 0 : fabs(dist / tanDist * 4.);
394 
395  tanDistNextCheck = fabs(tanDist) * 0.5 - 0.5; //need a better guess (to-do)
396  //reasonable limit
397  if (tanDistNextCheck > defaultStep_ * 20.)
398  tanDistNextCheck = defaultStep_ * 20.;
399  oldRefDirection = refDirection;
400  } else {
401  tanDist = tanDist > 0. ? tanDist - oldDStep : tanDist + oldDStep;
402  refDirection = oldRefDirection;
403  if (debug_)
404  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
405  << "Skipped refToDest: guess tanDist = " << tanDist << " next check at " << tanDistNextCheck
406  << std::endl;
407  }
409  double fastSkipDist = fabs(dist) > fabs(tanDist) ? tanDist : dist;
410 
411  if (propagationDirection() == anyDirection) {
412  dir = refDirection;
413  } else {
414  dir = propagationDirection();
415  if (fabs(tanDist) < 0.1 && refDirection != dir) {
416  //how did it get here? nOsc++;
417  dir = refDirection;
418  if (debug_)
419  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
420  << "NOTE: overstepped last time: switch direction (can do it if within 1 mm)" << std::endl;
421  }
422  }
423 
424  if (useMagVolumes_ && !(fabs(dist) < fabs(epsilon))) { //need to know the general direction
425  if (tanDistMagNextCheck < 0) {
426  resultToMag = refToMagVolume(
427  (*svCurrent), dir, distMag, tanDistMag, fabs(fastSkipDist), expectNewMagVolume, fabs(tanDist));
428  // constrain allowed path for a tangential approach
429  if (fabs(tanDistMag) > 4. * (fabs(distMag)))
430  tanDistMag *= tanDistMag == 0 ? 0 : fabs(distMag / tanDistMag * 4.);
431 
432  tanDistMagNextCheck = fabs(tanDistMag) * 0.5 - 0.5; //need a better guess (to-do)
433  //reasonable limit; "turn off" checking if bounds are further than the destination
434  if (tanDistMagNextCheck > defaultStep_ * 20. || fabs(dist) < fabs(distMag) ||
435  resultToMag == SteppingHelixStateInfo::INACC)
436  tanDistMagNextCheck = defaultStep_ * 20 > fabs(fastSkipDist) ? fabs(fastSkipDist) : defaultStep_ * 20;
437  if (resultToMag != SteppingHelixStateInfo::INACC && resultToMag != SteppingHelixStateInfo::OK)
438  tanDistMagNextCheck = -1;
439  } else {
440  // resultToMag = SteppingHelixStateInfo::OK;
441  tanDistMag = tanDistMag > 0. ? tanDistMag - oldDStep : tanDistMag + oldDStep;
442  if (debug_)
443  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
444  << "Skipped refToMag: guess tanDistMag = " << tanDistMag << " next check at "
445  << tanDistMagNextCheck;
446  }
447  }
448 
449  if (useMatVolumes_ && !(fabs(dist) < fabs(epsilon))) { //need to know the general direction
450  if (tanDistMatNextCheck < 0) {
451  resultToMat = refToMatVolume((*svCurrent), dir, distMat, tanDistMat, fabs(fastSkipDist));
452  // constrain allowed path for a tangential approach
453  if (fabs(tanDistMat) > 4. * (fabs(distMat)))
454  tanDistMat *= tanDistMat == 0 ? 0 : fabs(distMat / tanDistMat * 4.);
455 
456  tanDistMatNextCheck = fabs(tanDistMat) * 0.5 - 0.5; //need a better guess (to-do)
457  //reasonable limit; "turn off" checking if bounds are further than the destination
458  if (tanDistMatNextCheck > defaultStep_ * 20. || fabs(dist) < fabs(distMat) ||
459  resultToMat == SteppingHelixStateInfo::INACC)
460  tanDistMatNextCheck = defaultStep_ * 20 > fabs(fastSkipDist) ? fabs(fastSkipDist) : defaultStep_ * 20;
461  if (resultToMat != SteppingHelixStateInfo::INACC && resultToMat != SteppingHelixStateInfo::OK)
462  tanDistMatNextCheck = -1;
463  } else {
464  // resultToMat = SteppingHelixStateInfo::OK;
465  tanDistMat = tanDistMat > 0. ? tanDistMat - oldDStep : tanDistMat + oldDStep;
466  if (debug_)
467  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
468  << "Skipped refToMat: guess tanDistMat = " << tanDistMat << " next check at "
469  << tanDistMatNextCheck;
470  }
471  }
472 
473  double rDotP = svCurrent->r3.dot(svCurrent->p3);
474  if ((fabs(curZ) > 1.5e3 || curR > 800.) &&
475  ((dir == alongMomentum && rDotP > 0) || (dir == oppositeToMomentum && rDotP < 0))) {
476  dStep = fabs(tanDist) - 1e-12;
477  }
478  double tanDistMin = fabs(tanDist);
479  if (tanDistMin > fabs(tanDistMag) + 0.05 &&
480  (resultToMag == SteppingHelixStateInfo::OK || resultToMag == SteppingHelixStateInfo::WRONG_VOLUME)) {
481  tanDistMin = fabs(tanDistMag) + 0.05; //try to step into the next volume
482  expectNewMagVolume = true;
483  } else
484  expectNewMagVolume = false;
485 
486  if (tanDistMin > fabs(tanDistMat) + 0.05 && resultToMat == SteppingHelixStateInfo::OK) {
487  tanDistMin = fabs(tanDistMat) + 0.05; //try to step into the next volume
488  if (expectNewMagVolume)
489  expectNewMagVolume = false;
490  }
491 
492  if (tanDistMin * fabs(svCurrent->dEdx) > 0.5 * svCurrent->p3.mag()) {
493  tanDistMin = 0.5 * svCurrent->p3.mag() / fabs(svCurrent->dEdx);
494  if (expectNewMagVolume)
495  expectNewMagVolume = false;
496  }
497 
498  double tanDistMinLazy = fabs(tanDistMin);
499  if ((type == POINT_PCA_DT || type == LINE_PCA_DT) && fabs(tanDist) < 2. * fabs(tanDistMin)) {
500  //being lazy here; the best is to take into account the curvature
501  tanDistMinLazy = fabs(tanDistMin) * 0.5;
502  }
503 
504  if (fabs(tanDistMinLazy) < dStep) {
505  dStep = fabs(tanDistMinLazy);
506  }
507 
508  //keep this path length for the next step
509  oldDStep = dStep;
510 
511  if (dStep > 1e-10 && !(fabs(dist) < fabs(epsilon))) {
512  StateInfo* svNext = &svBuf[cIndex_(nPoints)];
513  makeAtomStep((*svCurrent), (*svNext), dStep, dir, HEL_AS_F);
514  // if (useMatVolumes_ && expectNewMagVolume
515  // && svCurrent->magVol == svNext->magVol){
516  // double tmpDist=0;
517  // double tmpDistMag = 0;
518  // if (refToMagVolume((*svNext), dir, tmpDist, tmpDistMag, fabs(dist)) != SteppingHelixStateInfo::OK){
519  // //the point appears to be outside, but findVolume claims the opposite
520  // dStep += 0.05*fabs(tanDistMag/distMag); oldDStep = dStep; //do it again with a bigger step
521  // if (debug_) LogTrace(metname)
522  // <<"Failed to get into new mag volume: will try with new bigger step "<<dStep<<std::endl;
523  // makeAtomStep((*svCurrent), (*svNext), dStep, dir, HEL_AS_F);
524  // }
525  // }
526  nPoints++;
527  svCurrent = &svBuf[cIndex_(nPoints - 1)];
528  if (oldDir != dir) {
529  nOsc++;
530  tanDistNextCheck = -1; //check dist after osc
531  tanDistMagNextCheck = -1;
532  tanDistMatNextCheck = -1;
533  }
534  oldDir = dir;
535  }
536 
537  if (nOsc > 1 && fabs(dStep) > epsilon) {
538  if (debug_)
539  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Ooops" << std::endl;
540  }
541 
542  if (fabs(dist) < fabs(epsilon))
543  result = SteppingHelixStateInfo::OK;
544 
545  if ((type == POINT_PCA_DT || type == LINE_PCA_DT) && fabs(dStep) < fabs(epsilon)) {
546  //now check if it's not a branch point (peek ahead at 1 cm)
547  double nextDist = 0;
548  double nextTanDist = 0;
549  PropagationDirection nextRefDirection = anyDirection;
550  StateInfo* svNext = &svBuf[cIndex_(nPoints)];
551  makeAtomStep((*svCurrent), (*svNext), 1., dir, HEL_AS_F);
552  nPoints++;
553  svCurrent = &svBuf[cIndex_(nPoints - 1)];
554  refToDest(type, (*svCurrent), pars, nextDist, nextTanDist, nextRefDirection);
555  if (fabs(nextDist) > fabs(dist)) {
556  nPoints--;
557  svCurrent = &svBuf[cIndex_(nPoints - 1)];
558  result = SteppingHelixStateInfo::OK;
559  if (debug_) {
560  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
561  << "Found real local minimum in PCA" << std::endl;
562  }
563  } else {
564  //keep this trial point and continue
565  dStep = defaultStep_;
566  if (debug_) {
567  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Found branch point in PCA"
568  << std::endl;
569  }
570  }
571  }
572 
573  if (nPoints > MAX_STEPS * 1. / defaultStep_ || loopCount > MAX_STEPS * 100 || nOsc > 6)
574  result = SteppingHelixStateInfo::FAULT;
575 
576  if (svCurrent->p3.mag() < 0.1)
577  result = SteppingHelixStateInfo::RANGEOUT;
578 
579  curZ = svCurrent->r3.z();
580  curR = svCurrent->r3.perp();
581  if (curR > 20000 || fabs(curZ) > 20000)
582  result = SteppingHelixStateInfo::INACC;
583 
584  makeNextStep = result == SteppingHelixStateInfo::UNDEFINED;
585  svCurrent->status_ = result;
586  svCurrent->isValid_ = (result == SteppingHelixStateInfo::OK || makeNextStep);
587  svCurrent->field = field_;
588 
589  isFirstStep = false;
590  loopCount++;
591  }
592 
593  if (sendLogWarning_ && result != SteppingHelixStateInfo::OK) {
594  edm::LogWarning(metname) << std::setprecision(17) << std::setw(20) << std::scientific
595  << " Propagation failed with status " << SteppingHelixStateInfo::ResultName[result]
596  << std::endl;
597  if (result == SteppingHelixStateInfo::RANGEOUT)
598  edm::LogWarning(metname) << std::setprecision(17) << std::setw(20) << std::scientific
599  << " Momentum at last point is too low (<0.1) p_last = " << svCurrent->p3.mag()
600  << std::endl;
601  if (result == SteppingHelixStateInfo::INACC)
602  edm::LogWarning(metname) << std::setprecision(17) << std::setw(20) << std::scientific
603  << " Went too far: the last point is at " << svCurrent->r3 << std::endl;
604  if (result == SteppingHelixStateInfo::FAULT && nOsc > 6)
605  edm::LogWarning(metname) << std::setprecision(17) << std::setw(20) << std::scientific
606  << " Infinite loop condidtion detected: going in cycles. Break after 6 cycles"
607  << std::endl;
608  if (result == SteppingHelixStateInfo::FAULT && nPoints > MAX_STEPS * 1. / defaultStep_)
609  edm::LogWarning(metname) << std::setprecision(17) << std::setw(20) << std::scientific
610  << " Tired to go farther. Made too many steps: more than "
611  << MAX_STEPS * 1. / defaultStep_ << std::endl;
612  }
613 
614  if (debug_) {
615  switch (type) {
616  case RADIUS_DT:
617  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "going to radius "
618  << pars[RADIUS_P] << std::endl;
619  break;
620  case Z_DT:
621  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "going to z " << pars[Z_P]
622  << std::endl;
623  break;
624  case PATHL_DT:
625  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "going to pathL "
626  << pars[PATHL_P] << std::endl;
627  break;
628  case PLANE_DT: {
629  Point rPlane(pars[0], pars[1], pars[2]);
630  Vector nPlane(pars[3], pars[4], pars[5]);
631  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "going to plane r0:" << rPlane
632  << " n:" << nPlane << std::endl;
633  } break;
634  case POINT_PCA_DT: {
635  Point rDest(pars[0], pars[1], pars[2]);
636  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "going to PCA to point "
637  << rDest << std::endl;
638  } break;
639  case LINE_PCA_DT: {
640  Point rDest1(pars[0], pars[1], pars[2]);
641  Point rDest2(pars[3], pars[4], pars[5]);
642  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "going to PCA to line "
643  << rDest1 << " - " << rDest2 << std::endl;
644  } break;
645  default:
646  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "going to NOT IMPLEMENTED"
647  << std::endl;
648  break;
649  }
650  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Made " << nPoints - 1
651  << " steps and stopped at(cur step) " << svCurrent->r3 << " nOsc " << nOsc << std::endl;
652  }
653 
654  return result;
655 }
SteppingHelixPropagator::refToDest
Result refToDest(DestType dest, const SteppingHelixPropagator::StateInfo &sv, const double pars[6], double &dist, double &tanDist, PropagationDirection &refDirection, double fastSkipDist=1e12) const
(Internals) determine distance and direction from the current position to the plane ...
Definition: SteppingHelixPropagator.cc:1781
metname
const std::string metname
Definition: MuonSeedOrcaPatternRecognition.cc:40
Vector
ROOT::Math::Plane3D::Vector Vector
Definition: EcalHitMaker.cc:29
SteppingHelixPropagator::StateInfo
SteppingHelixStateInfo StateInfo
Definition: SteppingHelixPropagator.h:41
Propagator::propagationDirection
virtual PropagationDirection propagationDirection() const final
Definition: Propagator.h:139
SteppingHelixPropagator::cIndex_
int cIndex_(int ind) const
(Internals) circular index for array of transients
Definition: SteppingHelixPropagator.cc:1773
PropagationDirection
PropagationDirection
Definition: PropagationDirection.h:4
SteppingHelixPropagator::refToMatVolume
Result refToMatVolume(const SteppingHelixPropagator::StateInfo &sv, PropagationDirection dir, double &dist, double &tanDist, double fastSkipDist=1e12) const
Definition: SteppingHelixPropagator.cc:2308
type
type
Definition: SiPixelVCal_PayloadInspector.cc:39
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
LogTrace
#define LogTrace(id)
Definition: MessageLogger.h:234
mps_fire.result
tuple result
Definition: mps_fire.py:311
SteppingHelixPropagator::Result
SteppingHelixStateInfo::Result Result
Definition: SteppingHelixPropagator.h:42
DeadROC_duringRun.dir
string dir
Definition: DeadROC_duringRun.py:23
SteppingHelixPropagator::makeAtomStep
bool makeAtomStep(SteppingHelixPropagator::StateInfo &svCurrent, SteppingHelixPropagator::StateInfo &svNext, double dS, PropagationDirection dir, SteppingHelixPropagator::Fancy fancy) const
main stepping function: compute next state vector after a step of length dS
Definition: SteppingHelixPropagator.cc:879
SteppingHelixPropagator::refToMagVolume
Result refToMagVolume(const SteppingHelixPropagator::StateInfo &sv, PropagationDirection dir, double &dist, double &tanDist, double fastSkipDist=1e12, bool expectNewMagVolume=false, double maxStep=1e12) const
Definition: SteppingHelixPropagator.cc:2049
Point
Structure Point Contains parameters of Gaussian fits to DMRs.
Definition: DMRtrends.cc:57
edm::LogWarning
Log< level::Warning, false > LogWarning
Definition: MessageLogger.h:122
SteppingHelixPropagator::field_
const MagneticField * field_
Definition: SteppingHelixPropagator.h:265
SteppingHelixStateInfo::ResultName
static const std::string ResultName[MAX_RESULT]
Definition: SteppingHelixStateInfo.h:36
std::pair< TrajectoryStateOnSurface, double > SteppingHelixPropagator::propagateWithPath ( const FreeTrajectoryState &  ftsStart,
const Plane &  pDest 
) const
overridevirtual

Propagate to Plane given a starting point: return final TrajectoryState and path length from start to this point

Implements Propagator.

Definition at line 79 of file SteppingHelixPropagator.cc.

References SteppingHelixStateInfo::getStateOnSurface(), initStateArraySHPSpecific(), SteppingHelixStateInfo::path(), propagate(), and setIState().

Referenced by PropagateToCal::propagate(), and propagateWithPath().

80  {
81  StateArray svBuf;
82  initStateArraySHPSpecific(svBuf, true);
83  int nPoints = 0;
84  setIState(SteppingHelixStateInfo(ftsStart), svBuf, nPoints);
85 
86  StateInfo svCurrent;
87  propagate(svBuf[0], pDest, svCurrent);
88 
89  return TsosPP(svCurrent.getStateOnSurface(pDest), svCurrent.path());
90 }
SteppingHelixPropagator::StateArray
StateInfo StateArray[MAX_POINTS+1]
Definition: SteppingHelixPropagator.h:172
SteppingHelixPropagator::StateInfo
SteppingHelixStateInfo StateInfo
Definition: SteppingHelixPropagator.h:41
SteppingHelixPropagator::TsosPP
std::pair< TrajectoryStateOnSurface, double > TsosPP
Definition: SteppingHelixPropagator.h:257
SteppingHelixPropagator::initStateArraySHPSpecific
void initStateArraySHPSpecific(StateArray &svBuf, bool flagsOnly) const
(Internals) set defaults needed for states used inside propagate methods
Definition: SteppingHelixPropagator.cc:38
SteppingHelixPropagator::propagate
void propagate(const SteppingHelixStateInfo &ftsStart, const Surface &sDest, SteppingHelixStateInfo &out) const
Propagate to Plane given a starting point.
Definition: SteppingHelixPropagator.cc:152
SteppingHelixPropagator::setIState
void setIState(const SteppingHelixStateInfo &sStart, StateArray &svBuff, int &nPoints) const
(Internals) Init starting point
Definition: SteppingHelixPropagator.cc:296
std::pair< TrajectoryStateOnSurface, double > SteppingHelixPropagator::propagateWithPath ( const FreeTrajectoryState &  ftsStart,
const Cylinder &  cDest 
) const
overridevirtual

Propagate to Cylinder given a starting point: return final TrajectoryState and path length from start to this point

Implements Propagator.

Definition at line 92 of file SteppingHelixPropagator.cc.

References SteppingHelixStateInfo::getStateOnSurface(), initStateArraySHPSpecific(), SteppingHelixStateInfo::path(), propagate(), returnTangentPlane_, and setIState().

93  {
94  StateArray svBuf;
95  initStateArraySHPSpecific(svBuf, true);
96  int nPoints = 0;
97  setIState(SteppingHelixStateInfo(ftsStart), svBuf, nPoints);
98 
99  StateInfo svCurrent;
100  propagate(svBuf[0], cDest, svCurrent);
101 
102  return TsosPP(svCurrent.getStateOnSurface(cDest, returnTangentPlane_), svCurrent.path());
103 }
SteppingHelixPropagator::StateArray
StateInfo StateArray[MAX_POINTS+1]
Definition: SteppingHelixPropagator.h:172
SteppingHelixPropagator::StateInfo
SteppingHelixStateInfo StateInfo
Definition: SteppingHelixPropagator.h:41
SteppingHelixPropagator::TsosPP
std::pair< TrajectoryStateOnSurface, double > TsosPP
Definition: SteppingHelixPropagator.h:257
SteppingHelixPropagator::initStateArraySHPSpecific
void initStateArraySHPSpecific(StateArray &svBuf, bool flagsOnly) const
(Internals) set defaults needed for states used inside propagate methods
Definition: SteppingHelixPropagator.cc:38
SteppingHelixPropagator::propagate
void propagate(const SteppingHelixStateInfo &ftsStart, const Surface &sDest, SteppingHelixStateInfo &out) const
Propagate to Plane given a starting point.
Definition: SteppingHelixPropagator.cc:152
SteppingHelixPropagator::setIState
void setIState(const SteppingHelixStateInfo &sStart, StateArray &svBuff, int &nPoints) const
(Internals) Init starting point
Definition: SteppingHelixPropagator.cc:296
std::pair< FreeTrajectoryState, double > SteppingHelixPropagator::propagateWithPath ( const FreeTrajectoryState &  ftsStart,
const GlobalPoint &  pDest 
) const
overridevirtual

Propagate to PCA to point given a starting point.

Reimplemented from Propagator.

Definition at line 105 of file SteppingHelixPropagator.cc.

References SteppingHelixStateInfo::getFreeState(), initStateArraySHPSpecific(), SteppingHelixStateInfo::path(), propagate(), and setIState().

106  {
107  StateArray svBuf;
108  initStateArraySHPSpecific(svBuf, true);
109  int nPoints = 0;
110  setIState(SteppingHelixStateInfo(ftsStart), svBuf, nPoints);
111 
112  StateInfo svCurrent;
113  propagate(svBuf[0], pDest, svCurrent);
114 
115  FreeTrajectoryState ftsDest;
116  svCurrent.getFreeState(ftsDest);
117 
118  return FtsPP(ftsDest, svCurrent.path());
119 }
SteppingHelixPropagator::StateArray
StateInfo StateArray[MAX_POINTS+1]
Definition: SteppingHelixPropagator.h:172
SteppingHelixPropagator::StateInfo
SteppingHelixStateInfo StateInfo
Definition: SteppingHelixPropagator.h:41
SteppingHelixPropagator::FtsPP
std::pair< FreeTrajectoryState, double > FtsPP
Definition: SteppingHelixPropagator.h:258
SteppingHelixPropagator::initStateArraySHPSpecific
void initStateArraySHPSpecific(StateArray &svBuf, bool flagsOnly) const
(Internals) set defaults needed for states used inside propagate methods
Definition: SteppingHelixPropagator.cc:38
SteppingHelixPropagator::propagate
void propagate(const SteppingHelixStateInfo &ftsStart, const Surface &sDest, SteppingHelixStateInfo &out) const
Propagate to Plane given a starting point.
Definition: SteppingHelixPropagator.cc:152
SteppingHelixPropagator::setIState
void setIState(const SteppingHelixStateInfo &sStart, StateArray &svBuff, int &nPoints) const
(Internals) Init starting point
Definition: SteppingHelixPropagator.cc:296
std::pair< FreeTrajectoryState, double > SteppingHelixPropagator::propagateWithPath ( const FreeTrajectoryState &  ftsStart,
const GlobalPoint &  pDest1,
const GlobalPoint &  pDest2 
) const
overridevirtual

Propagate to PCA to a line (given by 2 points) given a starting point.

Reimplemented from Propagator.

Definition at line 121 of file SteppingHelixPropagator.cc.

References alignCSCRings::e, SteppingHelixStateInfo::getFreeState(), initStateArraySHPSpecific(), mag(), SteppingHelixStateInfo::path(), propagate(), sendLogWarning_, and setIState().

123  {
124  if ((pDest1 - pDest2).mag() < 1e-10) {
125  if (sendLogWarning_) {
126  edm::LogWarning("SteppingHelixPropagator")
127  << "Can't propagate: the points should be at a bigger distance" << std::endl;
128  }
129  return FtsPP();
130  }
131  StateArray svBuf;
132  initStateArraySHPSpecific(svBuf, true);
133  int nPoints = 0;
134  setIState(SteppingHelixStateInfo(ftsStart), svBuf, nPoints);
135 
136  StateInfo svCurrent;
137  propagate(svBuf[0], pDest1, pDest2, svCurrent);
138 
139  FreeTrajectoryState ftsDest;
140  svCurrent.getFreeState(ftsDest);
141 
142  return FtsPP(ftsDest, svCurrent.path());
143 }
SteppingHelixPropagator::StateArray
StateInfo StateArray[MAX_POINTS+1]
Definition: SteppingHelixPropagator.h:172
mag
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
Definition: Basic3DVectorLD.h:127
SteppingHelixPropagator::StateInfo
SteppingHelixStateInfo StateInfo
Definition: SteppingHelixPropagator.h:41
SteppingHelixPropagator::FtsPP
std::pair< FreeTrajectoryState, double > FtsPP
Definition: SteppingHelixPropagator.h:258
SteppingHelixPropagator::initStateArraySHPSpecific
void initStateArraySHPSpecific(StateArray &svBuf, bool flagsOnly) const
(Internals) set defaults needed for states used inside propagate methods
Definition: SteppingHelixPropagator.cc:38
SteppingHelixPropagator::propagate
void propagate(const SteppingHelixStateInfo &ftsStart, const Surface &sDest, SteppingHelixStateInfo &out) const
Propagate to Plane given a starting point.
Definition: SteppingHelixPropagator.cc:152
SteppingHelixPropagator::setIState
void setIState(const SteppingHelixStateInfo &sStart, StateArray &svBuff, int &nPoints) const
(Internals) Init starting point
Definition: SteppingHelixPropagator.cc:296
edm::LogWarning
Log< level::Warning, false > LogWarning
Definition: MessageLogger.h:122
std::pair< FreeTrajectoryState, double > SteppingHelixPropagator::propagateWithPath ( const FreeTrajectoryState &  ftsStart,
const reco::BeamSpot &  beamSpot 
) const
overridevirtual

Propagate to PCA to a line (given by beamSpot position and slope) given a starting point.

Reimplemented from Propagator.

Definition at line 145 of file SteppingHelixPropagator.cc.

References reco::BeamSpot::dxdz(), reco::BeamSpot::dydz(), propagateWithPath(), reco::BeamSpot::x0(), reco::BeamSpot::y0(), and reco::BeamSpot::z0().

146  {
147  GlobalPoint pDest1(beamSpot.x0(), beamSpot.y0(), beamSpot.z0());
148  GlobalPoint pDest2(pDest1.x() + beamSpot.dxdz() * 1e3, pDest1.y() + beamSpot.dydz() * 1e3, pDest1.z() + 1e3);
149  return propagateWithPath(ftsStart, pDest1, pDest2);
150 }
reco::BeamSpot::z0
double z0() const
z coordinate
Definition: BeamSpot.h:65
reco::BeamSpot::dydz
double dydz() const
dydz slope
Definition: BeamSpot.h:80
reco::BeamSpot::dxdz
double dxdz() const
dxdz slope
Definition: BeamSpot.h:78
reco::BeamSpot::y0
double y0() const
y coordinate
Definition: BeamSpot.h:63
SteppingHelixPropagator::propagateWithPath
std::pair< TrajectoryStateOnSurface, double > propagateWithPath(const FreeTrajectoryState &ftsStart, const Plane &pDest) const override
Definition: SteppingHelixPropagator.cc:79
reco::BeamSpot::x0
double x0() const
x coordinate
Definition: BeamSpot.h:61
SteppingHelixPropagator::Result SteppingHelixPropagator::refToDest ( SteppingHelixPropagator::DestType  dest,
const SteppingHelixPropagator::StateInfo &  sv,
const double  pars[6],
double &  dist,
double &  tanDist,
PropagationDirection &  refDirection,
double  fastSkipDist = 1e12 
) const
protected

(Internals) determine distance and direction from the current position to the plane

Definition at line 1781 of file SteppingHelixPropagator.cc.

References alongMomentum, anyDirection, b0, SteppingHelixStateInfo::bf, CONE_DT, funct::cos(), debug_, alignCSCRings::e, mps_fire::i, SteppingHelixStateInfo::INACC, LINE_PCA_DT, LogTrace, mag(), metname, SteppingHelixStateInfo::NOT_IMPLEMENTED, SteppingHelixStateInfo::OK, oppositeToMomentum, SteppingHelixStateInfo::p3, SteppingHelixStateInfo::path(), PATHL_DT, PATHL_P, Geom::pi(), PLANE_DT, POINT_PCA_DT, SteppingHelixStateInfo::q, SteppingHelixStateInfo::r3, RADIUS_DT, RADIUS_P, mps_fire::result, funct::sin(), mathSSE::sqrt(), AlCaHLTBitMon_QueryRunRegistry::string, metsig::tau, theta(), Z_DT, and Z_P.

Referenced by propagate(), refToMagVolume(), and refToMatVolume().

1787  {
1788  static const std::string metname = "SteppingHelixPropagator";
1789  Result result = SteppingHelixStateInfo::NOT_IMPLEMENTED;
1790 
1791  switch (dest) {
1792  case RADIUS_DT: {
1793  double curR = sv.r3.perp();
1794  dist = pars[RADIUS_P] - curR;
1795  if (fabs(dist) > fastSkipDist) {
1796  result = SteppingHelixStateInfo::INACC;
1797  break;
1798  }
1799  double curP2 = sv.p3.mag2();
1800  double curPtPos2 = sv.p3.perp2();
1801  if (curPtPos2 < 1e-16)
1802  curPtPos2 = 1e-16;
1803 
1804  double cosDPhiPR =
1805  (sv.r3.x() * sv.p3.x() + sv.r3.y() * sv.p3.y()); //only the sign is needed cos((sv.r3.deltaPhi(sv.p3)));
1806  refDirection = dist * cosDPhiPR > 0 ? alongMomentum : oppositeToMomentum;
1807  tanDist = dist * sqrt(curP2 / curPtPos2);
1808  result = SteppingHelixStateInfo::OK;
1809  } break;
1810  case Z_DT: {
1811  double curZ = sv.r3.z();
1812  dist = pars[Z_P] - curZ;
1813  if (fabs(dist) > fastSkipDist) {
1814  result = SteppingHelixStateInfo::INACC;
1815  break;
1816  }
1817  double curP = sv.p3.mag();
1818  refDirection = sv.p3.z() * dist > 0. ? alongMomentum : oppositeToMomentum;
1819  tanDist = dist / sv.p3.z() * curP;
1820  result = SteppingHelixStateInfo::OK;
1821  } break;
1822  case PLANE_DT: {
1823  Point rPlane(pars[0], pars[1], pars[2]);
1824  Vector nPlane(pars[3], pars[4], pars[5]);
1825 
1826  // unfortunately this doesn't work: the numbers are too large
1827  // bool pVertical = fabs(pars[5])>0.9999;
1828  // double dRDotN = pVertical? (sv.r3.z() - rPlane.z())*nPlane.z() :(sv.r3 - rPlane).dot(nPlane);
1829  double dRDotN = (sv.r3.x() - rPlane.x()) * nPlane.x() + (sv.r3.y() - rPlane.y()) * nPlane.y() +
1830  (sv.r3.z() - rPlane.z()) * nPlane.z(); //(sv.r3 - rPlane).dot(nPlane);
1831 
1832  dist = fabs(dRDotN);
1833  if (dist > fastSkipDist) {
1834  result = SteppingHelixStateInfo::INACC;
1835  break;
1836  }
1837  double curP = sv.p3.mag();
1838  double p0 = curP;
1839  double p0Inv = 1. / p0;
1840  Vector tau(sv.p3);
1841  tau *= p0Inv;
1842  double tN = tau.dot(nPlane);
1843  refDirection = tN * dRDotN < 0. ? alongMomentum : oppositeToMomentum;
1844  double b0 = sv.bf.mag();
1845  if (fabs(tN) > 1e-24) {
1846  tanDist = -dRDotN / tN;
1847  } else {
1848  tN = 1e-24;
1849  if (fabs(dRDotN) > 1e-24)
1850  tanDist = 1e6;
1851  else
1852  tanDist = 1;
1853  }
1854  if (fabs(tanDist) > 1e4)
1855  tanDist = 1e4;
1856  if (b0 > 1.5e-6) {
1857  double b0Inv = 1. / b0;
1858  double tNInv = 1. / tN;
1859  Vector bHat(sv.bf);
1860  bHat *= b0Inv;
1861  double bHatN = bHat.dot(nPlane);
1862  double cosPB = bHat.dot(tau);
1863  double kVal = 2.99792458e-3 * sv.q * p0Inv * b0;
1864  double aVal = tanDist * kVal;
1865  Vector lVec = bHat.cross(tau);
1866  double bVal = lVec.dot(nPlane) * tNInv;
1867  if (fabs(aVal * bVal) < 0.3) {
1868  double cVal =
1869  1. - bHatN * cosPB * tNInv; // - sv.bf.cross(lVec).dot(nPlane)*b0Inv*tNInv; //1- bHat_n*bHat_tau/tau_n;
1870  double aacVal = cVal * aVal * aVal;
1871  if (fabs(aacVal) < 1) {
1872  double tanDCorr = bVal / 2. + (bVal * bVal / 2. + cVal / 6) * aVal;
1873  tanDCorr *= aVal;
1874  //+ (-bVal/24. + 0.625*bVal*bVal*bVal + 5./12.*bVal*cVal)*aVal*aVal*aVal
1875  if (debug_)
1876  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << tanDist << " vs "
1877  << tanDist * (1. + tanDCorr) << " corr " << tanDist * tanDCorr << std::endl;
1878  tanDist *= (1. + tanDCorr);
1879  } else {
1880  if (debug_)
1881  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "AACVal "
1882  << fabs(aacVal) << " = " << aVal << "**2 * " << cVal << " too large:: will not converge"
1883  << std::endl;
1884  }
1885  } else {
1886  if (debug_)
1887  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "ABVal "
1888  << fabs(aVal * bVal) << " = " << aVal << " * " << bVal << " too large:: will not converge"
1889  << std::endl;
1890  }
1891  }
1892  result = SteppingHelixStateInfo::OK;
1893  } break;
1894  case CONE_DT: {
1895  Point cVertex(pars[0], pars[1], pars[2]);
1896  if (cVertex.perp2() < 1e-10) {
1897  //assumes the cone axis/vertex is along z
1898  Vector relV3 = sv.r3 - cVertex;
1899  double relV3mag = relV3.mag();
1900  // double relV3Theta = relV3.theta();
1901  double theta(pars[3]);
1902  // double dTheta = theta-relV3Theta;
1903  double sinTheta = sin(theta);
1904  double cosTheta = cos(theta);
1905  double cosV3Theta = relV3.z() / relV3mag;
1906  if (cosV3Theta > 1)
1907  cosV3Theta = 1;
1908  if (cosV3Theta < -1)
1909  cosV3Theta = -1;
1910  double sinV3Theta = sqrt(1. - cosV3Theta * cosV3Theta);
1911 
1912  double sinDTheta = sinTheta * cosV3Theta - cosTheta * sinV3Theta; //sin(dTheta);
1913  double cosDTheta = cosTheta * cosV3Theta + sinTheta * sinV3Theta; //cos(dTheta);
1914  bool isInside = sinTheta > sinV3Theta && cosTheta * cosV3Theta > 0;
1915  dist = isInside || cosDTheta > 0 ? relV3mag * sinDTheta : relV3mag;
1916  if (fabs(dist) > fastSkipDist) {
1917  result = SteppingHelixStateInfo::INACC;
1918  break;
1919  }
1920 
1921  double relV3phi = relV3.phi();
1922  double normPhi = isInside ? Geom::pi() + relV3phi : relV3phi;
1923  double normTheta = theta > Geom::pi() / 2. ? (isInside ? 1.5 * Geom::pi() - theta : theta - Geom::pi() / 2.)
1924  : (isInside ? Geom::pi() / 2. - theta : theta + Geom::pi() / 2);
1925  //this is a normVector from the cone to the point
1926  Vector norm;
1927  norm.setRThetaPhi(fabs(dist), normTheta, normPhi);
1928  double curP = sv.p3.mag();
1929  double cosp3theta = sv.p3.z() / curP;
1930  if (cosp3theta > 1)
1931  cosp3theta = 1;
1932  if (cosp3theta < -1)
1933  cosp3theta = -1;
1934  double sineConeP = sinTheta * cosp3theta - cosTheta * sqrt(1. - cosp3theta * cosp3theta);
1935 
1936  double sinSolid = norm.dot(sv.p3) / (fabs(dist) * curP);
1937  tanDist = fabs(sinSolid) > fabs(sineConeP) ? dist / fabs(sinSolid) : dist / fabs(sineConeP);
1938 
1939  refDirection = sinSolid > 0 ? oppositeToMomentum : alongMomentum;
1940  if (debug_) {
1941  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Cone pars: theta " << theta
1942  << " normTheta " << norm.theta() << " p3Theta " << sv.p3.theta() << " sinD: " << sineConeP
1943  << " sinSolid " << sinSolid;
1944  }
1945  if (debug_) {
1946  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
1947  << "refToDest:toCone the point is " << (isInside ? "in" : "out") << "side the cone"
1948  << std::endl;
1949  }
1950  }
1951  result = SteppingHelixStateInfo::OK;
1952  } break;
1953  // case CYLINDER_DT:
1954  // break;
1955  case PATHL_DT: {
1956  double curS = fabs(sv.path());
1957  dist = pars[PATHL_P] - curS;
1958  if (fabs(dist) > fastSkipDist) {
1959  result = SteppingHelixStateInfo::INACC;
1960  break;
1961  }
1962  refDirection = pars[PATHL_P] > 0 ? alongMomentum : oppositeToMomentum;
1963  tanDist = dist;
1964  result = SteppingHelixStateInfo::OK;
1965  } break;
1966  case POINT_PCA_DT: {
1967  Point pDest(pars[0], pars[1], pars[2]);
1968  double curP = sv.p3.mag();
1969  dist = (sv.r3 - pDest).mag() + 1e-24; //add a small number to avoid 1/0
1970  tanDist = (sv.r3.dot(sv.p3) - pDest.dot(sv.p3)) / curP;
1971  //account for bending in magnetic field (quite approximate)
1972  double b0 = sv.bf.mag();
1973  if (b0 > 1.5e-6) {
1974  double p0 = curP;
1975  double kVal = 2.99792458e-3 * sv.q / p0 * b0;
1976  double aVal = fabs(dist * kVal);
1977  tanDist *= 1. / (1. + aVal);
1978  if (debug_)
1979  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "corrected by aVal " << aVal
1980  << " to " << tanDist;
1981  }
1982  refDirection = tanDist < 0 ? alongMomentum : oppositeToMomentum;
1983  result = SteppingHelixStateInfo::OK;
1984  } break;
1985  case LINE_PCA_DT: {
1986  Point rLine(pars[0], pars[1], pars[2]);
1987  Vector dLine(pars[3], pars[4], pars[5]);
1988  dLine = (dLine - rLine);
1989  dLine *= 1. / dLine.mag();
1990  if (debug_)
1991  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "dLine " << dLine;
1992 
1993  Vector dR = sv.r3 - rLine;
1994  double curP = sv.p3.mag();
1995  Vector dRPerp = dR - dLine * (dR.dot(dLine));
1996  if (debug_)
1997  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "dRperp " << dRPerp;
1998 
1999  dist = dRPerp.mag() + 1e-24; //add a small number to avoid 1/0
2000  tanDist = dRPerp.dot(sv.p3) / curP;
2001  //angle wrt line
2002  double cosAlpha2 = dLine.dot(sv.p3) / curP;
2003  cosAlpha2 *= cosAlpha2;
2004  tanDist *= 1. / sqrt(fabs(1. - cosAlpha2) + 1e-96);
2005  //correct for dPhi in magnetic field: this isn't made quite right here
2006  //(the angle between the line and the trajectory plane is neglected .. conservative)
2007  double b0 = sv.bf.mag();
2008  if (b0 > 1.5e-6) {
2009  double p0 = curP;
2010  double kVal = 2.99792458e-3 * sv.q / p0 * b0;
2011  double aVal = fabs(dist * kVal);
2012  tanDist *= 1. / (1. + aVal);
2013  if (debug_)
2014  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "corrected by aVal " << aVal
2015  << " to " << tanDist;
2016  }
2017  refDirection = tanDist < 0 ? alongMomentum : oppositeToMomentum;
2018  result = SteppingHelixStateInfo::OK;
2019  } break;
2020  default: {
2021  //some large number
2022  dist = 1e12;
2023  tanDist = 1e12;
2024  refDirection = anyDirection;
2025  result = SteppingHelixStateInfo::NOT_IMPLEMENTED;
2026  } break;
2027  }
2028 
2029  double tanDistConstrained = tanDist;
2030  if (fabs(tanDist) > 4. * fabs(dist))
2031  tanDistConstrained *= tanDist == 0 ? 0 : fabs(dist / tanDist * 4.);
2032 
2033  if (debug_) {
2034  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "refToDest input: dest" << dest
2035  << " pars[]: ";
2036  for (int i = 0; i < 6; i++) {
2037  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << ", " << i << " " << pars[i];
2038  }
2039  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << std::endl;
2040  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "refToDest output: "
2041  << "\t dist" << dist << "\t tanDist " << tanDist << "\t tanDistConstr " << tanDistConstrained
2042  << "\t refDirection " << refDirection << std::endl;
2043  }
2044  tanDist = tanDistConstrained;
2045 
2046  return result;
2047 }
mps_fire.i
i
Definition: mps_fire.py:428
metname
const std::string metname
Definition: MuonSeedOrcaPatternRecognition.cc:40
Vector
ROOT::Math::Plane3D::Vector Vector
Definition: EcalHitMaker.cc:29
mag
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
Definition: Basic3DVectorLD.h:127
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
theta
Geom::Theta< T > theta() const
Definition: Basic3DVectorLD.h:150
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
LogTrace
#define LogTrace(id)
Definition: MessageLogger.h:234
mps_fire.result
tuple result
Definition: mps_fire.py:311
SteppingHelixPropagator::Result
SteppingHelixStateInfo::Result Result
Definition: SteppingHelixPropagator.h:42
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
Point
Structure Point Contains parameters of Gaussian fits to DMRs.
Definition: DMRtrends.cc:57
b0
static constexpr float b0
Definition: L1EGammaCrystalsEmulatorProducer.cc:83
Geom::pi
constexpr double pi()
Definition: Pi.h:31
SteppingHelixPropagator::Result SteppingHelixPropagator::refToMagVolume ( const SteppingHelixPropagator::StateInfo &  sv,
PropagationDirection  dir,
double &  dist,
double &  tanDist,
double  fastSkipDist = 1e12,
bool  expectNewMagVolume = false,
double  maxStep = 1e12 
) const
protected

(Internals) determine distance and direction from the current position to the boundary of current mag volume

Definition at line 2049 of file SteppingHelixPropagator.cc.

References alongMomentum, anyDirection, CONE_DT, debug_, alignCSCRings::e, MagVolume::faces(), mps_fire::i, SteppingHelixStateInfo::INACC, MagVolume::inside(), dqmiolumiharvest::j, LogTrace, SteppingHelixStateInfo::magVol, metname, SteppingHelixStateInfo::NOT_IMPLEMENTED, SteppingHelixStateInfo::OK, Cone::openingAngle(), SteppingHelixStateInfo::p3, PLANE_DT, GloballyPositioned< T >::position(), SteppingHelixStateInfo::r3, Cylinder::radius(), RADIUS_DT, RADIUS_P, refToDest(), mps_fire::result, GloballyPositioned< T >::rotation(), jetcorrextractor::sign(), AlCaHLTBitMon_QueryRunRegistry::string, SteppingHelixStateInfo::UNDEFINED, UNDEFINED_DT, Cone::vertex(), SteppingHelixStateInfo::WRONG_VOLUME, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), PV3DBase< T, PVType, FrameType >::z(), TkRotation< T >::zx(), TkRotation< T >::zy(), and TkRotation< T >::zz().

Referenced by propagate().

2055  {
2056  static const std::string metname = "SteppingHelixPropagator";
2057  Result result = SteppingHelixStateInfo::NOT_IMPLEMENTED;
2058  const MagVolume* cVol = sv.magVol;
2059 
2060  if (cVol == nullptr)
2061  return result;
2062  const std::vector<VolumeSide>& cVolFaces(cVol->faces());
2063 
2064  double distToFace[6] = {0, 0, 0, 0, 0, 0};
2065  double tanDistToFace[6] = {0, 0, 0, 0, 0, 0};
2066  PropagationDirection refDirectionToFace[6] = {
2067  anyDirection, anyDirection, anyDirection, anyDirection, anyDirection, anyDirection};
2068  Result resultToFace[6] = {result, result, result, result, result, result};
2069  int iFDest = -1;
2070  int iDistMin = -1;
2071 
2072  unsigned int iFDestSorted[6] = {0, 0, 0, 0, 0, 0};
2073  int nDestSorted = 0;
2074  unsigned int nearParallels = 0;
2075 
2076  double curP = sv.p3.mag();
2077 
2078  if (debug_) {
2079  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Trying volume "
2080  << " with " << cVolFaces.size() << " faces" << std::endl;
2081  }
2082 
2083  unsigned int nFaces = cVolFaces.size();
2084  for (unsigned int iFace = 0; iFace < nFaces; ++iFace) {
2085  if (iFace > 5) {
2086  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Too many faces" << std::endl;
2087  }
2088  if (debug_) {
2089  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Start with face " << iFace
2090  << std::endl;
2091  }
2092  // const Plane* cPlane = dynamic_cast<const Plane*>(&cVolFaces[iFace].surface());
2093  // const Cylinder* cCyl = dynamic_cast<const Cylinder*>(&cVolFaces[iFace].surface());
2094  // const Cone* cCone = dynamic_cast<const Cone*>(&cVolFaces[iFace].surface());
2095  const Surface* cPlane = nullptr; //only need to know the loc->glob transform
2096  const Cylinder* cCyl = nullptr;
2097  const Cone* cCone = nullptr;
2098  auto& iSurface = cVolFaces[iFace].surface();
2099  if (typeid(iSurface) == typeid(const Plane&)) {
2100  cPlane = &cVolFaces[iFace].surface();
2101  } else if (typeid(iSurface) == typeid(const Cylinder&)) {
2102  cCyl = dynamic_cast<const Cylinder*>(&cVolFaces[iFace].surface());
2103  } else if (typeid(iSurface) == typeid(const Cone&)) {
2104  cCone = dynamic_cast<const Cone*>(&cVolFaces[iFace].surface());
2105  } else {
2106  edm::LogWarning(metname) << std::setprecision(17) << std::setw(20) << std::scientific
2107  << "Could not cast a volume side surface to a known type" << std::endl;
2108  }
2109 
2110  if (debug_) {
2111  if (cPlane != nullptr)
2112  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "The face is a plane at "
2113  << cPlane << std::endl;
2114  if (cCyl != nullptr)
2115  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "The face is a cylinder at "
2116  << cCyl << std::endl;
2117  }
2118 
2119  double pars[6] = {0, 0, 0, 0, 0, 0};
2120  DestType dType = UNDEFINED_DT;
2121  if (cPlane != nullptr) {
2122  Point rPlane(cPlane->position().x(), cPlane->position().y(), cPlane->position().z());
2123  // = cPlane->toGlobal(LocalVector(0,0,1.)); nPlane = nPlane.unit();
2124  Vector nPlane(cPlane->rotation().zx(), cPlane->rotation().zy(), cPlane->rotation().zz());
2125  nPlane /= nPlane.mag();
2126 
2127  pars[0] = rPlane.x();
2128  pars[1] = rPlane.y();
2129  pars[2] = rPlane.z();
2130  pars[3] = nPlane.x();
2131  pars[4] = nPlane.y();
2132  pars[5] = nPlane.z();
2133  dType = PLANE_DT;
2134  } else if (cCyl != nullptr) {
2135  if (debug_) {
2136  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Cylinder at "
2137  << cCyl->position() << " rorated by " << cCyl->rotation() << std::endl;
2138  }
2139  pars[RADIUS_P] = cCyl->radius();
2140  dType = RADIUS_DT;
2141  } else if (cCone != nullptr) {
2142  if (debug_) {
2143  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Cone at "
2144  << cCone->position() << " rorated by " << cCone->rotation() << " vertex at "
2145  << cCone->vertex() << " angle of " << cCone->openingAngle() << std::endl;
2146  }
2147  pars[0] = cCone->vertex().x();
2148  pars[1] = cCone->vertex().y();
2149  pars[2] = cCone->vertex().z();
2150  pars[3] = cCone->openingAngle();
2151  dType = CONE_DT;
2152  } else {
2153  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Unknown surface" << std::endl;
2154  resultToFace[iFace] = SteppingHelixStateInfo::UNDEFINED;
2155  continue;
2156  }
2157  resultToFace[iFace] =
2158  refToDest(dType, sv, pars, distToFace[iFace], tanDistToFace[iFace], refDirectionToFace[iFace], fastSkipDist);
2159 
2160  if (resultToFace[iFace] != SteppingHelixStateInfo::OK) {
2161  if (resultToFace[iFace] == SteppingHelixStateInfo::INACC)
2162  result = SteppingHelixStateInfo::INACC;
2163  }
2164 
2165  //keep those in right direction for later use
2166  if (resultToFace[iFace] == SteppingHelixStateInfo::OK) {
2167  double invDTFPosiF = 1. / (1e-32 + fabs(tanDistToFace[iFace]));
2168  double dSlope = fabs(distToFace[iFace]) * invDTFPosiF;
2169  double maxStepL = maxStep > 100 ? 100 : maxStep;
2170  if (maxStepL < 10)
2171  maxStepL = 10.;
2172  bool isNearParallel =
2173  fabs(tanDistToFace[iFace]) + 100. * curP * dSlope < maxStepL //
2174  //a better choice is to use distance to next check of mag volume instead of 100cm; the last is for ~1.5arcLength(4T)+tandistance< maxStep
2175  && dSlope < 0.15; //
2176  if (refDirectionToFace[iFace] == dir || isNearParallel) {
2177  if (isNearParallel)
2178  nearParallels++;
2179  iFDestSorted[nDestSorted] = iFace;
2180  nDestSorted++;
2181  }
2182  }
2183 
2184  //pick a shortest distance here (right dir only for now)
2185  if (refDirectionToFace[iFace] == dir) {
2186  if (iDistMin == -1)
2187  iDistMin = iFace;
2188  else if (fabs(distToFace[iFace]) < fabs(distToFace[iDistMin]))
2189  iDistMin = iFace;
2190  }
2191  if (debug_)
2192  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << cVol << " " << iFace << " "
2193  << tanDistToFace[iFace] << " " << distToFace[iFace] << " " << refDirectionToFace[iFace] << " "
2194  << dir << std::endl;
2195  }
2196 
2197  for (int i = 0; i < nDestSorted; ++i) {
2198  int iMax = nDestSorted - i - 1;
2199  for (int j = 0; j < nDestSorted - i; ++j) {
2200  if (fabs(tanDistToFace[iFDestSorted[j]]) > fabs(tanDistToFace[iFDestSorted[iMax]])) {
2201  iMax = j;
2202  }
2203  }
2204  int iTmp = iFDestSorted[nDestSorted - i - 1];
2205  iFDestSorted[nDestSorted - i - 1] = iFDestSorted[iMax];
2206  iFDestSorted[iMax] = iTmp;
2207  }
2208 
2209  if (debug_) {
2210  for (int i = 0; i < nDestSorted; ++i) {
2211  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << cVol << " " << i << " "
2212  << iFDestSorted[i] << " " << tanDistToFace[iFDestSorted[i]] << std::endl;
2213  }
2214  }
2215 
2216  //now go from the shortest to the largest distance hoping to get a point in the volume.
2217  //other than in case of a near-parallel travel this should stop after the first try
2218 
2219  for (int i = 0; i < nDestSorted; ++i) {
2220  iFDest = iFDestSorted[i];
2221 
2222  double sign = dir == alongMomentum ? 1. : -1.;
2223  Point gPointEst(sv.r3);
2224  Vector lDelta(sv.p3);
2225  lDelta *= sign / curP * fabs(distToFace[iFDest]);
2226  gPointEst += lDelta;
2227  if (debug_) {
2228  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Linear est point " << gPointEst
2229  << " for iFace " << iFDest << std::endl;
2230  }
2231  GlobalPoint gPointEstNorZ(gPointEst.x(), gPointEst.y(), gPointEst.z());
2232  if (cVol->inside(gPointEstNorZ)) {
2233  if (debug_) {
2234  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
2235  << "The point is inside the volume" << std::endl;
2236  }
2237 
2238  result = SteppingHelixStateInfo::OK;
2239  dist = distToFace[iFDest];
2240  tanDist = tanDistToFace[iFDest];
2241  if (debug_) {
2242  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
2243  << "Got a point near closest boundary -- face " << iFDest << std::endl;
2244  }
2245  break;
2246  }
2247  }
2248 
2249  if (result != SteppingHelixStateInfo::OK && expectNewMagVolume) {
2250  double sign = dir == alongMomentum ? 1. : -1.;
2251 
2252  //check if it's a wrong volume situation
2253  if (nDestSorted - nearParallels > 0)
2254  result = SteppingHelixStateInfo::WRONG_VOLUME;
2255  else {
2256  //get here if all faces in the corr direction were skipped
2257  Point gPointEst(sv.r3);
2258  double lDist = iDistMin == -1 ? fastSkipDist : fabs(distToFace[iDistMin]);
2259  if (lDist > fastSkipDist)
2260  lDist = fastSkipDist;
2261  Vector lDelta(sv.p3);
2262  lDelta *= sign / curP * lDist;
2263  gPointEst += lDelta;
2264  if (debug_) {
2265  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
2266  << "Linear est point to shortest dist " << gPointEst << " for iFace " << iDistMin
2267  << " at distance " << lDist * sign << std::endl;
2268  }
2269  GlobalPoint gPointEstNorZ(gPointEst.x(), gPointEst.y(), gPointEst.z());
2270  if (cVol->inside(gPointEstNorZ)) {
2271  if (debug_) {
2272  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
2273  << "The point is inside the volume" << std::endl;
2274  }
2275 
2276  } else {
2277  result = SteppingHelixStateInfo::WRONG_VOLUME;
2278  }
2279  }
2280 
2281  if (result == SteppingHelixStateInfo::WRONG_VOLUME) {
2282  dist = sign * 0.05;
2283  tanDist = dist * 1.01;
2284  if (debug_) {
2285  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
2286  << "Wrong volume located: return small dist, tandist" << std::endl;
2287  }
2288  }
2289  }
2290 
2291  if (result == SteppingHelixStateInfo::INACC) {
2292  if (debug_)
2293  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "All faces are too far"
2294  << std::endl;
2295  } else if (result == SteppingHelixStateInfo::WRONG_VOLUME) {
2296  if (debug_)
2297  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Appear to be in a wrong volume"
2298  << std::endl;
2299  } else if (result != SteppingHelixStateInfo::OK) {
2300  if (debug_)
2301  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Something else went wrong"
2302  << std::endl;
2303  }
2304 
2305  return result;
2306 }
Cone
Definition: Cone.h:17
mps_fire.i
i
Definition: mps_fire.py:428
SteppingHelixPropagator::refToDest
Result refToDest(DestType dest, const SteppingHelixPropagator::StateInfo &sv, const double pars[6], double &dist, double &tanDist, PropagationDirection &refDirection, double fastSkipDist=1e12) const
(Internals) determine distance and direction from the current position to the plane ...
Definition: SteppingHelixPropagator.cc:1781
metname
const std::string metname
Definition: MuonSeedOrcaPatternRecognition.cc:40
Vector
ROOT::Math::Plane3D::Vector Vector
Definition: EcalHitMaker.cc:29
jetcorrextractor::sign
double sign(double x)
Definition: JetCorrExtractorT.h:31
PV3DBase::y
T y() const
Definition: PV3DBase.h:60
PropagationDirection
PropagationDirection
Definition: PropagationDirection.h:4
Plane
Definition: Plane.h:16
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
LogTrace
#define LogTrace(id)
Definition: MessageLogger.h:234
mps_fire.result
tuple result
Definition: mps_fire.py:311
TkRotation::zx
T zx() const
Definition: extTkRotation.h:199
Cone::vertex
GlobalPoint vertex() const
Global position of the cone vertex.
Definition: Cone.h:47
TkRotation::zz
T zz() const
Definition: extTkRotation.h:201
MagVolume::inside
virtual bool inside(const GlobalPoint &gp, double tolerance=0.) const =0
SteppingHelixPropagator::Result
SteppingHelixStateInfo::Result Result
Definition: SteppingHelixPropagator.h:42
Cylinder::radius
Scalar radius() const
Radius of the cylinder.
Definition: Cylinder.h:64
DeadROC_duringRun.dir
string dir
Definition: DeadROC_duringRun.py:23
PV3DBase::z
T z() const
Definition: PV3DBase.h:61
TkRotation::zy
T zy() const
Definition: extTkRotation.h:200
MagVolume::faces
virtual const std::vector< VolumeSide > & faces() const =0
Access to volume faces.
Point
Structure Point Contains parameters of Gaussian fits to DMRs.
Definition: DMRtrends.cc:57
GloballyPositioned::rotation
const RotationType & rotation() const
Definition: GloballyPositioned.h:38
edm::LogWarning
Log< level::Warning, false > LogWarning
Definition: MessageLogger.h:122
PV3DBase::x
T x() const
Definition: PV3DBase.h:59
GloballyPositioned::position
const PositionType & position() const
Definition: GloballyPositioned.h:36
Cone::openingAngle
Geom::Theta< float > openingAngle() const
Angle of the cone.
Definition: Cone.h:50
SteppingHelixPropagator::Result SteppingHelixPropagator::refToMatVolume ( const SteppingHelixPropagator::StateInfo &  sv,
PropagationDirection  dir,
double &  dist,
double &  tanDist,
double  fastSkipDist = 1e12 
) const
protected

Definition at line 2308 of file SteppingHelixPropagator.cc.

References alongMomentum, anyDirection, CONE_DT, debug_, alignCSCRings::e, SteppingHelixStateInfo::INACC, LogTrace, metname, SteppingHelixStateInfo::NOT_IMPLEMENTED, SteppingHelixStateInfo::OK, SteppingHelixStateInfo::p3, Geom::pi(), PLANE_DT, SteppingHelixStateInfo::r3, RADIUS_DT, RADIUS_P, refToDest(), mps_fire::result, SteppingHelixStateInfo::VolumeBounds::rMax, SteppingHelixStateInfo::VolumeBounds::rMin, SteppingHelixStateInfo::rzLims, jetcorrextractor::sign(), AlCaHLTBitMon_QueryRunRegistry::string, funct::tan(), SteppingHelixStateInfo::VolumeBounds::th1, SteppingHelixStateInfo::VolumeBounds::th2, UNDEFINED_DT, SteppingHelixStateInfo::VolumeBounds::zMax, and SteppingHelixStateInfo::VolumeBounds::zMin.

Referenced by propagate().

2312  {
2313  static const std::string metname = "SteppingHelixPropagator";
2314  Result result = SteppingHelixStateInfo::NOT_IMPLEMENTED;
2315 
2316  double parLim[6] = {sv.rzLims.rMin, sv.rzLims.rMax, sv.rzLims.zMin, sv.rzLims.zMax, sv.rzLims.th1, sv.rzLims.th2};
2317 
2318  double distToFace[4] = {0, 0, 0, 0};
2319  double tanDistToFace[4] = {0, 0, 0, 0};
2320  PropagationDirection refDirectionToFace[4] = {anyDirection, anyDirection, anyDirection, anyDirection};
2321  Result resultToFace[4] = {result, result, result, result};
2322  int iFDest = -1;
2323 
2324  double curP = sv.p3.mag();
2325 
2326  for (unsigned int iFace = 0; iFace < 4; iFace++) {
2327  if (debug_) {
2328  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Start with mat face " << iFace
2329  << std::endl;
2330  }
2331 
2332  double pars[6] = {0, 0, 0, 0, 0, 0};
2333  DestType dType = UNDEFINED_DT;
2334  if (iFace > 1) {
2335  if (fabs(parLim[iFace + 2]) < 1e-6) { //plane
2336  if (sv.r3.z() < 0) {
2337  pars[0] = 0;
2338  pars[1] = 0;
2339  pars[2] = -parLim[iFace];
2340  pars[3] = 0;
2341  pars[4] = 0;
2342  pars[5] = 1;
2343  } else {
2344  pars[0] = 0;
2345  pars[1] = 0;
2346  pars[2] = parLim[iFace];
2347  pars[3] = 0;
2348  pars[4] = 0;
2349  pars[5] = 1;
2350  }
2351  dType = PLANE_DT;
2352  } else {
2353  if (sv.r3.z() > 0) {
2354  pars[0] = 0;
2355  pars[1] = 0;
2356  pars[2] = parLim[iFace];
2357  pars[3] = Geom::pi() / 2. - parLim[iFace + 2];
2358  } else {
2359  pars[0] = 0;
2360  pars[1] = 0;
2361  pars[2] = -parLim[iFace];
2362  pars[3] = Geom::pi() / 2. + parLim[iFace + 2];
2363  }
2364  dType = CONE_DT;
2365  }
2366  } else {
2367  pars[RADIUS_P] = parLim[iFace];
2368  dType = RADIUS_DT;
2369  }
2370 
2371  resultToFace[iFace] =
2372  refToDest(dType, sv, pars, distToFace[iFace], tanDistToFace[iFace], refDirectionToFace[iFace], fastSkipDist);
2373 
2374  if (resultToFace[iFace] != SteppingHelixStateInfo::OK) {
2375  if (resultToFace[iFace] == SteppingHelixStateInfo::INACC)
2376  result = SteppingHelixStateInfo::INACC;
2377  continue;
2378  }
2379  if (refDirectionToFace[iFace] == dir || fabs(distToFace[iFace]) < 2e-2 * fabs(tanDistToFace[iFace])) {
2380  double sign = dir == alongMomentum ? 1. : -1.;
2381  Point gPointEst(sv.r3);
2382  Vector lDelta(sv.p3);
2383  lDelta *= sign * fabs(distToFace[iFace]) / curP;
2384  gPointEst += lDelta;
2385  if (debug_) {
2386  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific << "Linear est point "
2387  << gPointEst << std::endl;
2388  }
2389  double lZ = fabs(gPointEst.z());
2390  double lR = gPointEst.perp();
2391  double tan4 = parLim[4] == 0 ? 0 : tan(parLim[4]);
2392  double tan5 = parLim[5] == 0 ? 0 : tan(parLim[5]);
2393  if ((lZ - parLim[2]) > lR * tan4 && (lZ - parLim[3]) < lR * tan5 && lR > parLim[0] && lR < parLim[1]) {
2394  if (debug_) {
2395  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
2396  << "The point is inside the volume" << std::endl;
2397  }
2398  //OK, guessed a point still inside the volume
2399  if (iFDest == -1) {
2400  iFDest = iFace;
2401  } else {
2402  if (fabs(tanDistToFace[iFDest]) > fabs(tanDistToFace[iFace])) {
2403  iFDest = iFace;
2404  }
2405  }
2406  } else {
2407  if (debug_) {
2408  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
2409  << "The point is NOT inside the volume" << std::endl;
2410  }
2411  }
2412  }
2413  }
2414  if (iFDest != -1) {
2415  result = SteppingHelixStateInfo::OK;
2416  dist = distToFace[iFDest];
2417  tanDist = tanDistToFace[iFDest];
2418  if (debug_) {
2419  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
2420  << "Got a point near closest boundary -- face " << iFDest << std::endl;
2421  }
2422  } else {
2423  if (debug_) {
2424  LogTrace(metname) << std::setprecision(17) << std::setw(20) << std::scientific
2425  << "Failed to find a dest point inside the volume" << std::endl;
2426  }
2427  }
2428 
2429  return result;
2430 }
SteppingHelixPropagator::refToDest
Result refToDest(DestType dest, const SteppingHelixPropagator::StateInfo &sv, const double pars[6], double &dist, double &tanDist, PropagationDirection &refDirection, double fastSkipDist=1e12) const
(Internals) determine distance and direction from the current position to the plane ...
Definition: SteppingHelixPropagator.cc:1781
metname
const std::string metname
Definition: MuonSeedOrcaPatternRecognition.cc:40
Vector
ROOT::Math::Plane3D::Vector Vector
Definition: EcalHitMaker.cc:29
jetcorrextractor::sign
double sign(double x)
Definition: JetCorrExtractorT.h:31
PropagationDirection
PropagationDirection
Definition: PropagationDirection.h:4
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
LogTrace
#define LogTrace(id)
Definition: MessageLogger.h:234
mps_fire.result
tuple result
Definition: mps_fire.py:311
SteppingHelixPropagator::Result
SteppingHelixStateInfo::Result Result
Definition: SteppingHelixPropagator.h:42
DeadROC_duringRun.dir
string dir
Definition: DeadROC_duringRun.py:23
funct::tan
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
Point
Structure Point Contains parameters of Gaussian fits to DMRs.
Definition: DMRtrends.cc:57
Geom::pi
constexpr double pi()
Definition: Pi.h:31
void SteppingHelixPropagator::setDebug ( bool  debug)
inline

Switch debug printouts (to cout) .. very verbose.

Definition at line 121 of file SteppingHelixPropagator.h.

References debug, and debug_.

121 { debug_ = debug; }
debug
#define debug
Definition: HDRShower.cc:19
void SteppingHelixPropagator::setEndcapShiftsInZPosNeg ( double  valPos,
double  valNeg 
)
inline

set shifts in Z for endcap pieces (includes EE, HE, ME, YE)

Definition at line 164 of file SteppingHelixPropagator.h.

References ecShiftNeg_, and ecShiftPos_.

164  {
165  ecShiftPos_ = valPos;
166  ecShiftNeg_ = valNeg;
167  }
void SteppingHelixPropagator::setIState ( const SteppingHelixStateInfo &  sStart,
StateArray &  svBuff,
int &  nPoints 
) const
protected

(Internals) Init starting point

Definition at line 296 of file SteppingHelixPropagator.cc.

References cIndex_(), SteppingHelixStateInfo::covCurv, SteppingHelixStateInfo::hasErrorPropagated_, SteppingHelixStateInfo::isComplete, loadState(), noErrorPropagation_, SteppingHelixStateInfo::p3, Propagator::propagationDirection(), SteppingHelixStateInfo::q, and SteppingHelixStateInfo::r3.

Referenced by propagate(), and propagateWithPath().

296  {
297  nPoints = 0;
298  svBuf[cIndex_(nPoints)] = sStart; //do this anyways to have a fresh start
299  if (sStart.isComplete) {
300  svBuf[cIndex_(nPoints)] = sStart;
301  nPoints++;
302  } else {
303  loadState(svBuf[cIndex_(nPoints)], sStart.p3, sStart.r3, sStart.q, propagationDirection(), sStart.covCurv);
304  nPoints++;
305  }
306  svBuf[cIndex_(0)].hasErrorPropagated_ = sStart.hasErrorPropagated_ & !noErrorPropagation_;
307 }
SteppingHelixStateInfo::covCurv
AlgebraicSymMatrix55 covCurv
Definition: SteppingHelixStateInfo.h:87
Propagator::propagationDirection
virtual PropagationDirection propagationDirection() const final
Definition: Propagator.h:139
SteppingHelixPropagator::cIndex_
int cIndex_(int ind) const
(Internals) circular index for array of transients
Definition: SteppingHelixPropagator.cc:1773
SteppingHelixPropagator::loadState
void loadState(SteppingHelixPropagator::StateInfo &svCurrent, const SteppingHelixPropagator::Vector &p3, const SteppingHelixPropagator::Point &r3, int charge, PropagationDirection dir, const AlgebraicSymMatrix55 &covCurv) const
Definition: SteppingHelixPropagator.cc:657
void SteppingHelixPropagator::setMaterialMode ( bool  noMaterial)
inline

Switch for material effects mode: no material effects if true.

Definition at line 124 of file SteppingHelixPropagator.h.

References noMaterialMode_.

Referenced by MuonSimHitProducer::beginRun(), AlCaHOCalibProducer::fillHOStore(), and PropagateToCal::propagate().

124 { noMaterialMode_ = noMaterial; }
void SteppingHelixPropagator::setNoErrorPropagation ( bool  noErrorPropagation)
inline

Force no error propagation.

Definition at line 127 of file SteppingHelixPropagator.h.

References noErrorPropagation_.

Referenced by PropagateToCal::propagate().

127 { noErrorPropagation_ = noErrorPropagation; }
void SteppingHelixPropagator::setRep ( SteppingHelixPropagator::Basis &  rep,
const SteppingHelixPropagator::Vector &  tau 
) const
protected

Set/compute basis vectors for local coordinates at current step (called by incrementState)

Definition at line 870 of file SteppingHelixPropagator.cc.

References SteppingHelixPropagator::Basis::lX, SteppingHelixPropagator::Basis::lY, SteppingHelixPropagator::Basis::lZ, and metsig::tau.

Referenced by makeAtomStep().

871  {
872  Vector zRep(0., 0., 1.);
873  rep.lX = tau;
874  rep.lY = zRep.cross(tau);
875  rep.lY *= 1. / tau.perp();
876  rep.lZ = rep.lX.cross(rep.lY);
877 }
Vector
ROOT::Math::Plane3D::Vector Vector
Definition: EcalHitMaker.cc:29
void SteppingHelixPropagator::setReturnTangentPlane ( bool  val)
inline

flag to return tangent plane for non-plane input

Definition at line 142 of file SteppingHelixPropagator.h.

References returnTangentPlane_, and hgcalPerformanceValidation::val.

142 { returnTangentPlane_ = val; }
void SteppingHelixPropagator::setSendLogWarning ( bool  val)
inline

flag to send LogWarning on failures

Definition at line 145 of file SteppingHelixPropagator.h.

References sendLogWarning_, and hgcalPerformanceValidation::val.

145 { sendLogWarning_ = val; }
void SteppingHelixPropagator::setUseInTeslaFromMagField ( bool  val)
inline

force getting field value from MagneticField, not the geometric one

Definition at line 161 of file SteppingHelixPropagator.h.

References useInTeslaFromMagField_, and hgcalPerformanceValidation::val.

161 { useInTeslaFromMagField_ = val; }
void SteppingHelixPropagator::setUseIsYokeFlag ( bool  val)
inline

(temporary?) flag to identify if the volume is yoke based on MagVolume internals works if useMatVolumes_ is also true

Definition at line 149 of file SteppingHelixPropagator.h.

References useIsYokeFlag_, and hgcalPerformanceValidation::val.

149 { useIsYokeFlag_ = val; }
void SteppingHelixPropagator::setUseMagVolumes ( bool  val)
inline

Switch to using MagneticField Volumes .. as in VolumeBasedMagneticField.

Definition at line 136 of file SteppingHelixPropagator.h.

References useMagVolumes_, and hgcalPerformanceValidation::val.

136 { useMagVolumes_ = val; }
void SteppingHelixPropagator::setUseMatVolumes ( bool  val)
inline

Switch to using Material Volumes .. internally defined for now.

Definition at line 139 of file SteppingHelixPropagator.h.

References useMatVolumes_, and hgcalPerformanceValidation::val.

139 { useMatVolumes_ = val; }
void SteppingHelixPropagator::setUseTuningForL2Speed ( bool  val)
inline

will use bigger steps ~tuned for good-enough L2/Standalone reco use this in hope of a speed-up

Definition at line 153 of file SteppingHelixPropagator.h.

References useTuningForL2Speed_, and hgcalPerformanceValidation::val.

153 { useTuningForL2Speed_ = val; }
void SteppingHelixPropagator::setVBFPointer ( const VolumeBasedMagneticField *  val)
inline

set VolumBasedField pointer allows to have geometry description in uniformField scenario only important/relevant in the barrel region

Definition at line 158 of file SteppingHelixPropagator.h.

References hgcalPerformanceValidation::val, and vbField_.

158 { vbField_ = val; }
SteppingHelixPropagator::vbField_
const VolumeBasedMagneticField * vbField_
Definition: SteppingHelixPropagator.h:266

Member Data Documentation

bool SteppingHelixPropagator::applyRadX0Correction_
private

Definition at line 271 of file SteppingHelixPropagator.h.

Referenced by applyRadX0Correction(), makeAtomStep(), and SteppingHelixPropagator().

bool SteppingHelixPropagator::debug_
private

Definition at line 268 of file SteppingHelixPropagator.h.

Referenced by getNextState(), isYokeVolume(), loadState(), makeAtomStep(), propagate(), refToDest(), refToMagVolume(), refToMatVolume(), setDebug(), and SteppingHelixPropagator().

double SteppingHelixPropagator::defaultStep_
private

Definition at line 280 of file SteppingHelixPropagator.h.

Referenced by propagate(), and SteppingHelixPropagator().

double SteppingHelixPropagator::ecShiftNeg_
private

Definition at line 283 of file SteppingHelixPropagator.h.

Referenced by getDeDx(), setEndcapShiftsInZPosNeg(), and SteppingHelixPropagator().

double SteppingHelixPropagator::ecShiftPos_
private

Definition at line 282 of file SteppingHelixPropagator.h.

Referenced by getDeDx(), setEndcapShiftsInZPosNeg(), and SteppingHelixPropagator().

const MagneticField* SteppingHelixPropagator::field_
private

Definition at line 265 of file SteppingHelixPropagator.h.

Referenced by getNextState(), loadState(), magneticField(), makeAtomStep(), propagate(), and SteppingHelixPropagator().

StateInfo SteppingHelixPropagator::invalidState_
private

Definition at line 263 of file SteppingHelixPropagator.h.

Referenced by propagate().

const int SteppingHelixPropagator::MAX_POINTS = 7
staticprotected

Definition at line 171 of file SteppingHelixPropagator.h.

Referenced by cIndex_(), and initStateArraySHPSpecific().

const int SteppingHelixPropagator::MAX_STEPS = 10000
staticprivate

Definition at line 259 of file SteppingHelixPropagator.h.

Referenced by propagate().

bool SteppingHelixPropagator::noErrorPropagation_
private

Definition at line 270 of file SteppingHelixPropagator.h.

Referenced by initStateArraySHPSpecific(), setIState(), setNoErrorPropagation(), and SteppingHelixPropagator().

bool SteppingHelixPropagator::noMaterialMode_
private

Definition at line 269 of file SteppingHelixPropagator.h.

Referenced by getDeDx(), setMaterialMode(), and SteppingHelixPropagator().

bool SteppingHelixPropagator::returnTangentPlane_
private

Definition at line 276 of file SteppingHelixPropagator.h.

Referenced by propagateWithPath(), setReturnTangentPlane(), and SteppingHelixPropagator().

bool SteppingHelixPropagator::sendLogWarning_
private

Definition at line 277 of file SteppingHelixPropagator.h.

Referenced by propagate(), propagateWithPath(), setSendLogWarning(), and SteppingHelixPropagator().

const AlgebraicSymMatrix55 SteppingHelixPropagator::unit55_
private

Definition at line 267 of file SteppingHelixPropagator.h.

Referenced by makeAtomStep().

bool SteppingHelixPropagator::useInTeslaFromMagField_
private

Definition at line 275 of file SteppingHelixPropagator.h.

Referenced by getNextState(), loadState(), makeAtomStep(), setUseInTeslaFromMagField(), and SteppingHelixPropagator().

bool SteppingHelixPropagator::useIsYokeFlag_
private

Definition at line 273 of file SteppingHelixPropagator.h.

Referenced by getDeDx(), getNextState(), loadState(), propagate(), setUseIsYokeFlag(), and SteppingHelixPropagator().

bool SteppingHelixPropagator::useMagVolumes_
private

Definition at line 272 of file SteppingHelixPropagator.h.

Referenced by getDeDx(), getNextState(), loadState(), makeAtomStep(), propagate(), setUseMagVolumes(), and SteppingHelixPropagator().

bool SteppingHelixPropagator::useMatVolumes_
private

Definition at line 274 of file SteppingHelixPropagator.h.

Referenced by propagate(), setUseMatVolumes(), and SteppingHelixPropagator().

bool SteppingHelixPropagator::useTuningForL2Speed_
private

Definition at line 278 of file SteppingHelixPropagator.h.

Referenced by propagate(), setUseTuningForL2Speed(), and SteppingHelixPropagator().

const VolumeBasedMagneticField* SteppingHelixPropagator::vbField_
private

Definition at line 266 of file SteppingHelixPropagator.h.

Referenced by getNextState(), loadState(), setVBFPointer(), and SteppingHelixPropagator().

Generated for CMSSW Reference Manual by doxygen 1.8.5