|
|
|
#include <SteppingHelixPropagator.h>
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) |
| virtual SteppingHelixPropagator * | clone () const |
| virtual const MagneticField * | magneticField () const |
| virtual FreeTrajectoryState | propagate (const FreeTrajectoryState &ftsStart, const GlobalPoint &pDest) const |
| Propagate to PCA to point given a starting point. | |
| const SteppingHelixStateInfo & | propagate (const SteppingHelixStateInfo &ftsStart, const Surface &sDest) const |
| Propagate to Plane given a starting point. | |
| const SteppingHelixStateInfo & | propagate (const SteppingHelixStateInfo &ftsStart, const Plane &pDest) const |
| virtual FreeTrajectoryState | propagate (const FreeTrajectoryState &ftsStart, const GlobalPoint &pDest1, const GlobalPoint &pDest2) const |
| Propagate to PCA to a line (given by 2 points) given a starting point. | |
| const SteppingHelixStateInfo & | propagate (const SteppingHelixStateInfo &ftsStart, const Cylinder &cDest) const |
| Propagate to Cylinder given a starting point (a Cylinder is assumed to be positioned at 0,0,0) | |
| const SteppingHelixStateInfo & | propagate (const SteppingHelixStateInfo &ftsStart, const GlobalPoint &pDest) const |
| Propagate to PCA to point given a starting point. | |
| virtual FreeTrajectoryState | propagate (const FreeTrajectoryState &ftsStart, const reco::BeamSpot &beamSpot) const |
| Propagate to PCA to a line determined by BeamSpot position and slope given a starting point. | |
| const SteppingHelixStateInfo & | propagate (const SteppingHelixStateInfo &ftsStart, const GlobalPoint &pDest1, const GlobalPoint &pDest2) const |
| Propagate to PCA to a line (given by 2 points) given a starting point. | |
| virtual TrajectoryStateOnSurface | propagate (const FreeTrajectoryState &ftsStart, const Plane &pDest) const |
| Propagate to Plane given a starting point. | |
| virtual TrajectoryStateOnSurface | propagate (const FreeTrajectoryState &ftsStart, const Cylinder &cDest) const |
| Propagate to Cylinder given a starting point (a Cylinder is assumed to be positioned at 0,0,0) | |
| virtual std::pair < TrajectoryStateOnSurface, double > | propagateWithPath (const FreeTrajectoryState &ftsStart, const Plane &pDest) const |
| virtual std::pair < TrajectoryStateOnSurface, double > | propagateWithPath (const FreeTrajectoryState &ftsStart, const Cylinder &cDest) const |
| virtual std::pair < FreeTrajectoryState, double > | propagateWithPath (const FreeTrajectoryState &ftsStart, const GlobalPoint &pDest) const |
| Propagate to PCA to point given a starting point. | |
| virtual std::pair < FreeTrajectoryState, double > | propagateWithPath (const FreeTrajectoryState &ftsStart, const GlobalPoint &pDest1, const GlobalPoint &pDest2) const |
| Propagate to PCA to a line (given by 2 points) given a starting point. | |
| virtual std::pair < FreeTrajectoryState, double > | propagateWithPath (const FreeTrajectoryState &ftsStart, const reco::BeamSpot &beamSpot) const |
| Propagate to PCA to a line (given by beamSpot position and slope) given a starting point. | |
| void | setDebug (bool debug) |
| Switch debug printouts (to cout) .. very verbose. | |
| void | setEndcapShiftsInZPosNeg (double valPos, double valNeg) |
| set shifts in Z for endcap pieces (includes EE, HE, ME, YE) | |
| void | setMaterialMode (bool noMaterial) |
| Switch for material effects mode: no material effects if true. | |
| void | setNoErrorPropagation (bool noErrorPropagation) |
| Force no error propagation. | |
| void | setReturnTangentPlane (bool val) |
| flag to return tangent plane for non-plane input | |
| void | setSendLogWarning (bool val) |
| flag to send LogWarning on failures | |
| void | setUseInTeslaFromMagField (bool val) |
| force getting field value from MagneticField, not the geometric one | |
| void | setUseIsYokeFlag (bool val) |
| void | setUseMagVolumes (bool val) |
| Switch to using MagneticField Volumes .. as in VolumeBasedMagneticField. | |
| void | setUseMatVolumes (bool val) |
| Switch to using Material Volumes .. internally defined for now. | |
| void | setUseTuningForL2Speed (bool val) |
| void | setVBFPointer (const VolumeBasedMagneticField *val) |
| SteppingHelixPropagator (const MagneticField *field, PropagationDirection dir=alongMomentum) | |
| SteppingHelixPropagator () | |
| Constructors. | |
| ~SteppingHelixPropagator () | |
| Destructor. | |
Protected Types | |
| typedef SteppingHelixStateInfo::VolumeBounds | MatBounds |
Protected Member Functions | |
| int | cIndex_ (int ind) const |
| (Internals) circular index for array of transients | |
| 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 | |
| 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 |
| bool | isYokeVolume (const MagVolume *vol) const |
| check if it's a yoke/iron based on this MagVol internals | |
| 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 | |
| Result | propagate (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 | |
| 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) const |
| (Internals) Init starting point | |
| void | setRep (SteppingHelixPropagator::Basis &rep, const SteppingHelixPropagator::Vector &tau) const |
| Set/compute basis vectors for local coordinates at current step (called by incrementState) | |
Private Types | |
| typedef std::pair < FreeTrajectoryState, double > | FtsPP |
| typedef std::pair < TrajectoryStateOnSurface, double > | TsosPP |
Private Attributes | |
| bool | applyRadX0Correction_ |
| AlgebraicMatrix55 | covCurvRot_ |
| AlgebraicMatrix55 | dCCurvTransform_ |
| bool | debug_ |
| double | defaultStep_ |
| double | ecShiftNeg_ |
| double | ecShiftPos_ |
| const MagneticField * | field_ |
| StateInfo | invalidState_ |
| bool | noErrorPropagation_ |
| bool | noMaterialMode_ |
| int | nPoints_ |
| bool | returnTangentPlane_ |
| bool | sendLogWarning_ |
| StateInfo | svBuf_ [MAX_POINTS+1] |
| const AlgebraicSymMatrix55 | unit55_ |
| bool | useInTeslaFromMagField_ |
| bool | useIsYokeFlag_ |
| bool | useMagVolumes_ |
| bool | useMatVolumes_ |
| bool | useTuningForL2Speed_ |
| const VolumeBasedMagneticField * | vbField_ |
Static Private Attributes | |
| static const int | MAX_POINTS = 7 |
| static const int | MAX_STEPS = 10000 |
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.
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.
Definition at line 44 of file SteppingHelixPropagator.h.
typedef std::pair<FreeTrajectoryState, double> SteppingHelixPropagator::FtsPP [private] |
Definition at line 268 of file SteppingHelixPropagator.h.
typedef SteppingHelixStateInfo::VolumeBounds SteppingHelixPropagator::MatBounds [protected] |
Definition at line 200 of file SteppingHelixPropagator.h.
| typedef CLHEP::Hep3Vector SteppingHelixPropagator::Point |
Definition at line 47 of file SteppingHelixPropagator.h.
Definition at line 50 of file SteppingHelixPropagator.h.
Definition at line 49 of file SteppingHelixPropagator.h.
typedef std::pair<TrajectoryStateOnSurface, double> SteppingHelixPropagator::TsosPP [private] |
Definition at line 267 of file SteppingHelixPropagator.h.
| typedef CLHEP::Hep3Vector SteppingHelixPropagator::Vector |
Definition at line 46 of file SteppingHelixPropagator.h.
| RADIUS_DT | |
| Z_DT | |
| PLANE_DT | |
| CONE_DT | |
| CYLINDER_DT | |
| PATHL_DT | |
| POINT_PCA_DT | |
| LINE_PCA_DT | |
| UNDEFINED_DT |
Definition at line 64 of file SteppingHelixPropagator.h.
{
RADIUS_DT=0,
Z_DT,
PLANE_DT,
CONE_DT,
CYLINDER_DT,
PATHL_DT,
POINT_PCA_DT,
LINE_PCA_DT,
UNDEFINED_DT
};
Definition at line 76 of file SteppingHelixPropagator.h.
Definition at line 58 of file SteppingHelixPropagator.h.
| SteppingHelixPropagator::SteppingHelixPropagator | ( | ) |
Constructors.
Definition at line 42 of file SteppingHelixPropagator.cc.
References field_.
Referenced by clone().
: Propagator(anyDirection) { field_ = 0; }
| SteppingHelixPropagator::SteppingHelixPropagator | ( | const MagneticField * | field, |
| PropagationDirection | dir = alongMomentum |
||
| ) |
Definition at line 48 of file SteppingHelixPropagator.cc.
References applyRadX0Correction_, SteppingHelixStateInfo::bf, SteppingHelixStateInfo::bfGradLoc, SteppingHelixStateInfo::covCurv, covCurvRot_, dCCurvTransform_, debug_, defaultStep_, ecShiftNeg_, ecShiftPos_, field_, SteppingHelixStateInfo::hasErrorPropagated_, i, SteppingHelixStateInfo::isComplete, SteppingHelixStateInfo::isValid_, SteppingHelixStateInfo::matDCovCurv, MAX_POINTS, noErrorPropagation_, noMaterialMode_, SteppingHelixStateInfo::p3, SteppingHelixStateInfo::r3, returnTangentPlane_, sendLogWarning_, svBuf_, unit55_, useInTeslaFromMagField_, useIsYokeFlag_, useMagVolumes_, useMatVolumes_, useTuningForL2Speed_, and vbField_.
: Propagator(dir), unit55_(AlgebraicMatrixID()) { field_ = field; vbField_ = dynamic_cast<const VolumeBasedMagneticField*>(field_); covCurvRot_ = AlgebraicMatrix55(); dCCurvTransform_ = unit55_; debug_ = false; noMaterialMode_ = false; noErrorPropagation_ = false; applyRadX0Correction_ = true; useMagVolumes_ = true; useIsYokeFlag_ = true; useMatVolumes_ = true; useInTeslaFromMagField_ = false; //overrides behavior only if true returnTangentPlane_ = true; sendLogWarning_ = false; useTuningForL2Speed_ = false; for (int i = 0; i <= MAX_POINTS; i++){ svBuf_[i].p3 = Vector(0,0,0); svBuf_[i].r3 = Point(0,0,0); svBuf_[i].bf = Vector(0,0,0); svBuf_[i].bfGradLoc = Vector(0,0,0); svBuf_[i].covCurv = AlgebraicSymMatrix55(); svBuf_[i].matDCovCurv = AlgebraicSymMatrix55(); svBuf_[i].isComplete = true; svBuf_[i].isValid_ = true; svBuf_[i].hasErrorPropagated_ = !noErrorPropagation_; } defaultStep_ = 5.; ecShiftPos_ = 0; ecShiftNeg_ = 0; }
| SteppingHelixPropagator::~SteppingHelixPropagator | ( | ) | [inline] |
| 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 164 of file SteppingHelixPropagator.h.
References applyRadX0Correction(), and applyRadX0Correction_.
Referenced by applyRadX0Correction(), TrackDetectorAssociator::init(), HTrackAssociator::init(), SteppingHelixPropagatorESProducer::produce(), and AlCaHOCalibProducer::produce().
| int SteppingHelixPropagator::cIndex_ | ( | int | ind | ) | const [protected] |
(Internals) circular index for array of transients
Definition at line 1682 of file SteppingHelixPropagator.cc.
References MAX_POINTS, and query::result.
Referenced by propagate(), and setIState().
{
int result = ind%MAX_POINTS;
if (ind != 0 && result == 0){
result = MAX_POINTS;
}
return result;
}
| virtual SteppingHelixPropagator* SteppingHelixPropagator::clone | ( | void | ) | const [inline, virtual] |
Implements Propagator.
Definition at line 88 of file SteppingHelixPropagator.h.
References SteppingHelixPropagator().
{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 1322 of file SteppingHelixPropagator.cc.
References SteppingHelixStateInfo::bf, ecShiftNeg_, ecShiftPos_, SteppingHelixStateInfo::isYokeVol, funct::log(), SteppingHelixStateInfo::magVol, noMaterialMode_, SteppingHelixStateInfo::p3, SteppingHelixStateInfo::r3, mathSSE::sqrt(), useIsYokeFlag_, and useMagVolumes_.
Referenced by getNextState(), loadState(), and makeAtomStep().
{
radX0 = 1.e24;
dEdXPrime = 0.;
rzLims = MatBounds();
if (noMaterialMode_) return 0;
double dEdx = 0.;
double lR = sv.r3.perp();
double lZ = fabs(sv.r3.z());
//assume "Iron" .. seems to be quite the same for brass/iron/PbW04
//good for Fe within 3% for 0.2 GeV to 10PeV
double dEdX_HCal = 0.95; //extracted from sim
double dEdX_ECal = 0.45;
double dEdX_coil = 0.35; //extracted from sim .. closer to 40% in fact
double dEdX_Fe = 1;
double dEdX_MCh = 0.053; //chambers on average
double dEdX_Trk = 0.0114;
double dEdX_Air = 2E-4;
double dEdX_Vac = 0.0;
double radX0_HCal = 1.44/0.8; //guessing
double radX0_ECal = 0.89/0.7;
double radX0_coil = 4.; //
double radX0_Fe = 1.76;
double radX0_MCh = 1e3; //
double radX0_Trk = 320.;
double radX0_Air = 3.e4;
double radX0_Vac = 3.e9; //"big" number for vacuum
//not all the boundaries are set below: this will be a default
if (! (lR < 380 && lZ < 785)){
if (lZ > 785 ) rzLims = MatBounds(0, 1e4, 785, 1e4);
if (lZ < 785 ) rzLims = MatBounds(380, 1e4, 0, 785);
}
//this def makes sense assuming we don't jump into endcap volume from the other z-side in one step
//also, it is a positive shift only (at least for now): can't move ec further into the detector
double ecShift = sv.r3.z() > 0 ? fabs(ecShiftPos_) : fabs(ecShiftNeg_);
//this should roughly figure out where things are
//(numbers taken from Fig1.1.2 TDR and from geom xmls)
if (lR < 2.9){ //inside beampipe
dEdx = dEdX_Vac; radX0 = radX0_Vac;
rzLims = MatBounds(0, 2.9, 0, 1E4);
}
else if (lR < 129){
if (lZ < 294){
rzLims = MatBounds(2.9, 129, 0, 294); //tracker boundaries
dEdx = dEdX_Trk; radX0 = radX0_Trk;
//somewhat empirical formula that ~ matches the average if going from 0,0,0
//assuming "uniform" tracker material
//doesn't really track material layer to layer
double lEtaDet = fabs(sv.r3.eta());
double scaleRadX = lEtaDet > 1.5 ? 0.7724 : 1./cosh(0.5*lEtaDet);//sin(2.*atan(exp(-0.5*lEtaDet)));
scaleRadX *= scaleRadX;
if (lEtaDet > 2 && lZ > 20) scaleRadX *= (lEtaDet-1.);
if (lEtaDet > 2.5 && lZ > 20) scaleRadX *= (lEtaDet-1.);
radX0 *= scaleRadX;
}
//endcap part begins here
else if ( lZ < 294 + ecShift ){
//gap in front of EE here, piece inside 2.9<R<129
rzLims = MatBounds(2.9, 129, 294, 294 + ecShift);
dEdx = dEdX_Air; radX0 = radX0_Air;
}
else if (lZ < 372 + ecShift){
rzLims = MatBounds(2.9, 129, 294 + ecShift, 372 + ecShift); //EE here, piece inside 2.9<R<129
dEdx = dEdX_ECal; radX0 = radX0_ECal;
}//EE averaged out over a larger space
else if (lZ < 398 + ecShift){
rzLims = MatBounds(2.9, 129, 372 + ecShift, 398 + ecShift); //whatever goes behind EE 2.9<R<129 is air up to Z=398
dEdx = dEdX_HCal*0.05; radX0 = radX0_Air;
}//betw EE and HE
else if (lZ < 555 + ecShift){
rzLims = MatBounds(2.9, 129, 398 + ecShift, 555 + ecShift); //HE piece 2.9<R<129;
dEdx = dEdX_HCal*0.96; radX0 = radX0_HCal/0.96;
} //HE calor abit less dense
else {
// rzLims = MatBounds(2.9, 129, 555, 785);
// set the boundaries first: they serve as stop-points here
// the material is set below
if (lZ < 568 + ecShift) rzLims = MatBounds(2.9, 129, 555 + ecShift, 568 + ecShift); //a piece of HE support R<129, 555<Z<568
else if (lZ < 625 + ecShift){
if (lR < 85 + ecShift) rzLims = MatBounds(2.9, 85, 568 + ecShift, 625 + ecShift);
else rzLims = MatBounds(85, 129, 568 + ecShift, 625 + ecShift);
} else if (lZ < 785 + ecShift) rzLims = MatBounds(2.9, 129, 625 + ecShift, 785 + ecShift);
else if (lZ < 1500 + ecShift) rzLims = MatBounds(2.9, 129, 785 + ecShift, 1500 + ecShift);
else rzLims = MatBounds(2.9, 129, 1500 + ecShift, 1E4);
//iron .. don't care about no material in front of HF (too forward)
if (! (lZ > 568 + ecShift && lZ < 625 + ecShift && lR > 85 ) // HE support
&& ! (lZ > 785 + ecShift && lZ < 850 + ecShift && lR > 118)) {dEdx = dEdX_Fe; radX0 = radX0_Fe; }
else { dEdx = dEdX_MCh; radX0 = radX0_MCh; } //ME at eta > 2.2
}
}
else if (lR < 287){
if (lZ < 372 + ecShift && lR < 177){ // 129<<R<177
if (lZ < 304) rzLims = MatBounds(129, 177, 0, 304); //EB 129<<R<177 0<Z<304
else if (lZ < 343){ // 129<<R<177 304<Z<343
if (lR < 135 ) rzLims = MatBounds(129, 135, 304, 343);// tk piece 129<<R<135 304<Z<343
else if (lR < 172 ){ //
if (lZ < 311 ) rzLims = MatBounds(135, 172, 304, 311);
else rzLims = MatBounds(135, 172, 311, 343);
} else {
if (lZ < 328) rzLims = MatBounds(172, 177, 304, 328);
else rzLims = MatBounds(172, 177, 328, 343);
}
}
else if ( lZ < 343 + ecShift){
rzLims = MatBounds(129, 177, 343, 343 + ecShift); //gap
}
else {
if (lR < 156 ) rzLims = MatBounds(129, 156, 343 + ecShift, 372 + ecShift);
else if ( (lZ - 343 - ecShift) > (lR - 156)*1.38 )
rzLims = MatBounds(156, 177, 127.72 + ecShift, 372 + ecShift, atan(1.38), 0);
else rzLims = MatBounds(156, 177, 343 + ecShift, 127.72 + ecShift, 0, atan(1.38));
}
if (!(lR > 135 && lZ <343 + ecShift && lZ > 304 )
&& ! (lR > 156 && lZ < 372 + ecShift && lZ > 343 + ecShift && ((lZ-343. - ecShift )< (lR-156.)*1.38)))
{
//the crystals are the same length, but they are not 100% of material
double cosThetaEquiv = 0.8/sqrt(1.+lZ*lZ/lR/lR) + 0.2;
if (lZ > 343) cosThetaEquiv = 1.;
dEdx = dEdX_ECal*cosThetaEquiv; radX0 = radX0_ECal/cosThetaEquiv;
} //EB
else {
if ( (lZ > 304 && lZ < 328 && lR < 177 && lR > 135)
&& ! (lZ > 311 && lR < 172) ) {dEdx = dEdX_Fe; radX0 = radX0_Fe; } //Tk_Support
else if ( lZ > 343 && lZ < 343 + ecShift) { dEdx = dEdX_Air; radX0 = radX0_Air; }
else {dEdx = dEdX_ECal*0.2; radX0 = radX0_Air;} //cables go here <-- will be abit too dense for ecShift > 0
}
}
else if (lZ < 554 + ecShift){ // 129<R<177 372<Z<554 AND 177<R<287 0<Z<554
if (lR < 177){ // 129<R<177 372<Z<554
if ( lZ > 372 + ecShift && lZ < 398 + ecShift )rzLims = MatBounds(129, 177, 372 + ecShift, 398 + ecShift); // EE 129<R<177 372<Z<398
else if (lZ < 548 + ecShift) rzLims = MatBounds(129, 177, 398 + ecShift, 548 + ecShift); // HE 129<R<177 398<Z<548
else rzLims = MatBounds(129, 177, 548 + ecShift, 554 + ecShift); // HE gap 129<R<177 548<Z<554
}
else if (lR < 193){ // 177<R<193 0<Z<554
if ((lZ - 307) < (lR - 177.)*1.739) rzLims = MatBounds(177, 193, 0, -0.803, 0, atan(1.739));
else if ( lZ < 389) rzLims = MatBounds(177, 193, -0.803, 389, atan(1.739), 0.);
else if ( lZ < 389 + ecShift) rzLims = MatBounds(177, 193, 389, 389 + ecShift); // air insert
else if ( lZ < 548 + ecShift ) rzLims = MatBounds(177, 193, 389 + ecShift, 548 + ecShift);// HE 177<R<193 389<Z<548
else rzLims = MatBounds(177, 193, 548 + ecShift, 554 + ecShift); // HE gap 177<R<193 548<Z<554
}
else if (lR < 264){ // 193<R<264 0<Z<554
double anApex = 375.7278 - 193./1.327; // 230.28695599096
if ( (lZ - 375.7278) < (lR - 193.)/1.327) rzLims = MatBounds(193, 264, 0, anApex, 0, atan(1./1.327)); //HB
else if ( (lZ - 392.7278 ) < (lR - 193.)/1.327)
rzLims = MatBounds(193, 264, anApex, anApex+17., atan(1./1.327), atan(1./1.327)); // HB-HE gap
else if ( (lZ - 392.7278 - ecShift ) < (lR - 193.)/1.327)
rzLims = MatBounds(193, 264, anApex+17., anApex+17. + ecShift, atan(1./1.327), atan(1./1.327)); // HB-HE gap air insert
// HE (372,193)-(517,193)-(517,264)-(417.8,264)
else if ( lZ < 517 + ecShift ) rzLims = MatBounds(193, 264, anApex+17. + ecShift, 517 + ecShift, atan(1./1.327), 0);
else if (lZ < 548 + ecShift){ // HE+gap 193<R<264 517<Z<548
if (lR < 246 ) rzLims = MatBounds(193, 246, 517 + ecShift, 548 + ecShift); // HE 193<R<246 517<Z<548
else rzLims = MatBounds(246, 264, 517 + ecShift, 548 + ecShift); // HE gap 246<R<264 517<Z<548
}
else rzLims = MatBounds(193, 264, 548 + ecShift, 554 + ecShift); // HE gap 193<R<246 548<Z<554
}
else if ( lR < 275){ // 264<R<275 0<Z<554
if (lZ < 433) rzLims = MatBounds(264, 275, 0, 433); //HB
else if (lZ < 554 ) rzLims = MatBounds(264, 275, 433, 554); // HE gap 264<R<275 433<Z<554
else rzLims = MatBounds(264, 275, 554, 554 + ecShift); // HE gap 264<R<275 554<Z<554 air insert
}
else { // 275<R<287 0<Z<554
if (lZ < 402) rzLims = MatBounds(275, 287, 0, 402);// HB 275<R<287 0<Z<402
else if (lZ < 554) rzLims = MatBounds(275, 287, 402, 554);// //HE gap 275<R<287 402<Z<554
else rzLims = MatBounds(275, 287, 554, 554 + ecShift); //HE gap 275<R<287 554<Z<554 air insert
}
if ((lZ < 433 || lR < 264) && (lZ < 402 || lR < 275) && (lZ < 517 + ecShift || lR < 246) //notches
//I should've made HE and HF different .. now need to shorten HE to match
&& lZ < 548 + ecShift
&& ! (lZ < 389 + ecShift && lZ > 335 && lR < 193 ) //not a gap EB-EE 129<R<193
&& ! (lZ > 307 && lZ < 335 && lR < 193 && ((lZ - 307) > (lR - 177.)*1.739)) //not a gap
&& ! (lR < 177 && lZ < 398 + ecShift) //under the HE nose
&& ! (lR < 264 && lR > 193 && fabs(441.5+0.5*ecShift - lZ + (lR - 269.)/1.327) < (8.5 + ecShift*0.5)) ) //not a gap
{ dEdx = dEdX_HCal; radX0 = radX0_HCal; //hcal
}
else if( (lR < 193 && lZ > 389 && lZ < 389 + ecShift ) || (lR > 264 && lR < 287 && lZ > 554 && lZ < 554 + ecShift)
||(lR < 264 && lR > 193 && fabs(441.5+8.5+0.5*ecShift - lZ + (lR - 269.)/1.327) < ecShift*0.5) ) {
dEdx = dEdX_Air; radX0 = radX0_Air;
}
else {dEdx = dEdX_HCal*0.05; radX0 = radX0_Air; }//endcap gap
}
//the rest is a tube of endcap volume 129<R<251 554<Z<largeValue
else if (lZ < 564 + ecShift){ // 129<R<287 554<Z<564
if (lR < 251) {
rzLims = MatBounds(129, 251, 554 + ecShift, 564 + ecShift); // HE support 129<R<251 554<Z<564
dEdx = dEdX_Fe; radX0 = radX0_Fe;
}//HE support
else {
rzLims = MatBounds(251, 287, 554 + ecShift, 564 + ecShift); //HE/ME gap 251<R<287 554<Z<564
dEdx = dEdX_MCh; radX0 = radX0_MCh;
}
}
else if (lZ < 625 + ecShift){ // ME/1/1 129<R<287 564<Z<625
rzLims = MatBounds(129, 287, 564 + ecShift, 625 + ecShift);
dEdx = dEdX_MCh; radX0 = radX0_MCh;
}
else if (lZ < 785 + ecShift){ //129<R<287 625<Z<785
if (lR < 275) rzLims = MatBounds(129, 275, 625 + ecShift, 785 + ecShift); //YE/1 129<R<275 625<Z<785
else { // 275<R<287 625<Z<785
if (lZ < 720 + ecShift) rzLims = MatBounds(275, 287, 625 + ecShift, 720 + ecShift); // ME/1/2 275<R<287 625<Z<720
else rzLims = MatBounds(275, 287, 720 + ecShift, 785 + ecShift); // YE/1 275<R<287 720<Z<785
}
if (! (lR > 275 && lZ < 720 + ecShift)) { dEdx = dEdX_Fe; radX0 = radX0_Fe; }//iron
else { dEdx = dEdX_MCh; radX0 = radX0_MCh; }
}
else if (lZ < 850 + ecShift){
rzLims = MatBounds(129, 287, 785 + ecShift, 850 + ecShift);
dEdx = dEdX_MCh; radX0 = radX0_MCh;
}
else if (lZ < 910 + ecShift){
rzLims = MatBounds(129, 287, 850 + ecShift, 910 + ecShift);
dEdx = dEdX_Fe; radX0 = radX0_Fe;
}//iron
else if (lZ < 975 + ecShift){
rzLims = MatBounds(129, 287, 910 + ecShift, 975 + ecShift);
dEdx = dEdX_MCh; radX0 = radX0_MCh;
}
else if (lZ < 1000 + ecShift){
rzLims = MatBounds(129, 287, 975 + ecShift, 1000 + ecShift);
dEdx = dEdX_Fe; radX0 = radX0_Fe;
}//iron
else if (lZ < 1063 + ecShift){
rzLims = MatBounds(129, 287, 1000 + ecShift, 1063 + ecShift);
dEdx = dEdX_MCh; radX0 = radX0_MCh;
}
else if ( lZ < 1073 + ecShift){
rzLims = MatBounds(129, 287, 1063 + ecShift, 1073 + ecShift);
dEdx = dEdX_Fe; radX0 = radX0_Fe;
}
else if (lZ < 1E4 ) {
rzLims = MatBounds(129, 287, 1073 + ecShift, 1E4);
dEdx = dEdX_Air; radX0 = radX0_Air;
}
else {
dEdx = dEdX_Air; radX0 = radX0_Air;
}
}
else if (lR <380 && lZ < 667){
if (lZ < 630){
if (lR < 315) rzLims = MatBounds(287, 315, 0, 630);
else if (lR < 341 ) rzLims = MatBounds(315, 341, 0, 630); //b-field ~linear rapid fall here
else rzLims = MatBounds(341, 380, 0, 630);
} else rzLims = MatBounds(287, 380, 630, 667);
if (lZ < 630) { dEdx = dEdX_coil; radX0 = radX0_coil; }//a guess for the solenoid average
else {dEdx = dEdX_Air; radX0 = radX0_Air; }//endcap gap
}
else {
if (lZ < 667) {
if (lR < 850){
bool isIron = false;
//sanity check in addition to flags
if (useIsYokeFlag_ && useMagVolumes_ && sv.magVol != 0){
isIron = sv.isYokeVol;
} else {
double bMag = sv.bf.mag();
isIron = (bMag > 0.75 && ! (lZ > 500 && lR <500 && bMag < 1.15)
&& ! (lZ < 450 && lR > 420 && bMag < 1.15 ) );
}
//tell the material stepper where mat bounds are
rzLims = MatBounds(380, 850, 0, 667);
if (isIron) { dEdx = dEdX_Fe; radX0 = radX0_Fe; }//iron
else { dEdx = dEdX_MCh; radX0 = radX0_MCh; }
} else {
rzLims = MatBounds(850, 1E4, 0, 667);
dEdx = dEdX_Air; radX0 = radX0_Air;
}
}
else if (lR > 750 ){
rzLims = MatBounds(750, 1E4, 667, 1E4);
dEdx = dEdX_Air; radX0 = radX0_Air;
}
else if (lZ < 667 + ecShift){
rzLims = MatBounds(287, 750, 667, 667 + ecShift);
dEdx = dEdX_Air; radX0 = radX0_Air;
}
//the rest is endcap piece with 287<R<750 Z>667
else if (lZ < 724 + ecShift){
if (lR < 380 ) rzLims = MatBounds(287, 380, 667 + ecShift, 724 + ecShift);
else rzLims = MatBounds(380, 750, 667 + ecShift, 724 + ecShift);
dEdx = dEdX_MCh; radX0 = radX0_MCh;
}
else if (lZ < 785 + ecShift){
if (lR < 380 ) rzLims = MatBounds(287, 380, 724 + ecShift, 785 + ecShift);
else rzLims = MatBounds(380, 750, 724 + ecShift, 785 + ecShift);
dEdx = dEdX_Fe; radX0 = radX0_Fe;
}//iron
else if (lZ < 850 + ecShift){
rzLims = MatBounds(287, 750, 785 + ecShift, 850 + ecShift);
dEdx = dEdX_MCh; radX0 = radX0_MCh;
}
else if (lZ < 910 + ecShift){
rzLims = MatBounds(287, 750, 850 + ecShift, 910 + ecShift);
dEdx = dEdX_Fe; radX0 = radX0_Fe;
}//iron
else if (lZ < 975 + ecShift){
rzLims = MatBounds(287, 750, 910 + ecShift, 975 + ecShift);
dEdx = dEdX_MCh; radX0 = radX0_MCh;
}
else if (lZ < 1000 + ecShift){
rzLims = MatBounds(287, 750, 975 + ecShift, 1000 + ecShift);
dEdx = dEdX_Fe; radX0 = radX0_Fe;
}//iron
else if (lZ < 1063 + ecShift){
if (lR < 360){
rzLims = MatBounds(287, 360, 1000 + ecShift, 1063 + ecShift);
dEdx = dEdX_MCh; radX0 = radX0_MCh;
}
//put empty air where me4/2 should be (future)
else {
rzLims = MatBounds(360, 750, 1000 + ecShift, 1063 + ecShift);
dEdx = dEdX_Air; radX0 = radX0_Air;
}
}
else if ( lZ < 1073 + ecShift){
rzLims = MatBounds(287, 750, 1063 + ecShift, 1073 + ecShift);
//this plate does not exist: air
dEdx = dEdX_Air; radX0 = radX0_Air;
}
else if (lZ < 1E4 ) {
rzLims = MatBounds(287, 750, 1073 + ecShift, 1E4);
dEdx = dEdX_Air; radX0 = radX0_Air;
}
else {dEdx = dEdX_Air; radX0 = radX0_Air; }//air
}
//dEdx so far is a relative number (relative to iron)
//scale by what's expected for iron (the function was fit from pdg table)
//0.065 (PDG) --> 0.044 to better match with MPV
double p0 = sv.p3.mag();
double logp0 = log(p0);
double p0powN33 = 0;
if (p0>3.) {
// p0powN33 = exp(-0.33*logp0); //calculate for p>3GeV
double xx=1./p0; xx=sqrt(sqrt(sqrt(sqrt(xx)))); xx=xx*xx*xx*xx*xx; // this is (p0)**(-5/16), close enough to -0.33
p0powN33 = xx;
}
double dEdX_mat = -(11.4 + 0.96*fabs(logp0+log(2.8)) + 0.033*p0*(1.0 - p0powN33) )*1e-3;
//in GeV/cm .. 0.8 to get closer to the median or MPV
dEdXPrime = dEdx == 0 ? 0 : -dEdx*(0.96/p0 + 0.033 - 0.022*p0powN33)*1e-3; //== d(dEdX)/dp
dEdx *= dEdX_mat;
return dEdx;
}
| 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 780 of file SteppingHelixPropagator.cc.
References SteppingHelixStateInfo::bf, SteppingHelixStateInfo::covCurv, debug_, SteppingHelixStateInfo::dEdx, SteppingHelixStateInfo::dEdXPrime, SteppingHelixStateInfo::dir, field_, VolumeBasedMagneticField::findVolume(), getDeDx(), SteppingHelixStateInfo::hasErrorPropagated_, MagVolume::inTesla(), MagneticField::inTesla(), SteppingHelixStateInfo::isYokeVol, isYokeVolume(), VolumeBasedMagneticField::isZSymmetric(), LogTrace, SteppingHelixStateInfo::magVol, SteppingHelixStateInfo::matDCovCurv, metname, SteppingHelixStateInfo::p3, SteppingHelixStateInfo::path(), SteppingHelixStateInfo::path_, SteppingHelixStateInfo::q, SteppingHelixStateInfo::r3, SteppingHelixStateInfo::radPath(), SteppingHelixStateInfo::radPath_, SteppingHelixStateInfo::radX0, SteppingHelixStateInfo::rzLims, metsig::tau, tmp, useInTeslaFromMagField_, useIsYokeFlag_, useMagVolumes_, and vbField_.
Referenced by makeAtomStep().
{
static const std::string metname = "SteppingHelixPropagator";
svNext.q = svPrevious.q;
svNext.dir = dS > 0.0 ? 1.: -1.;
svNext.p3 = tau; svNext.p3*=(svPrevious.p3.mag() - svNext.dir*fabs(dP));
svNext.r3 = svPrevious.r3; svNext.r3 += drVec;
svNext.path_ = svPrevious.path() + dS;
svNext.radPath_ = svPrevious.radPath() + dX0;
GlobalPoint gPointNegZ(svNext.r3.x(), svNext.r3.y(), -fabs(svNext.r3.z()));
GlobalPoint gPointNorZ(svNext.r3.x(), svNext.r3.y(), svNext.r3.z());
GlobalVector bf(0,0,0);
if (useMagVolumes_){
if (vbField_ != 0){
if (vbField_->isZSymmetric()){
svNext.magVol = vbField_->findVolume(gPointNegZ);
} else {
svNext.magVol = vbField_->findVolume(gPointNorZ);
}
if (useIsYokeFlag_){
double curRad = svNext.r3.perp();
if (curRad > 380 && curRad < 850 && fabs(svNext.r3.z()) < 667){
svNext.isYokeVol = isYokeVolume(svNext.magVol);
} else {
svNext.isYokeVol = false;
}
}
} else {
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Failed to cast into VolumeBasedMagneticField"<<std::endl;
svNext.magVol = 0;
}
if (debug_){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Got volume at "<<svNext.magVol<<std::endl;
}
}
if (useMagVolumes_ && svNext.magVol != 0 && ! useInTeslaFromMagField_){
if (vbField_->isZSymmetric()){
bf = svNext.magVol->inTesla(gPointNegZ);
} else {
bf = svNext.magVol->inTesla(gPointNorZ);
}
if (svNext.r3.z() > 0 && vbField_->isZSymmetric() ){
svNext.bf.set(-bf.x(), -bf.y(), bf.z());
} else {
svNext.bf.set(bf.x(), bf.y(), bf.z());
}
} else {
GlobalPoint gPoint(svNext.r3.x(), svNext.r3.y(), svNext.r3.z());
bf = field_->inTesla(gPoint);
svNext.bf.set(bf.x(), bf.y(), bf.z());
}
if (svNext.bf.mag2() < 1e-32) svNext.bf.set(0., 0., 1e-16);
double dEdXPrime = 0;
double dEdx = 0;
double radX0 = 0;
MatBounds rzLims;
dEdx = getDeDx(svNext, dEdXPrime, radX0, rzLims);
svNext.dEdx = dEdx; svNext.dEdXPrime = dEdXPrime;
svNext.radX0 = radX0;
svNext.rzLims = rzLims;
//update Emat only if it's valid
svNext.hasErrorPropagated_ = svPrevious.hasErrorPropagated_;
if (svPrevious.hasErrorPropagated_){
{
AlgebraicMatrix55 tmp = dCovCurvTransform*svPrevious.covCurv;
ROOT::Math::AssignSym::Evaluate(svNext.covCurv, tmp*ROOT::Math::Transpose(dCovCurvTransform));
svNext.covCurv += svPrevious.matDCovCurv;
}
} else {
//could skip dragging along the unprop. cov later .. now
// svNext.cov = svPrevious.cov;
}
if (debug_){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Now at path: "<<svNext.path()<<" radPath: "<<svNext.radPath()
<<" p3 "<<" pt: "<<svNext.p3.perp()<<" phi: "<<svNext.p3.phi()
<<" eta: "<<svNext.p3.eta()
<<" "<<svNext.p3
<<" r3: "<<svNext.r3
<<" dPhi: "<<acos(svNext.p3.unit().dot(svPrevious.p3.unit()))
<<" bField: "<<svNext.bf.mag()
<<" dedx: "<<svNext.dEdx
<<std::endl;
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"New Covariance "<<svNext.covCurv<<std::endl;
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Transf by dCovTransform "<<dCovCurvTransform<<std::endl;
}
}
| bool SteppingHelixPropagator::isYokeVolume | ( | const MagVolume * | vol | ) | const [protected] |
check if it's a yoke/iron based on this MagVol internals
Definition at line 2318 of file SteppingHelixPropagator.cc.
References debug_, MagVolume::isIron(), LogTrace, metname, GloballyPositioned< T >::position(), and query::result.
Referenced by getNextState(), and loadState().
{
static const std::string metname = "SteppingHelixPropagator";
if (vol == 0) return false;
/*
const MFGrid* mGrid = reinterpret_cast<const MFGrid*>(vol->provider());
std::vector<int> dims(mGrid->dimensions());
LocalVector lVCen(mGrid->nodeValue(dims[0]/2, dims[1]/2, dims[2]/2));
LocalVector lVZLeft(mGrid->nodeValue(dims[0]/2, dims[1]/2, dims[2]/5));
LocalVector lVZRight(mGrid->nodeValue(dims[0]/2, dims[1]/2, (dims[2]*4)/5));
double mag2VCen = lVCen.mag2();
double mag2VZLeft = lVZLeft.mag2();
double mag2VZRight = lVZRight.mag2();
bool result = false;
if (mag2VCen > 0.6 && mag2VZLeft > 0.6 && mag2VZRight > 0.6){
if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Volume is magnetic, located at "<<vol->position()<<std::endl;
result = true;
} else {
if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Volume is not magnetic, located at "<<vol->position()<<std::endl;
}
*/
bool result = vol->isIron();
if (result){
if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Volume is magnetic, located at "<<vol->position()<<std::endl;
} else {
if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Volume is not magnetic, located at "<<vol->position()<<std::endl;
}
return result;
}
| 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 651 of file SteppingHelixPropagator.cc.
References alongMomentum, SteppingHelixStateInfo::bf, DeDxDiscriminatorTools::charge(), SteppingHelixStateInfo::covCurv, debug_, SteppingHelixStateInfo::dEdx, SteppingHelixStateInfo::dEdXPrime, SteppingHelixStateInfo::dir, field_, MagVolume::fieldInTesla(), VolumeBasedMagneticField::findVolume(), getDeDx(), MagVolume::inTesla(), MagneticField::inTesla(), SteppingHelixStateInfo::isComplete, SteppingHelixStateInfo::isValid_, SteppingHelixStateInfo::isYokeVol, isYokeVolume(), VolumeBasedMagneticField::isZSymmetric(), LogTrace, PV3DBase< T, PVType, FrameType >::mag2(), SteppingHelixStateInfo::magVol, metname, p3, SteppingHelixStateInfo::p3, SteppingHelixStateInfo::path(), SteppingHelixStateInfo::path_, SteppingHelixStateInfo::q, SteppingHelixStateInfo::r3, SteppingHelixStateInfo::radPath(), SteppingHelixStateInfo::radPath_, SteppingHelixStateInfo::radX0, SteppingHelixStateInfo::rzLims, 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().
{
static const std::string metname = "SteppingHelixPropagator";
svCurrent.q = charge;
svCurrent.p3 = p3;
svCurrent.r3 = r3;
svCurrent.dir = dir == alongMomentum ? 1.: -1.;
svCurrent.path_ = 0; // this could've held the initial path
svCurrent.radPath_ = 0;
GlobalPoint gPointNegZ(svCurrent.r3.x(), svCurrent.r3.y(), -fabs(svCurrent.r3.z()));
GlobalPoint gPointNorZ(svCurrent.r3.x(), svCurrent.r3.y(), svCurrent.r3.z());
float gpmag = gPointNegZ.mag2();
float pmag2 = p3.mag2();
if (gpmag > 1e20f ) {
LogTrace(metname)<<"Initial point is too far";
svCurrent.isValid_ = false;
return;
}
if (! (gpmag == gpmag) ) {
LogTrace(metname)<<"Initial point is a nan";
edm::LogWarning("SteppingHelixPropagatorNAN")<<"Initial point is a nan";
svCurrent.isValid_ = false;
return;
}
if (! (pmag2 == pmag2) ) {
LogTrace(metname)<<"Initial momentum is a nan";
edm::LogWarning("SteppingHelixPropagatorNAN")<<"Initial momentum is a nan";
svCurrent.isValid_ = false;
return;
}
GlobalVector bf(0,0,0);
// = field_->inTesla(gPoint);
if (useMagVolumes_){
if (vbField_ ){
if (vbField_->isZSymmetric()){
svCurrent.magVol = vbField_->findVolume(gPointNegZ);
} else {
svCurrent.magVol = vbField_->findVolume(gPointNorZ);
}
if (useIsYokeFlag_){
double curRad = svCurrent.r3.perp();
if (curRad > 380 && curRad < 850 && fabs(svCurrent.r3.z()) < 667){
svCurrent.isYokeVol = isYokeVolume(svCurrent.magVol);
} else {
svCurrent.isYokeVol = false;
}
}
} else {
edm::LogWarning(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Failed to cast into VolumeBasedMagneticField: fall back to the default behavior"<<std::endl;
svCurrent.magVol = 0;
}
if (debug_){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Got volume at "<<svCurrent.magVol<<std::endl;
}
}
if (useMagVolumes_ && svCurrent.magVol != 0 && ! useInTeslaFromMagField_){
if (vbField_->isZSymmetric()){
bf = svCurrent.magVol->inTesla(gPointNegZ);
if (debug_){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Loaded bfield float: "<<bf
<<" at global float "<< gPointNegZ<<" double "<< svCurrent.r3<<std::endl;
LocalPoint lPoint(svCurrent.magVol->toLocal(gPointNegZ));
LocalVector lbf = svCurrent.magVol->fieldInTesla(lPoint);
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"\t cf in local locF: "<<lbf<<" at "<<lPoint<<std::endl;
}
} else {
bf = svCurrent.magVol->inTesla(gPointNorZ);
if (debug_){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Loaded bfield float: "<<bf
<<" at global float "<< gPointNorZ<<" double "<< svCurrent.r3<<std::endl;
LocalPoint lPoint(svCurrent.magVol->toLocal(gPointNorZ));
LocalVector lbf = svCurrent.magVol->fieldInTesla(lPoint);
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"\t cf in local locF: "<<lbf<<" at "<<lPoint<<std::endl;
}
}
if (r3.z() > 0 && vbField_->isZSymmetric() ){
svCurrent.bf.set(-bf.x(), -bf.y(), bf.z());
} else {
svCurrent.bf.set(bf.x(), bf.y(), bf.z());
}
} else {
GlobalPoint gPoint(r3.x(), r3.y(), r3.z());
bf = field_->inTesla(gPoint);
svCurrent.bf.set(bf.x(), bf.y(), bf.z());
}
if (svCurrent.bf.mag2() < 1e-32) svCurrent.bf.set(0., 0., 1e-16);
if (debug_){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific
<<"Loaded bfield double: "<<svCurrent.bf<<" from float: "<<bf
<<" at float "<< gPointNorZ<<" "<<gPointNegZ<<" double "<< svCurrent.r3<<std::endl;
}
double dEdXPrime = 0;
double dEdx = 0;
double radX0 = 0;
MatBounds rzLims;
dEdx = getDeDx(svCurrent, dEdXPrime, radX0, rzLims);
svCurrent.dEdx = dEdx; svCurrent.dEdXPrime = dEdXPrime;
svCurrent.radX0 = radX0;
svCurrent.rzLims = rzLims;
svCurrent.covCurv =covCurv;
svCurrent.isComplete = true;
svCurrent.isValid_ = true;
if (debug_){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Loaded at path: "<<svCurrent.path()<<" radPath: "<<svCurrent.radPath()
<<" p3 "<<" pt: "<<svCurrent.p3.perp()<<" phi: "<<svCurrent.p3.phi()
<<" eta: "<<svCurrent.p3.eta()
<<" "<<svCurrent.p3
<<" r3: "<<svCurrent.r3
<<" bField: "<<svCurrent.bf.mag()
<<std::endl;
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Input Covariance in Curvilinear RF "<<covCurv<<std::endl;
}
}
| virtual const MagneticField* SteppingHelixPropagator::magneticField | ( | ) | const [inline, virtual] |
Implements Propagator.
Definition at line 93 of file SteppingHelixPropagator.h.
References field_.
{ return field_;}
| 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 889 of file SteppingHelixPropagator.cc.
References alongMomentum, applyRadX0Correction_, SteppingHelixStateInfo::bf, funct::cos(), prof2calltree::cost, dCCurvTransform_, debug_, SteppingHelixStateInfo::dEdx, SteppingHelixStateInfo::dEdXPrime, field_, getDeDx(), getNextState(), SteppingHelixStateInfo::hasErrorPropagated_, HEL_ALL_F, HEL_AS_F, MagVolume::inTesla(), MagneticField::inTesla(), SteppingHelixStateInfo::isYokeVol, VolumeBasedMagneticField::isZSymmetric(), funct::log(), LogTrace, SteppingHelixPropagator::Basis::lX, SteppingHelixPropagator::Basis::lY, SteppingHelixPropagator::Basis::lZ, PV3DBase< T, PVType, FrameType >::mag(), SteppingHelixStateInfo::magVol, SteppingHelixStateInfo::matDCovCurv, metname, oppositeToMomentum, SteppingHelixStateInfo::p3, SteppingHelixStateInfo::path(), phi, createTree::pp, SteppingHelixStateInfo::q, lumiQueryAPI::q, SteppingHelixStateInfo::r3, SteppingHelixStateInfo::radPath(), SteppingHelixStateInfo::radX0, setRep(), funct::sin(), mathSSE::sqrt(), metsig::tau, unit55_, useInTeslaFromMagField_, useMagVolumes_, vbField_, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().
Referenced by propagate().
{
static const std::string metname = "SteppingHelixPropagator";
if (debug_){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Make atom step "<<svCurrent.path()<<" with step "<<dS<<" in direction "<<dir<<std::endl;
}
double dP = 0;
double curP = svCurrent.p3.mag();
Vector tau = svCurrent.p3; tau *= 1./curP;
Vector tauNext(tau);
Vector drVec(0,0,0);
dS = dir == alongMomentum ? fabs(dS) : -fabs(dS);
double radX0 = 1e24;
switch (fancy){
case HEL_AS_F:
case HEL_ALL_F:{
double p0 = curP;// see above = svCurrent.p3.mag();
double p0Inv = 1./p0;
double b0 = svCurrent.bf.mag();
//get to the mid-point first
double phi = (2.99792458e-3*svCurrent.q*b0*p0Inv)*dS/2.;
bool phiSmall = fabs(phi) < 1e-4;
double cosPhi = 0;
double sinPhi = 0;
double oneLessCosPhi=0;
double oneLessCosPhiOPhi=0;
double sinPhiOPhi=0;
double phiLessSinPhiOPhi=0;
if (phiSmall){
double phi2 = phi*phi;
double phi3 = phi2*phi;
double phi4 = phi3*phi;
sinPhi = phi - phi3/6. + phi4*phi/120.;
cosPhi = 1. -phi2/2. + phi4/24.;
oneLessCosPhi = phi2/2. - phi4/24. + phi2*phi4/720.; // 0.5*phi*phi;//*(1.- phi*phi/12.);
oneLessCosPhiOPhi = 0.5*phi - phi3/24. + phi2*phi3/720.;//*(1.- phi*phi/12.);
sinPhiOPhi = 1. - phi*phi/6. + phi4/120.;
phiLessSinPhiOPhi = phi*phi/6. - phi4/120. + phi4*phi2/5040.;//*(1. - phi*phi/20.);
} else {
cosPhi = cos(phi);
sinPhi = sin(phi);
oneLessCosPhi = 1.-cosPhi;
oneLessCosPhiOPhi = oneLessCosPhi/phi;
sinPhiOPhi = sinPhi/phi;
phiLessSinPhiOPhi = 1 - sinPhiOPhi;
}
Vector bHat = svCurrent.bf; bHat *= 1./b0; //bHat.mag();
// bool bAlongZ = fabs(bHat.z()) > 0.9999;
Vector btVec(bHat.cross(tau)); // for balong z btVec.z()==0
double tauB = tau.dot(bHat);
Vector bbtVec(bHat*tauB - tau); // (-tau.x(), -tau.y(), 0)
//don't need it here tauNext = tau + bbtVec*oneLessCosPhi - btVec*sinPhi;
drVec = bbtVec*phiLessSinPhiOPhi; drVec -= btVec*oneLessCosPhiOPhi; drVec += tau;
drVec *= dS/2.;
double dEdx = svCurrent.dEdx;
double dEdXPrime = svCurrent.dEdXPrime;
radX0 = svCurrent.radX0;
dP = dEdx*dS;
//improve with above values:
drVec += svCurrent.r3;
GlobalVector bfGV(0,0,0);
Vector bf(0,0,0); //(bfGV.x(), bfGV.y(), bfGV.z());
// = svCurrent.magVol->inTesla(GlobalPoint(drVec.x(), drVec.y(), -fabs(drVec.z())));
if (useMagVolumes_ && svCurrent.magVol != 0 && ! useInTeslaFromMagField_){
// this negative-z business will break at some point
if (vbField_->isZSymmetric()){
bfGV = svCurrent.magVol->inTesla(GlobalPoint(drVec.x(), drVec.y(), -fabs(drVec.z())));
} else {
bfGV = svCurrent.magVol->inTesla(GlobalPoint(drVec.x(), drVec.y(), drVec.z()));
}
if (drVec.z() > 0 && vbField_->isZSymmetric()){
bf.set(-bfGV.x(), -bfGV.y(), bfGV.z());
} else {
bf.set(bfGV.x(), bfGV.y(), bfGV.z());
}
} else {
bfGV = field_->inTesla(GlobalPoint(drVec.x(), drVec.y(), drVec.z()));
bf.set(bfGV.x(), bfGV.y(), bfGV.z());
}
double b0Init = b0;
b0 = bf.mag();
if (b0 < 1e-16) {
b0 = 1e-16;
bf.set(0., 0., 1e-16);
}
if (debug_){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Improved b "<<b0
<<" at r3 "<<drVec<<std::endl;
}
if (fabs((b0-b0Init)*dS) > 1){
//missed the mag volume boundary?
if (debug_){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Large bf*dS change "<<fabs((b0-svCurrent.bf.mag())*dS)
<<" --> recalc dedx"<<std::endl;
}
svNext.r3 = drVec;
svNext.bf = bf;
svNext.p3 = svCurrent.p3;
svNext.isYokeVol = svCurrent.isYokeVol;
MatBounds rzTmp;
dEdx = getDeDx(svNext, dEdXPrime, radX0, rzTmp);
dP = dEdx*dS;
if (debug_){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"New dEdX= "<<dEdx
<<" dP= "<<dP
<<" for p0 "<<p0<<std::endl;
}
}
//p0 is mid-way and b0 from mid-point
p0 += dP/2.; if (p0 < 1e-2) p0 = 1e-2;
p0Inv = 1./p0;
phi = (2.99792458e-3*svCurrent.q*b0*p0Inv)*dS;
phiSmall = fabs(phi) < 1e-4;
if (phiSmall){
double phi2 = phi*phi;
double phi3 = phi2*phi;
double phi4 = phi3*phi;
sinPhi = phi - phi3/6. + phi4*phi/120.;
cosPhi = 1. -phi2/2. + phi4/24.;
oneLessCosPhi = phi2/2. - phi4/24. + phi2*phi4/720.; // 0.5*phi*phi;//*(1.- phi*phi/12.);
oneLessCosPhiOPhi = 0.5*phi - phi3/24. + phi2*phi3/720.;//*(1.- phi*phi/12.);
sinPhiOPhi = 1. - phi*phi/6. + phi4/120.;
phiLessSinPhiOPhi = phi*phi/6. - phi4/120. + phi4*phi2/5040.;//*(1. - phi*phi/20.);
}else {
cosPhi = cos(phi);
sinPhi = sin(phi);
oneLessCosPhi = 1.-cosPhi;
oneLessCosPhiOPhi = oneLessCosPhi/phi;
sinPhiOPhi = sinPhi/phi;
phiLessSinPhiOPhi = 1. - sinPhiOPhi;
}
bHat = bf; bHat *= 1./b0;// as above =1./bHat.mag();
// bAlongZ = fabs(bHat.z()) > 0.9999;
btVec = bHat.cross(tau); // for b||z (-tau.y(), tau.x() ,0)
tauB = tau.dot(bHat);
bbtVec = bHat*tauB - tau; //bHat.cross(btVec); for b||z: (-tau.x(), -tau.y(), 0)
tauNext = bbtVec*oneLessCosPhi; tauNext -= btVec*sinPhi; tauNext += tau; //for b||z tauNext.z() == tau.z()
double tauNextPerpInv = 1./tauNext.perp();
drVec = bbtVec*phiLessSinPhiOPhi; drVec -= btVec*oneLessCosPhiOPhi; drVec += tau;
drVec *= dS;
if (svCurrent.hasErrorPropagated_){
double theta02 = 0;
double dX0 = fabs(dS)/radX0;
if (applyRadX0Correction_){
// this provides the integrand for theta^2
// if summed up along the path, should result in
// theta_total^2 = Int_0^x0{ f(x)dX} = (13.6/p0)^2*x0*(1+0.036*ln(x0+1))
// x0+1 above is to make the result infrared safe.
double x0 = fabs(svCurrent.radPath());
double alphaX0 = 13.6e-3*p0Inv; alphaX0 *= alphaX0;
double betaX0 = 0.038;
double logx0p1 = log(x0+1);
theta02 = dX0*alphaX0*(1+betaX0*logx0p1)*(1 + betaX0*logx0p1 + 2.*betaX0*x0/(x0+1) );
} else {
theta02 = 196e-6* p0Inv * p0Inv * dX0; //14.e-3/p0*sqrt(fabs(dS)/radX0); // .. drop log term (this is non-additive)
}
double epsilonP0 = 1.+ dP/(p0-0.5*dP);
// double omegaP0 = -dP/(p0-0.5*dP) + dS*dEdXPrime;
// double dsp = dS/(p0-0.5*dP); //use the initial p0 (not the mid-point) to keep the transport properly additive
Vector tbtVec(bHat - tauB*tau); // for b||z tau.z()*(-tau.x(), -tau.y(), 1.-tau.z())
dCCurvTransform_ = unit55_;
{
//Slightly modified copy of the curvilinear jacobian (don't use the original just because it's in float precision
// and seems to have some assumptions about the field values
// notation changes: p1--> tau, p2-->tauNext
// theta --> phi
// Vector p1 = tau;
// Vector p2 = tauNext;
Point xStart = svCurrent.r3;
Vector dx = drVec;
//GlobalVector h = MagneticField::inInverseGeV(xStart);
// Martijn: field is now given as parameter.. GlobalVector h = globalParameters.magneticFieldInInverseGeV(xStart);
//double qbp = fts.signedInverseMomentum();
double qbp = svCurrent.q*p0Inv;
// double absS = dS;
// calculate transport matrix
// Origin: TRPRFN
double t11 = tau.x(); double t12 = tau.y(); double t13 = tau.z();
double t21 = tauNext.x(); double t22 = tauNext.y(); double t23 = tauNext.z();
double cosl0 = tau.perp();
// double cosl1 = 1./tauNext.perp(); //not quite a cos .. it's a cosec--> change to cosecl1 below
double cosecl1 = tauNextPerpInv;
//AlgebraicMatrix a(5,5,1);
// define average magnetic field and gradient
// at initial point - inlike TRPRFN
Vector hn = bHat;
// double qp = -2.99792458e-3*b0;
// double q = -h.mag()*qbp;
double q = -phi/dS; //qp*qbp; // -phi/dS
// double theta = -phi;
double sint = -sinPhi; double cost = cosPhi;
double hn1 = hn.x(); double hn2 = hn.y(); double hn3 = hn.z();
double dx1 = dx.x(); double dx2 = dx.y(); double dx3 = dx.z();
// double hnDt1 = hn1*t11 + hn2*t12 + hn3*t13;
double gamma = hn1*t21 + hn2*t22 + hn3*t23;
double an1 = hn2*t23 - hn3*t22;
double an2 = hn3*t21 - hn1*t23;
double an3 = hn1*t22 - hn2*t21;
// double auInv = sqrt(1.- t13*t13); double au = auInv>0 ? 1./auInv : 1e24;
double auInv = cosl0; double au = auInv>0 ? 1./auInv : 1e24;
// double auInv = sqrt(t11*t11 + t12*t12); double au = auInv>0 ? 1./auInv : 1e24;
double u11 = -au*t12; double u12 = au*t11;
double v11 = -t13*u12; double v12 = t13*u11; double v13 = auInv;//t11*u12 - t12*u11;
auInv = sqrt(1. - t23*t23); au = auInv>0 ? 1./auInv : 1e24;
// auInv = sqrt(t21*t21 + t22*t22); au = auInv>0 ? 1./auInv : 1e24;
double u21 = -au*t22; double u22 = au*t21;
double v21 = -t23*u22; double v22 = t23*u21; double v23 = auInv;//t21*u22 - t22*u21;
// now prepare the transport matrix
// pp only needed in high-p case (WA)
// double pp = 1./qbp;
// moved up (where -h.mag() is needed()
// double qp = q*pp;
double anv = -(hn1*u21 + hn2*u22 );
double anu = (hn1*v21 + hn2*v22 + hn3*v23);
double omcost = oneLessCosPhi; double tmsint = -phi*phiLessSinPhiOPhi;
double hu1 = - hn3*u12;
double hu2 = hn3*u11;
double hu3 = hn1*u12 - hn2*u11;
double hv1 = hn2*v13 - hn3*v12;
double hv2 = hn3*v11 - hn1*v13;
double hv3 = hn1*v12 - hn2*v11;
// 1/p - doesn't change since |tau| = |tauNext| ... not. It changes now
dCCurvTransform_(0,0) = 1./(epsilonP0*epsilonP0)*(1. + dS*dEdXPrime);
// lambda
dCCurvTransform_(1,0) = phi*p0/svCurrent.q*cosecl1*
(sinPhi*bbtVec.z() - cosPhi*btVec.z());
//was dCCurvTransform_(1,0) = -qp*anv*(t21*dx1 + t22*dx2 + t23*dx3); //NOTE (SK) this was found to have an opposite sign
//from independent re-calculation ... in fact the tauNext.dot.dR piece isnt reproduced
dCCurvTransform_(1,1) = cost*(v11*v21 + v12*v22 + v13*v23) +
sint*(hv1*v21 + hv2*v22 + hv3*v23) +
omcost*(hn1*v11 + hn2*v12 + hn3*v13) * (hn1*v21 + hn2*v22 + hn3*v23) +
anv*(-sint*(v11*t21 + v12*t22 + v13*t23) +
omcost*(v11*an1 + v12*an2 + v13*an3) -
tmsint*gamma*(hn1*v11 + hn2*v12 + hn3*v13) );
dCCurvTransform_(1,2) = cost*(u11*v21 + u12*v22 ) +
sint*(hu1*v21 + hu2*v22 + hu3*v23) +
omcost*(hn1*u11 + hn2*u12 ) * (hn1*v21 + hn2*v22 + hn3*v23) +
anv*(-sint*(u11*t21 + u12*t22 ) +
omcost*(u11*an1 + u12*an2 ) -
tmsint*gamma*(hn1*u11 + hn2*u12 ) );
dCCurvTransform_(1,2) *= cosl0;
// Commented out in part for reproducibility purposes: these terms are zero in cart->curv
// dCCurvTransform_(1,3) = -q*anv*(u11*t21 + u12*t22 ); //don't show up in cartesian setup-->curv
//why would lambdaNext depend explicitely on initial position ? any arbitrary init point can be chosen not
// affecting the final state's momentum direction ... is this the field gradient in curvilinear coord?
// dCCurvTransform_(1,4) = -q*anv*(v11*t21 + v12*t22 + v13*t23); //don't show up in cartesian setup-->curv
// phi
dCCurvTransform_(2,0) = - phi*p0/svCurrent.q*cosecl1*cosecl1*
(oneLessCosPhi*bHat.z()*btVec.mag2() + sinPhi*btVec.z() + cosPhi*tbtVec.z()) ;
//was dCCurvTransform_(2,0) = -qp*anu*(t21*dx1 + t22*dx2 + t23*dx3)*cosecl1;
dCCurvTransform_(2,1) = cost*(v11*u21 + v12*u22 ) +
sint*(hv1*u21 + hv2*u22 ) +
omcost*(hn1*v11 + hn2*v12 + hn3*v13) *
(hn1*u21 + hn2*u22 ) +
anu*(-sint*(v11*t21 + v12*t22 + v13*t23) +
omcost*(v11*an1 + v12*an2 + v13*an3) -
tmsint*gamma*(hn1*v11 + hn2*v12 + hn3*v13) );
dCCurvTransform_(2,1) *= cosecl1;
dCCurvTransform_(2,2) = cost*(u11*u21 + u12*u22 ) +
sint*(hu1*u21 + hu2*u22 ) +
omcost*(hn1*u11 + hn2*u12 ) *
(hn1*u21 + hn2*u22 ) +
anu*(-sint*(u11*t21 + u12*t22 ) +
omcost*(u11*an1 + u12*an2 ) -
tmsint*gamma*(hn1*u11 + hn2*u12 ) );
dCCurvTransform_(2,2) *= cosecl1*cosl0;
// Commented out in part for reproducibility purposes: these terms are zero in cart->curv
// dCCurvTransform_(2,3) = -q*anu*(u11*t21 + u12*t22 )*cosecl1;
//why would lambdaNext depend explicitely on initial position ? any arbitrary init point can be chosen not
// affecting the final state's momentum direction ... is this the field gradient in curvilinear coord?
// dCCurvTransform_(2,4) = -q*anu*(v11*t21 + v12*t22 + v13*t23)*cosecl1;
// yt
double pp = 1./qbp;
// (SK) these terms seem to consistently have a sign opp from private derivation
dCCurvTransform_(3,0) = - pp*(u21*dx1 + u22*dx2 ); //NB: modified from the original: changed the sign
dCCurvTransform_(4,0) = - pp*(v21*dx1 + v22*dx2 + v23*dx3);
dCCurvTransform_(3,1) = (sint*(v11*u21 + v12*u22 ) +
omcost*(hv1*u21 + hv2*u22 ) +
tmsint*(hn1*u21 + hn2*u22 ) *
(hn1*v11 + hn2*v12 + hn3*v13))/q;
dCCurvTransform_(3,2) = (sint*(u11*u21 + u12*u22 ) +
omcost*(hu1*u21 + hu2*u22 ) +
tmsint*(hn1*u21 + hn2*u22 ) *
(hn1*u11 + hn2*u12 ))*cosl0/q;
dCCurvTransform_(3,3) = (u11*u21 + u12*u22 );
dCCurvTransform_(3,4) = (v11*u21 + v12*u22 );
// zt
dCCurvTransform_(4,1) = (sint*(v11*v21 + v12*v22 + v13*v23) +
omcost*(hv1*v21 + hv2*v22 + hv3*v23) +
tmsint*(hn1*v21 + hn2*v22 + hn3*v23) *
(hn1*v11 + hn2*v12 + hn3*v13))/q;
dCCurvTransform_(4,2) = (sint*(u11*v21 + u12*v22 ) +
omcost*(hu1*v21 + hu2*v22 + hu3*v23) +
tmsint*(hn1*v21 + hn2*v22 + hn3*v23) *
(hn1*u11 + hn2*u12 ))*cosl0/q;
dCCurvTransform_(4,3) = (u11*v21 + u12*v22 );
dCCurvTransform_(4,4) = (v11*v21 + v12*v22 + v13*v23);
// end of TRPRFN
}
if (debug_){
Basis rep; setRep(rep, tauNext);
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"rep X: "<<rep.lX<<" "<<rep.lX.mag()
<<"\t Y: "<<rep.lY<<" "<<rep.lY.mag()
<<"\t Z: "<<rep.lZ<<" "<<rep.lZ.mag();
}
//mind the sign of dS and dP (dS*dP < 0 allways)
//covariance should grow no matter which direction you propagate
//==> take abs values.
//reset not needed: fill all below svCurrent.matDCov *= 0.;
double mulRR = theta02*dS*dS/3.;
double mulRP = theta02*fabs(dS)*p0/2.;
double mulPP = theta02*p0*p0;
double losPP = dP*dP*1.6/fabs(dS)*(1.0 + p0*1e-3);
//another guess .. makes sense for 1 cm steps 2./dS == 2 [cm] / dS [cm] at low pt
//double it by 1TeV
//not gaussian anyways
// derived from the fact that sigma_p/eLoss ~ 0.08 after ~ 200 steps
//curvilinear
double sinThetaInv = tauNextPerpInv;
double p0Mat = p0+ 0.5*dP; // FIXME change this to p0 after it's clear that there's agreement in everything else
double p0Mat2 = p0Mat*p0Mat;
// with 6x6 formulation
svCurrent.matDCovCurv*=0;
svCurrent.matDCovCurv(0,0) = losPP/(p0Mat2*p0Mat2);
// svCurrent.matDCovCurv(0,1) = 0;
// svCurrent.matDCovCurv(0,2) = 0;
// svCurrent.matDCovCurv(0,3) = 0;
// svCurrent.matDCovCurv(0,4) = 0;
// svCurrent.matDCovCurv(1,0) = 0;
svCurrent.matDCovCurv(1,1) = mulPP/p0Mat2;
// svCurrent.matDCovCurv(1,2) = 0;
// svCurrent.matDCovCurv(1,3) = 0;
svCurrent.matDCovCurv(1,4) = mulRP/p0Mat;
// svCurrent.matDCovCurv(2,0) = 0;
// svCurrent.matDCovCurv(2,1) = 0;
svCurrent.matDCovCurv(2,2) = mulPP/p0Mat2*(sinThetaInv*sinThetaInv);
svCurrent.matDCovCurv(2,3) = mulRP/p0Mat*sinThetaInv;
// svCurrent.matDCovCurv(2,4) = 0;
// svCurrent.matDCovCurv(3,0) = 0;
// svCurrent.matDCovCurv(3,1) = 0;
svCurrent.matDCovCurv(3,2) = mulRP/p0Mat*sinThetaInv;
// svCurrent.matDCovCurv(3,0) = 0;
svCurrent.matDCovCurv(3,3) = mulRR;
// svCurrent.matDCovCurv(3,4) = 0;
// svCurrent.matDCovCurv(4,0) = 0;
svCurrent.matDCovCurv(4,1) = mulRP/p0Mat;
// svCurrent.matDCovCurv(4,2) = 0;
// svCurrent.matDCovCurv(4,3) = 0;
svCurrent.matDCovCurv(4,4) = mulRR;
}
break;
}
default:
break;
}
double pMag = curP;//svCurrent.p3.mag();
if (dir == oppositeToMomentum) dP = fabs(dP);
else if( dP < 0) { //the case of negative dP
dP = -dP > pMag ? -pMag+1e-5 : dP;
}
getNextState(svCurrent, svNext, dP, tauNext, drVec, dS, dS/radX0,
dCCurvTransform_);
return true;
}
| const SteppingHelixStateInfo & SteppingHelixPropagator::propagate | ( | const SteppingHelixStateInfo & | ftsStart, |
| const GlobalPoint & | pDest | ||
| ) | const |
Propagate to PCA to point given a starting point.
Definition at line 267 of file SteppingHelixPropagator.cc.
References cIndex_(), invalidState_, SteppingHelixStateInfo::isValid(), nPoints_, POINT_PCA_DT, propagate(), sendLogWarning_, setIState(), svBuf_, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().
{
if (! sStart.isValid()){
if (sendLogWarning_){
edm::LogWarning("SteppingHelixPropagator")<<"Can't propagate: invalid input state"
<<std::endl;
}
return invalidState_;
}
setIState(sStart);
double pars[6] = {pDest.x(), pDest.y(), pDest.z(), 0, 0, 0};
propagate(POINT_PCA_DT, pars);
return svBuf_[cIndex_(nPoints_-1)];
}
| SteppingHelixPropagator::Result SteppingHelixPropagator::propagate | ( | 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 327 of file SteppingHelixPropagator.cc.
References alongMomentum, anyDirection, cIndex_(), debug_, SteppingHelixStateInfo::dEdx, defaultStep_, dir, SteppingHelixStateInfo::FAULT, SteppingHelixStateInfo::field, field_, HEL_AS_F, SteppingHelixStateInfo::INACC, SteppingHelixStateInfo::isValid_, LINE_PCA_DT, LogTrace, makeAtomStep(), MAX_STEPS, metname, nPoints_, 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(), query::result, SteppingHelixStateInfo::ResultName, sendLogWarning_, SteppingHelixStateInfo::status_, svBuf_, pftools::UNDEFINED, useIsYokeFlag_, useMagVolumes_, useMatVolumes_, useTuningForL2Speed_, SteppingHelixStateInfo::WRONG_VOLUME, Z_DT, and Z_P.
{
static const std::string metname = "SteppingHelixPropagator";
StateInfo* svCurrent = &svBuf_[cIndex_(nPoints_-1)];
//check if it's going to work at all
double tanDist = 0;
double dist = 0;
PropagationDirection refDirection = anyDirection;
Result result = refToDest(type, (*svCurrent), pars, dist, tanDist, refDirection);
if (result != SteppingHelixStateInfo::OK || fabs(dist) > 1e6){
svCurrent->status_ = result;
if (fabs(dist) > 1e6) svCurrent->status_ = SteppingHelixStateInfo::INACC;
svCurrent->isValid_ = false;
svCurrent->field = field_;
if (sendLogWarning_){
edm::LogWarning(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<" Failed after first refToDest check with status "
<<SteppingHelixStateInfo::ResultName[result]
<<std::endl;
}
return result;
}
result = SteppingHelixStateInfo::UNDEFINED;
bool makeNextStep = true;
double dStep = defaultStep_;
PropagationDirection dir,oldDir;
dir = propagationDirection();
oldDir = dir;
int nOsc = 0;
double distMag = 1e12;
double tanDistMag = 1e12;
double distMat = 1e12;
double tanDistMat = 1e12;
double tanDistNextCheck = -0.1;//just need a negative start val
double tanDistMagNextCheck = -0.1;
double tanDistMatNextCheck = -0.1;
double oldDStep = 0;
PropagationDirection oldRefDirection = propagationDirection();
Result resultToMat = SteppingHelixStateInfo::UNDEFINED;
Result resultToMag = SteppingHelixStateInfo::UNDEFINED;
bool isFirstStep = true;
bool expectNewMagVolume = false;
int loopCount = 0;
while (makeNextStep){
dStep = defaultStep_;
svCurrent = &svBuf_[cIndex_(nPoints_-1)];
double curZ = svCurrent->r3.z();
double curR = svCurrent->r3.perp();
if ( fabs(curZ) < 440 && curR < 260) dStep = defaultStep_*2;
//more such ifs might be scattered around
//even though dStep is large, it will still make stops at each volume boundary
if (useTuningForL2Speed_){
dStep = 100.;
if (! useIsYokeFlag_ && fabs(curZ) < 667 && curR > 380 && curR < 850){
dStep = 5*(1+0.2*svCurrent->p3.mag());
}
}
// refDirection = propagationDirection();
tanDistNextCheck -= oldDStep;
tanDistMagNextCheck -= oldDStep;
tanDistMatNextCheck -= oldDStep;
if (tanDistNextCheck < 0){
//use pre-computed values if it's the first step
if (! isFirstStep) refToDest(type, (*svCurrent), pars, dist, tanDist, refDirection);
// constrain allowed path for a tangential approach
if (fabs(tanDist) > 4.*(fabs(dist)) ) tanDist *= tanDist == 0 ? 0 :fabs(dist/tanDist*4.);
tanDistNextCheck = fabs(tanDist)*0.5 - 0.5; //need a better guess (to-do)
//reasonable limit
if (tanDistNextCheck > defaultStep_*20. ) tanDistNextCheck = defaultStep_*20.;
oldRefDirection = refDirection;
} else {
tanDist = tanDist > 0. ? tanDist - oldDStep : tanDist + oldDStep;
refDirection = oldRefDirection;
if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Skipped refToDest: guess tanDist = "<<tanDist
<<" next check at "<<tanDistNextCheck<<std::endl;
}
double fastSkipDist = fabs(dist) > fabs(tanDist) ? tanDist : dist;
if (propagationDirection() == anyDirection){
dir = refDirection;
} else {
dir = propagationDirection();
if (fabs(tanDist)<0.1 && refDirection != dir ){
//how did it get here? nOsc++;
dir = refDirection;
if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"NOTE: overstepped last time: switch direction (can do it if within 1 mm)"<<std::endl;
}
}
if (useMagVolumes_ && ! (fabs(dist) < fabs(epsilon))){//need to know the general direction
if (tanDistMagNextCheck < 0){
resultToMag = refToMagVolume((*svCurrent), dir, distMag, tanDistMag, fabs(fastSkipDist), expectNewMagVolume, fabs(tanDist));
// constrain allowed path for a tangential approach
if (fabs(tanDistMag) > 4.*(fabs(distMag)) ) tanDistMag *= tanDistMag == 0 ? 0 : fabs(distMag/tanDistMag*4.);
tanDistMagNextCheck = fabs(tanDistMag)*0.5-0.5; //need a better guess (to-do)
//reasonable limit; "turn off" checking if bounds are further than the destination
if (tanDistMagNextCheck > defaultStep_*20.
|| fabs(dist) < fabs(distMag)
|| resultToMag ==SteppingHelixStateInfo::INACC)
tanDistMagNextCheck = defaultStep_*20 > fabs(fastSkipDist) ? fabs(fastSkipDist) : defaultStep_*20;
if (resultToMag != SteppingHelixStateInfo::INACC
&& resultToMag != SteppingHelixStateInfo::OK) tanDistMagNextCheck = -1;
} else {
// resultToMag = SteppingHelixStateInfo::OK;
tanDistMag = tanDistMag > 0. ? tanDistMag - oldDStep : tanDistMag + oldDStep;
if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Skipped refToMag: guess tanDistMag = "<<tanDistMag
<<" next check at "<<tanDistMagNextCheck;
}
}
if (useMatVolumes_ && ! (fabs(dist) < fabs(epsilon))){//need to know the general direction
if (tanDistMatNextCheck < 0){
resultToMat = refToMatVolume((*svCurrent), dir, distMat, tanDistMat, fabs(fastSkipDist));
// constrain allowed path for a tangential approach
if (fabs(tanDistMat) > 4.*(fabs(distMat)) ) tanDistMat *= tanDistMat == 0 ? 0 : fabs(distMat/tanDistMat*4.);
tanDistMatNextCheck = fabs(tanDistMat)*0.5-0.5; //need a better guess (to-do)
//reasonable limit; "turn off" checking if bounds are further than the destination
if (tanDistMatNextCheck > defaultStep_*20.
|| fabs(dist) < fabs(distMat)
|| resultToMat ==SteppingHelixStateInfo::INACC )
tanDistMatNextCheck = defaultStep_*20 > fabs(fastSkipDist) ? fabs(fastSkipDist) : defaultStep_*20;
if (resultToMat != SteppingHelixStateInfo::INACC
&& resultToMat != SteppingHelixStateInfo::OK) tanDistMatNextCheck = -1;
} else {
// resultToMat = SteppingHelixStateInfo::OK;
tanDistMat = tanDistMat > 0. ? tanDistMat - oldDStep : tanDistMat + oldDStep;
if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Skipped refToMat: guess tanDistMat = "<<tanDistMat
<<" next check at "<<tanDistMatNextCheck;
}
}
double rDotP = svCurrent->r3.dot(svCurrent->p3);
if ((fabs(curZ) > 1.5e3 || curR >800.)
&& ((dir == alongMomentum && rDotP > 0)
|| (dir == oppositeToMomentum && rDotP < 0) )
){
dStep = fabs(tanDist) -1e-12;
}
double tanDistMin = fabs(tanDist);
if (tanDistMin > fabs(tanDistMag)+0.05 &&
(resultToMag == SteppingHelixStateInfo::OK || resultToMag == SteppingHelixStateInfo::WRONG_VOLUME)){
tanDistMin = fabs(tanDistMag)+0.05; //try to step into the next volume
expectNewMagVolume = true;
} else expectNewMagVolume = false;
if (tanDistMin > fabs(tanDistMat)+0.05 && resultToMat == SteppingHelixStateInfo::OK){
tanDistMin = fabs(tanDistMat)+0.05; //try to step into the next volume
if (expectNewMagVolume) expectNewMagVolume = false;
}
if (tanDistMin*fabs(svCurrent->dEdx) > 0.5*svCurrent->p3.mag()){
tanDistMin = 0.5*svCurrent->p3.mag()/fabs(svCurrent->dEdx);
if (expectNewMagVolume) expectNewMagVolume = false;
}
double tanDistMinLazy = fabs(tanDistMin);
if ((type == POINT_PCA_DT || type == LINE_PCA_DT)
&& fabs(tanDist) < 2.*fabs(tanDistMin) ){
//being lazy here; the best is to take into account the curvature
tanDistMinLazy = fabs(tanDistMin)*0.5;
}
if (fabs(tanDistMinLazy) < dStep){
dStep = fabs(tanDistMinLazy);
}
//keep this path length for the next step
oldDStep = dStep;
if (dStep > 1e-10 && ! (fabs(dist) < fabs(epsilon))){
StateInfo* svNext = &svBuf_[cIndex_(nPoints_)];
makeAtomStep((*svCurrent), (*svNext), dStep, dir, HEL_AS_F);
// if (useMatVolumes_ && expectNewMagVolume
// && svCurrent->magVol == svNext->magVol){
// double tmpDist=0;
// double tmpDistMag = 0;
// if (refToMagVolume((*svNext), dir, tmpDist, tmpDistMag, fabs(dist)) != SteppingHelixStateInfo::OK){
// //the point appears to be outside, but findVolume claims the opposite
// dStep += 0.05*fabs(tanDistMag/distMag); oldDStep = dStep; //do it again with a bigger step
// if (debug_) LogTrace(metname)
// <<"Failed to get into new mag volume: will try with new bigger step "<<dStep<<std::endl;
// makeAtomStep((*svCurrent), (*svNext), dStep, dir, HEL_AS_F);
// }
// }
nPoints_++; svCurrent = &svBuf_[cIndex_(nPoints_-1)];
if (oldDir != dir){
nOsc++;
tanDistNextCheck = -1;//check dist after osc
tanDistMagNextCheck = -1;
tanDistMatNextCheck = -1;
}
oldDir = dir;
}
if (nOsc>1 && fabs(dStep)>epsilon){
if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Ooops"<<std::endl;
}
if (fabs(dist) < fabs(epsilon) ) result = SteppingHelixStateInfo::OK;
if ((type == POINT_PCA_DT || type == LINE_PCA_DT )
&& fabs(dStep) < fabs(epsilon) ){
//now check if it's not a branch point (peek ahead at 1 cm)
double nextDist = 0;
double nextTanDist = 0;
PropagationDirection nextRefDirection = anyDirection;
StateInfo* svNext = &svBuf_[cIndex_(nPoints_)];
makeAtomStep((*svCurrent), (*svNext), 1., dir, HEL_AS_F);
nPoints_++; svCurrent = &svBuf_[cIndex_(nPoints_-1)];
refToDest(type, (*svCurrent), pars, nextDist, nextTanDist, nextRefDirection);
if ( fabs(nextDist) > fabs(dist)){
nPoints_--; svCurrent = &svBuf_[cIndex_(nPoints_-1)];
result = SteppingHelixStateInfo::OK;
if (debug_){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Found real local minimum in PCA"<<std::endl;
}
} else {
//keep this trial point and continue
dStep = defaultStep_;
if (debug_){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Found branch point in PCA"<<std::endl;
}
}
}
if (nPoints_ > MAX_STEPS*1./defaultStep_ || loopCount > MAX_STEPS*100
|| nOsc > 6 ) result = SteppingHelixStateInfo::FAULT;
if (svCurrent->p3.mag() < 0.1 ) result = SteppingHelixStateInfo::RANGEOUT;
curZ = svCurrent->r3.z();
curR = svCurrent->r3.perp();
if ( curR > 20000 || fabs(curZ) > 20000 ) result = SteppingHelixStateInfo::INACC;
makeNextStep = result == SteppingHelixStateInfo::UNDEFINED;
svCurrent->status_ = result;
svCurrent->isValid_ = (result == SteppingHelixStateInfo::OK || makeNextStep );
svCurrent->field = field_;
isFirstStep = false;
loopCount++;
}
if (sendLogWarning_ && result != SteppingHelixStateInfo::OK){
edm::LogWarning(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<" Propagation failed with status "
<<SteppingHelixStateInfo::ResultName[result]
<<std::endl;
if (result == SteppingHelixStateInfo::RANGEOUT)
edm::LogWarning(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<" Momentum at last point is too low (<0.1) p_last = "
<<svCurrent->p3.mag()
<<std::endl;
if (result == SteppingHelixStateInfo::INACC)
edm::LogWarning(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<" Went too far: the last point is at "<<svCurrent->r3
<<std::endl;
if (result == SteppingHelixStateInfo::FAULT && nOsc > 6)
edm::LogWarning(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<" Infinite loop condidtion detected: going in cycles. Break after 6 cycles"
<<std::endl;
if (result == SteppingHelixStateInfo::FAULT && nPoints_ > MAX_STEPS*1./defaultStep_)
edm::LogWarning(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<" Tired to go farther. Made too many steps: more than "
<<MAX_STEPS*1./defaultStep_
<<std::endl;
}
if (debug_){
switch (type) {
case RADIUS_DT:
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"going to radius "<<pars[RADIUS_P]<<std::endl;
break;
case Z_DT:
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"going to z "<<pars[Z_P]<<std::endl;
break;
case PATHL_DT:
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"going to pathL "<<pars[PATHL_P]<<std::endl;
break;
case PLANE_DT:
{
Point rPlane(pars[0], pars[1], pars[2]);
Vector nPlane(pars[3], pars[4], pars[5]);
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"going to plane r0:"<<rPlane<<" n:"<<nPlane<<std::endl;
}
break;
case POINT_PCA_DT:
{
Point rDest(pars[0], pars[1], pars[2]);
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"going to PCA to point "<<rDest<<std::endl;
}
break;
case LINE_PCA_DT:
{
Point rDest1(pars[0], pars[1], pars[2]);
Point rDest2(pars[3], pars[4], pars[5]);
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"going to PCA to line "<<rDest1<<" - "<<rDest2<<std::endl;
}
break;
default:
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"going to NOT IMPLEMENTED"<<std::endl;
break;
}
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Made "<<nPoints_-1<<" steps and stopped at(cur step) "<<svCurrent->r3<<" nOsc "<<nOsc<<std::endl;
}
return result;
}
| TrajectoryStateOnSurface SteppingHelixPropagator::propagate | ( | const FreeTrajectoryState & | ftsStart, |
| const Cylinder & | cDest | ||
| ) | const [virtual] |
Propagate to Cylinder given a starting point (a Cylinder is assumed to be positioned at 0,0,0)
Implements Propagator.
Definition at line 92 of file SteppingHelixPropagator.cc.
References propagateWithPath().
{
return propagateWithPath(ftsStart, cDest).first;
}
| const SteppingHelixStateInfo & SteppingHelixPropagator::propagate | ( | const SteppingHelixStateInfo & | ftsStart, |
| const GlobalPoint & | pDest1, | ||
| const GlobalPoint & | pDest2 | ||
| ) | const |
Propagate to PCA to a line (given by 2 points) given a starting point.
Definition at line 288 of file SteppingHelixPropagator.cc.
References cIndex_(), invalidState_, SteppingHelixStateInfo::isValid(), LINE_PCA_DT, mag(), nPoints_, propagate(), sendLogWarning_, setIState(), svBuf_, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().
{
if ((pDest1-pDest2).mag() < 1e-10 || !sStart.isValid()){
if (sendLogWarning_){
if ((pDest1-pDest2).mag() < 1e-10)
edm::LogWarning("SteppingHelixPropagator")<<"Can't propagate: points are too close"
<<std::endl;
if (!sStart.isValid())
edm::LogWarning("SteppingHelixPropagator")<<"Can't propagate: invalid input state"
<<std::endl;
}
return invalidState_;
}
setIState(sStart);
double pars[6] = {pDest1.x(), pDest1.y(), pDest1.z(),
pDest2.x(), pDest2.y(), pDest2.z()};
propagate(LINE_PCA_DT, pars);
return svBuf_[cIndex_(nPoints_-1)];
}
| FreeTrajectoryState SteppingHelixPropagator::propagate | ( | const FreeTrajectoryState & | ftsStart, |
| const GlobalPoint & | pDest | ||
| ) | const [virtual] |
Propagate to PCA to point given a starting point.
Definition at line 98 of file SteppingHelixPropagator.cc.
References propagateWithPath().
{
return propagateWithPath(ftsStart, pDest).first;
}
| TrajectoryStateOnSurface SteppingHelixPropagator::propagate | ( | const FreeTrajectoryState & | ftsStart, |
| const Plane & | pDest | ||
| ) | const [virtual] |
Propagate to Plane given a starting point.
Implements Propagator.
Definition at line 87 of file SteppingHelixPropagator.cc.
References propagateWithPath().
Referenced by AlCaHOCalibProducer::produce(), propagate(), and propagateWithPath().
{
return propagateWithPath(ftsStart, pDest).first;
}
| FreeTrajectoryState SteppingHelixPropagator::propagate | ( | const FreeTrajectoryState & | ftsStart, |
| const GlobalPoint & | pDest1, | ||
| const GlobalPoint & | pDest2 | ||
| ) | const [virtual] |
Propagate to PCA to a line (given by 2 points) given a starting point.
Definition at line 104 of file SteppingHelixPropagator.cc.
References propagateWithPath().
{
return propagateWithPath(ftsStart, pDest1, pDest2).first;
}
| const SteppingHelixStateInfo & SteppingHelixPropagator::propagate | ( | const SteppingHelixStateInfo & | ftsStart, |
| const Surface & | sDest | ||
| ) | const |
Propagate to Plane given a starting point.
Definition at line 187 of file SteppingHelixPropagator.cc.
References invalidState_, SteppingHelixStateInfo::isValid(), propagate(), and sendLogWarning_.
{
if (! sStart.isValid()){
if (sendLogWarning_){
edm::LogWarning("SteppingHelixPropagator")<<"Can't propagate: invalid input state"
<<std::endl;
}
return invalidState_;
}
const Plane* pDest = dynamic_cast<const Plane*>(&sDest);
if (pDest != 0) return propagate(sStart, *pDest);
const Cylinder* cDest = dynamic_cast<const Cylinder*>(&sDest);
if (cDest != 0) return propagate(sStart, *cDest);
throw PropagationException("The surface is neither Cylinder nor Plane");
}
| FreeTrajectoryState SteppingHelixPropagator::propagate | ( | const FreeTrajectoryState & | ftsStart, |
| const reco::BeamSpot & | beamSpot | ||
| ) | const [virtual] |
Propagate to PCA to a line determined by BeamSpot position and slope given a starting point.
Reimplemented from Propagator.
Definition at line 111 of file SteppingHelixPropagator.cc.
References propagateWithPath().
{
return propagateWithPath(ftsStart, beamSpot).first;
}
| const SteppingHelixStateInfo & SteppingHelixPropagator::propagate | ( | const SteppingHelixStateInfo & | ftsStart, |
| const Plane & | pDest | ||
| ) | const |
Definition at line 209 of file SteppingHelixPropagator.cc.
References cIndex_(), SteppingHelixStateInfo::dir, invalidState_, SteppingHelixStateInfo::isValid(), nPoints_, SteppingHelixStateInfo::path(), PLANE_DT, GloballyPositioned< T >::position(), propagate(), GloballyPositioned< T >::rotation(), sendLogWarning_, setIState(), and svBuf_.
{
if (! sStart.isValid()){
if (sendLogWarning_){
edm::LogWarning("SteppingHelixPropagator")<<"Can't propagate: invalid input state"
<<std::endl;
}
return invalidState_;
}
setIState(sStart);
Point rPlane(pDest.position().x(), pDest.position().y(), pDest.position().z());
Vector nPlane(pDest.rotation().zx(), pDest.rotation().zy(), pDest.rotation().zz()); nPlane /= nPlane.mag();
double pars[6] = { rPlane.x(), rPlane.y(), rPlane.z(),
nPlane.x(), nPlane.y(), nPlane.z() };
propagate(PLANE_DT, pars);
//(re)set it before leaving: dir =1 (-1) if path increased (decreased) and 0 if it didn't change
//need to implement this somewhere else as a separate function
double lDir = 0;
if (sStart.path() < svBuf_[cIndex_(nPoints_-1)].path()) lDir = 1.;
if (sStart.path() > svBuf_[cIndex_(nPoints_-1)].path()) lDir = -1.;
svBuf_[cIndex_(nPoints_-1)].dir = lDir;
return svBuf_[cIndex_(nPoints_-1)];
}
| const SteppingHelixStateInfo & SteppingHelixPropagator::propagate | ( | const SteppingHelixStateInfo & | ftsStart, |
| const Cylinder & | cDest | ||
| ) | const |
Propagate to Cylinder given a starting point (a Cylinder is assumed to be positioned at 0,0,0)
Definition at line 239 of file SteppingHelixPropagator.cc.
References cIndex_(), SteppingHelixStateInfo::dir, invalidState_, SteppingHelixStateInfo::isValid(), nPoints_, SteppingHelixStateInfo::path(), propagate(), Cylinder::radius(), RADIUS_DT, RADIUS_P, sendLogWarning_, setIState(), and svBuf_.
{
if (! sStart.isValid()){
if (sendLogWarning_){
edm::LogWarning("SteppingHelixPropagator")<<"Can't propagate: invalid input state"
<<std::endl;
}
return invalidState_;
}
setIState(sStart);
double pars[6] = {0,0,0,0,0,0};
pars[RADIUS_P] = cDest.radius();
propagate(RADIUS_DT, pars);
//(re)set it before leaving: dir =1 (-1) if path increased (decreased) and 0 if it didn't change
//need to implement this somewhere else as a separate function
double lDir = 0;
if (sStart.path() < svBuf_[cIndex_(nPoints_-1)].path()) lDir = 1.;
if (sStart.path() > svBuf_[cIndex_(nPoints_-1)].path()) lDir = -1.;
svBuf_[cIndex_(nPoints_-1)].dir = lDir;
return svBuf_[cIndex_(nPoints_-1)];
}
| std::pair< TrajectoryStateOnSurface, double > SteppingHelixPropagator::propagateWithPath | ( | const FreeTrajectoryState & | ftsStart, |
| const Plane & | pDest | ||
| ) | const [virtual] |
Propagate to Plane given a starting point: return final TrajectoryState and path length from start to this point
Implements Propagator.
Definition at line 119 of file SteppingHelixPropagator.cc.
References SteppingHelixStateInfo::getStateOnSurface(), SteppingHelixStateInfo::path(), propagate(), setIState(), and svBuf_.
Referenced by PropagateToCal::propagate(), propagate(), and propagateWithPath().
| std::pair< TrajectoryStateOnSurface, double > SteppingHelixPropagator::propagateWithPath | ( | const FreeTrajectoryState & | ftsStart, |
| const Cylinder & | cDest | ||
| ) | const [virtual] |
Propagate to Cylinder given a starting point: return final TrajectoryState and path length from start to this point
Implements Propagator.
Definition at line 130 of file SteppingHelixPropagator.cc.
References SteppingHelixStateInfo::getStateOnSurface(), SteppingHelixStateInfo::path(), propagate(), returnTangentPlane_, setIState(), and svBuf_.
{
setIState(SteppingHelixStateInfo(ftsStart));
const StateInfo& svCurrent = propagate(svBuf_[0], cDest);
return TsosPP(svCurrent.getStateOnSurface(cDest, returnTangentPlane_), svCurrent.path());
}
| std::pair< FreeTrajectoryState, double > SteppingHelixPropagator::propagateWithPath | ( | const FreeTrajectoryState & | ftsStart, |
| const GlobalPoint & | pDest | ||
| ) | const [virtual] |
Propagate to PCA to point given a starting point.
Definition at line 142 of file SteppingHelixPropagator.cc.
References SteppingHelixStateInfo::getFreeState(), SteppingHelixStateInfo::path(), propagate(), setIState(), and svBuf_.
{
setIState(SteppingHelixStateInfo(ftsStart));
const StateInfo& svCurrent = propagate(svBuf_[0], pDest);
FreeTrajectoryState ftsDest;
svCurrent.getFreeState(ftsDest);
return FtsPP(ftsDest, svCurrent.path());
}
| std::pair< FreeTrajectoryState, double > SteppingHelixPropagator::propagateWithPath | ( | const FreeTrajectoryState & | ftsStart, |
| const GlobalPoint & | pDest1, | ||
| const GlobalPoint & | pDest2 | ||
| ) | const [virtual] |
Propagate to PCA to a line (given by 2 points) given a starting point.
Reimplemented from Propagator.
Definition at line 155 of file SteppingHelixPropagator.cc.
References SteppingHelixStateInfo::getFreeState(), mag(), SteppingHelixStateInfo::path(), propagate(), sendLogWarning_, setIState(), and svBuf_.
{
if ((pDest1-pDest2).mag() < 1e-10){
if (sendLogWarning_){
edm::LogWarning("SteppingHelixPropagator")<<"Can't propagate: the points should be at a bigger distance"
<<std::endl;
}
return FtsPP();
}
setIState(SteppingHelixStateInfo(ftsStart));
const StateInfo& svCurrent = propagate(svBuf_[0], pDest1, pDest2);
FreeTrajectoryState ftsDest;
svCurrent.getFreeState(ftsDest);
return FtsPP(ftsDest, svCurrent.path());
}
| std::pair< FreeTrajectoryState, double > SteppingHelixPropagator::propagateWithPath | ( | const FreeTrajectoryState & | ftsStart, |
| const reco::BeamSpot & | beamSpot | ||
| ) | const [virtual] |
Propagate to PCA to a line (given by beamSpot position and slope) given a starting point.
Definition at line 177 of file SteppingHelixPropagator.cc.
References reco::BeamSpot::dxdz(), reco::BeamSpot::dydz(), propagateWithPath(), reco::BeamSpot::x0(), reco::BeamSpot::y0(), and reco::BeamSpot::z0().
{
GlobalPoint pDest1(beamSpot.x0(), beamSpot.y0(), beamSpot.z0());
GlobalPoint pDest2(pDest1.x() + beamSpot.dxdz()*1e3,
pDest1.y() + beamSpot.dydz()*1e3,
pDest1.z() + 1e3);
return propagateWithPath(ftsStart, pDest1, pDest2);
}
| 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 1691 of file SteppingHelixPropagator.cc.
References alongMomentum, anyDirection, SteppingHelixStateInfo::bf, CONE_DT, funct::cos(), debug_, i, SteppingHelixStateInfo::INACC, LINE_PCA_DT, LogTrace, mag(), metname, lumiNorm::norm, SteppingHelixStateInfo::NOT_IMPLEMENTED, SteppingHelixStateInfo::OK, oppositeToMomentum, SteppingHelixStateInfo::p3, SteppingHelixStateInfo::path(), PATHL_DT, PATHL_P, Geom::pi(), pi, PLANE_DT, POINT_PCA_DT, SteppingHelixStateInfo::q, SteppingHelixStateInfo::r3, RADIUS_DT, RADIUS_P, query::result, funct::sin(), mathSSE::sqrt(), metsig::tau, theta(), Z_DT, and Z_P.
Referenced by propagate(), refToMagVolume(), and refToMatVolume().
{
static const std::string metname = "SteppingHelixPropagator";
Result result = SteppingHelixStateInfo::NOT_IMPLEMENTED;
switch (dest){
case RADIUS_DT:
{
double curR = sv.r3.perp();
dist = pars[RADIUS_P] - curR;
if (fabs(dist) > fastSkipDist){
result = SteppingHelixStateInfo::INACC;
break;
}
double curP2 = sv.p3.mag2();
double curPtPos2 = sv.p3.perp2(); if(curPtPos2< 1e-16) curPtPos2=1e-16;
double cosDPhiPR = (sv.r3.x()*sv.p3.x()+sv.r3.y()*sv.p3.y());//only the sign is needed cos((sv.r3.deltaPhi(sv.p3)));
refDirection = dist*cosDPhiPR > 0 ?
alongMomentum : oppositeToMomentum;
tanDist = dist*sqrt(curP2/curPtPos2);
result = SteppingHelixStateInfo::OK;
}
break;
case Z_DT:
{
double curZ = sv.r3.z();
dist = pars[Z_P] - curZ;
if (fabs(dist) > fastSkipDist){
result = SteppingHelixStateInfo::INACC;
break;
}
double curP = sv.p3.mag();
refDirection = sv.p3.z()*dist > 0. ?
alongMomentum : oppositeToMomentum;
tanDist = dist/sv.p3.z()*curP;
result = SteppingHelixStateInfo::OK;
}
break;
case PLANE_DT:
{
Point rPlane(pars[0], pars[1], pars[2]);
Vector nPlane(pars[3], pars[4], pars[5]);
// unfortunately this doesn't work: the numbers are too large
// bool pVertical = fabs(pars[5])>0.9999;
// double dRDotN = pVertical? (sv.r3.z() - rPlane.z())*nPlane.z() :(sv.r3 - rPlane).dot(nPlane);
double dRDotN = (sv.r3.x()-rPlane.x())*nPlane.x() + (sv.r3.y()-rPlane.y())*nPlane.y() + (sv.r3.z()-rPlane.z())*nPlane.z();//(sv.r3 - rPlane).dot(nPlane);
dist = fabs(dRDotN);
if (dist > fastSkipDist){
result = SteppingHelixStateInfo::INACC;
break;
}
double curP = sv.p3.mag();
double p0 = curP;
double p0Inv = 1./p0;
Vector tau(sv.p3); tau *=p0Inv;
double tN = tau.dot(nPlane);
refDirection = tN*dRDotN < 0. ?
alongMomentum : oppositeToMomentum;
double b0 = sv.bf.mag();
if (fabs(tN)>1e-24){
tanDist = -dRDotN/tN;
} else {
tN = 1e-24;
if (fabs(dRDotN)>1e-24) tanDist = 1e6;
else tanDist = 1;
}
if (fabs(tanDist) > 1e4) tanDist = 1e4;
if (b0>1.5e-6){
double b0Inv = 1./b0;
double tNInv = 1./tN;
Vector bHat(sv.bf); bHat *=b0Inv;
double bHatN = bHat.dot(nPlane);
double cosPB = bHat.dot(tau);
double kVal = 2.99792458e-3*sv.q*p0Inv*b0;
double aVal = tanDist*kVal;
Vector lVec = bHat.cross(tau);
double bVal = lVec.dot(nPlane)*tNInv;
if (fabs(aVal*bVal)< 0.3){
double cVal = 1. - bHatN*cosPB*tNInv; // - sv.bf.cross(lVec).dot(nPlane)*b0Inv*tNInv; //1- bHat_n*bHat_tau/tau_n;
double aacVal = cVal*aVal*aVal;
if (fabs(aacVal)<1){
double tanDCorr = bVal/2. + (bVal*bVal/2. + cVal/6)*aVal;
tanDCorr *= aVal;
//+ (-bVal/24. + 0.625*bVal*bVal*bVal + 5./12.*bVal*cVal)*aVal*aVal*aVal
if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<tanDist<<" vs "
<<tanDist*(1.+tanDCorr)<<" corr "<<tanDist*tanDCorr<<std::endl;
tanDist *= (1.+tanDCorr);
} else {
if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"AACVal "<< fabs(aacVal)
<<" = "<<aVal<<"**2 * "<<cVal<<" too large:: will not converge"<<std::endl;
}
} else {
if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"ABVal "<< fabs(aVal*bVal)
<<" = "<<aVal<<" * "<<bVal<<" too large:: will not converge"<<std::endl;
}
}
result = SteppingHelixStateInfo::OK;
}
break;
case CONE_DT:
{
Point cVertex(pars[0], pars[1], pars[2]);
if (cVertex.perp2() < 1e-10){
//assumes the cone axis/vertex is along z
Vector relV3 = sv.r3 - cVertex;
double relV3mag = relV3.mag();
// double relV3Theta = relV3.theta();
double theta(pars[3]);
// double dTheta = theta-relV3Theta;
double sinTheta = sin(theta);
double cosTheta = cos(theta);
double cosV3Theta = relV3.z()/relV3mag;
if (cosV3Theta>1) cosV3Theta=1;
if (cosV3Theta<-1) cosV3Theta=-1;
double sinV3Theta = sqrt(1.-cosV3Theta*cosV3Theta);
double sinDTheta = sinTheta*cosV3Theta - cosTheta*sinV3Theta;//sin(dTheta);
double cosDTheta = cosTheta*cosV3Theta + sinTheta*sinV3Theta;//cos(dTheta);
bool isInside = sinTheta > sinV3Theta && cosTheta*cosV3Theta > 0;
dist = isInside || cosDTheta > 0 ? relV3mag*sinDTheta : relV3mag;
if (fabs(dist) > fastSkipDist){
result = SteppingHelixStateInfo::INACC;
break;
}
double relV3phi=relV3.phi();
double normPhi = isInside ?
Geom::pi() + relV3phi : relV3phi;
double normTheta = theta > Geom::pi()/2. ?
( isInside ? 1.5*Geom::pi() - theta : theta - Geom::pi()/2. )
: ( isInside ? Geom::pi()/2. - theta : theta + Geom::pi()/2 );
//this is a normVector from the cone to the point
Vector norm; norm.setRThetaPhi(fabs(dist), normTheta, normPhi);
double curP = sv.p3.mag(); double cosp3theta = sv.p3.z()/curP;
if (cosp3theta>1) cosp3theta=1;
if (cosp3theta<-1) cosp3theta=-1;
double sineConeP = sinTheta*cosp3theta - cosTheta*sqrt(1.-cosp3theta*cosp3theta);
double sinSolid = norm.dot(sv.p3)/(fabs(dist)*curP);
tanDist = fabs(sinSolid) > fabs(sineConeP) ? dist/fabs(sinSolid) : dist/fabs(sineConeP);
refDirection = sinSolid > 0 ?
oppositeToMomentum : alongMomentum;
if (debug_){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Cone pars: theta "<<theta
<<" normTheta "<<norm.theta()
<<" p3Theta "<<sv.p3.theta()
<<" sinD: "<< sineConeP
<<" sinSolid "<<sinSolid;
}
if (debug_){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"refToDest:toCone the point is "
<<(isInside? "in" : "out")<<"side the cone"
<<std::endl;
}
}
result = SteppingHelixStateInfo::OK;
}
break;
// case CYLINDER_DT:
// break;
case PATHL_DT:
{
double curS = fabs(sv.path());
dist = pars[PATHL_P] - curS;
if (fabs(dist) > fastSkipDist){
result = SteppingHelixStateInfo::INACC;
break;
}
refDirection = pars[PATHL_P] > 0 ?
alongMomentum : oppositeToMomentum;
tanDist = dist;
result = SteppingHelixStateInfo::OK;
}
break;
case POINT_PCA_DT:
{
Point pDest(pars[0], pars[1], pars[2]);
double curP = sv.p3.mag();
dist = (sv.r3 - pDest).mag()+ 1e-24;//add a small number to avoid 1/0
tanDist = (sv.r3.dot(sv.p3) - pDest.dot(sv.p3))/curP;
//account for bending in magnetic field (quite approximate)
double b0 = sv.bf.mag();
if (b0>1.5e-6){
double p0 = curP;
double kVal = 2.99792458e-3*sv.q/p0*b0;
double aVal = fabs(dist*kVal);
tanDist *= 1./(1.+ aVal);
if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"corrected by aVal "<<aVal<<" to "<<tanDist;
}
refDirection = tanDist < 0 ?
alongMomentum : oppositeToMomentum;
result = SteppingHelixStateInfo::OK;
}
break;
case LINE_PCA_DT:
{
Point rLine(pars[0], pars[1], pars[2]);
Vector dLine(pars[3], pars[4], pars[5]);
dLine = (dLine - rLine);
dLine *= 1./dLine.mag(); if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"dLine "<<dLine;
Vector dR = sv.r3 - rLine;
double curP = sv.p3.mag();
Vector dRPerp = dR - dLine*(dR.dot(dLine));
if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"dRperp "<<dRPerp;
dist = dRPerp.mag() + 1e-24;//add a small number to avoid 1/0
tanDist = dRPerp.dot(sv.p3)/curP;
//angle wrt line
double cosAlpha2 = dLine.dot(sv.p3)/curP; cosAlpha2 *= cosAlpha2;
tanDist *= 1./sqrt(fabs(1.-cosAlpha2)+1e-96);
//correct for dPhi in magnetic field: this isn't made quite right here
//(the angle between the line and the trajectory plane is neglected .. conservative)
double b0 = sv.bf.mag();
if (b0>1.5e-6){
double p0 = curP;
double kVal = 2.99792458e-3*sv.q/p0*b0;
double aVal = fabs(dist*kVal);
tanDist *= 1./(1.+ aVal);
if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"corrected by aVal "<<aVal<<" to "<<tanDist;
}
refDirection = tanDist < 0 ?
alongMomentum : oppositeToMomentum;
result = SteppingHelixStateInfo::OK;
}
break;
default:
{
//some large number
dist = 1e12;
tanDist = 1e12;
refDirection = anyDirection;
result = SteppingHelixStateInfo::NOT_IMPLEMENTED;
}
break;
}
double tanDistConstrained = tanDist;
if (fabs(tanDist) > 4.*fabs(dist) ) tanDistConstrained *= tanDist == 0 ? 0 : fabs(dist/tanDist*4.);
if (debug_){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"refToDest input: dest"<<dest<<" pars[]: ";
for (int i = 0; i < 6; i++){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<", "<<i<<" "<<pars[i];
}
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<std::endl;
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"refToDest output: "
<<"\t dist" << dist
<<"\t tanDist "<< tanDist
<<"\t tanDistConstr "<< tanDistConstrained
<<"\t refDirection "<< refDirection
<<std::endl;
}
tanDist = tanDistConstrained;
return result;
}
| 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 1960 of file SteppingHelixPropagator.cc.
References alongMomentum, anyDirection, CONE_DT, debug_, MagVolume::faces(), i, SteppingHelixStateInfo::INACC, MagVolume::inside(), VolumeBasedMagneticField::isZSymmetric(), j, LogTrace, SteppingHelixStateInfo::magVol, metname, mergeVDriftHistosByStation::name, SteppingHelixStateInfo::NOT_IMPLEMENTED, SteppingHelixStateInfo::OK, Cone::openingAngle(), SteppingHelixStateInfo::p3, pi, PLANE_DT, GloballyPositioned< T >::position(), SteppingHelixStateInfo::r3, Cylinder::radius(), RADIUS_DT, RADIUS_P, refToDest(), query::result, GloballyPositioned< T >::rotation(), MagVolume::shapeType(), pftools::UNDEFINED, UNDEFINED_DT, vbField_, Cone::vertex(), SteppingHelixStateInfo::WRONG_VOLUME, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().
Referenced by propagate().
{
static const std::string metname = "SteppingHelixPropagator";
Result result = SteppingHelixStateInfo::NOT_IMPLEMENTED;
const MagVolume* cVol = sv.magVol;
if (cVol == 0) return result;
const std::vector<VolumeSide>& cVolFaces(cVol->faces());
double distToFace[6] = {0,0,0,0,0,0};
double tanDistToFace[6] = {0,0,0,0,0,0};
PropagationDirection refDirectionToFace[6] = {anyDirection, anyDirection, anyDirection, anyDirection, anyDirection, anyDirection};
Result resultToFace[6] = {result, result, result, result, result, result};
int iFDest = -1;
int iDistMin = -1;
unsigned int iFDestSorted[6] = {0,0,0,0,0,0};
unsigned int nDestSorted =0;
unsigned int nearParallels = 0;
double curP = sv.p3.mag();
if (debug_){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Trying volume "<<DDSolidShapesName::name(cVol->shapeType())
<<" with "<<cVolFaces.size()<<" faces"<<std::endl;
}
unsigned int nFaces = cVolFaces.size();
for (unsigned int iFace = 0; iFace < nFaces; ++iFace){
if (iFace > 5){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Too many faces"<<std::endl;
}
if (debug_){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Start with face "<<iFace<<std::endl;
}
// const Plane* cPlane = dynamic_cast<const Plane*>(&cVolFaces[iFace].surface());
// const Cylinder* cCyl = dynamic_cast<const Cylinder*>(&cVolFaces[iFace].surface());
// const Cone* cCone = dynamic_cast<const Cone*>(&cVolFaces[iFace].surface());
const Surface* cPlane = 0; //only need to know the loc->glob transform
const Cylinder* cCyl = 0;
const Cone* cCone = 0;
if (typeid(cVolFaces[iFace].surface()) == typeid(const Plane&)){
cPlane = &cVolFaces[iFace].surface();
} else if (typeid(cVolFaces[iFace].surface()) == typeid(const Cylinder&)){
cCyl = dynamic_cast<const Cylinder*>(&cVolFaces[iFace].surface());
} else if (typeid(cVolFaces[iFace].surface()) == typeid(const Cone&)){
cCone = dynamic_cast<const Cone*>(&cVolFaces[iFace].surface());
} else {
edm::LogWarning(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Could not cast a volume side surface to a known type"<<std::endl;
}
if (debug_){
if (cPlane!=0) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"The face is a plane at "<<cPlane<<std::endl;
if (cCyl!=0) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"The face is a cylinder at "<<cCyl<<std::endl;
}
double pars[6] = {0,0,0,0,0,0};
DestType dType = UNDEFINED_DT;
if (cPlane != 0){
Point rPlane(cPlane->position().x(),cPlane->position().y(),cPlane->position().z());
// = cPlane->toGlobal(LocalVector(0,0,1.)); nPlane = nPlane.unit();
Vector nPlane(cPlane->rotation().zx(), cPlane->rotation().zy(), cPlane->rotation().zz()); nPlane /= nPlane.mag();
if (sv.r3.z() < 0 || !vbField_->isZSymmetric() ){
pars[0] = rPlane.x(); pars[1] = rPlane.y(); pars[2] = rPlane.z();
pars[3] = nPlane.x(); pars[4] = nPlane.y(); pars[5] = nPlane.z();
} else {
pars[0] = rPlane.x(); pars[1] = rPlane.y(); pars[2] = -rPlane.z();
pars[3] = nPlane.x(); pars[4] = nPlane.y(); pars[5] = -nPlane.z();
}
dType = PLANE_DT;
} else if (cCyl != 0){
if (debug_){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Cylinder at "<<cCyl->position()
<<" rorated by "<<cCyl->rotation()
<<std::endl;
}
pars[RADIUS_P] = cCyl->radius();
dType = RADIUS_DT;
} else if (cCone != 0){
if (debug_){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Cone at "<<cCone->position()
<<" rorated by "<<cCone->rotation()
<<" vertex at "<<cCone->vertex()
<<" angle of "<<cCone->openingAngle()
<<std::endl;
}
if (sv.r3.z() < 0 || !vbField_->isZSymmetric()){
pars[0] = cCone->vertex().x(); pars[1] = cCone->vertex().y();
pars[2] = cCone->vertex().z();
pars[3] = cCone->openingAngle();
} else {
pars[0] = cCone->vertex().x(); pars[1] = cCone->vertex().y();
pars[2] = -cCone->vertex().z();
pars[3] = Geom::pi() - cCone->openingAngle();
}
dType = CONE_DT;
} else {
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Unknown surface"<<std::endl;
resultToFace[iFace] = SteppingHelixStateInfo::UNDEFINED;
continue;
}
resultToFace[iFace] =
refToDest(dType, sv, pars,
distToFace[iFace], tanDistToFace[iFace], refDirectionToFace[iFace], fastSkipDist);
if (resultToFace[iFace] != SteppingHelixStateInfo::OK){
if (resultToFace[iFace] == SteppingHelixStateInfo::INACC) result = SteppingHelixStateInfo::INACC;
}
//keep those in right direction for later use
if (resultToFace[iFace] == SteppingHelixStateInfo::OK){
double invDTFPosiF = 1./(1e-32+fabs(tanDistToFace[iFace]));
double dSlope = fabs(distToFace[iFace])*invDTFPosiF;
double maxStepL = maxStep> 100 ? 100 : maxStep; if (maxStepL < 10) maxStepL = 10.;
bool isNearParallel = fabs(tanDistToFace[iFace]) + 100.*curP*dSlope < maxStepL //
//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
&& dSlope < 0.15 ; //
if (refDirectionToFace[iFace] == dir || isNearParallel){
if (isNearParallel) nearParallels++;
iFDestSorted[nDestSorted] = iFace;
nDestSorted++;
}
}
//pick a shortest distance here (right dir only for now)
if (refDirectionToFace[iFace] == dir){
if (iDistMin == -1) iDistMin = iFace;
else if (fabs(distToFace[iFace]) < fabs(distToFace[iDistMin])) iDistMin = iFace;
}
if (debug_)
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<cVol<<" "<<iFace<<" "
<<tanDistToFace[iFace]<<" "<<distToFace[iFace]<<" "<<refDirectionToFace[iFace]<<" "<<dir<<std::endl;
}
for (unsigned int i = 0;i<nDestSorted; ++i){
unsigned int iMax = nDestSorted-i-1;
for (unsigned int j=0;j<nDestSorted-i; ++j){
if (fabs(tanDistToFace[iFDestSorted[j]]) > fabs(tanDistToFace[iFDestSorted[iMax]]) ){
iMax = j;
}
}
unsigned int iTmp = iFDestSorted[nDestSorted-i-1];
iFDestSorted[nDestSorted-i-1] = iFDestSorted[iMax];
iFDestSorted[iMax] = iTmp;
}
if (debug_){
for (unsigned int i=0;i<nDestSorted;++i){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<cVol<<" "<<i<<" "<<iFDestSorted[i]<<" "<<tanDistToFace[iFDestSorted[i]]<<std::endl;
}
}
//now go from the shortest to the largest distance hoping to get a point in the volume.
//other than in case of a near-parallel travel this should stop after the first try
for (unsigned int i=0; i<nDestSorted;++i){
iFDest = iFDestSorted[i];
double sign = dir == alongMomentum ? 1. : -1.;
Point gPointEst(sv.r3);
Vector lDelta(sv.p3); lDelta *= sign/curP*fabs(distToFace[iFDest]);
gPointEst += lDelta;
if (debug_){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Linear est point "<<gPointEst
<<" for iFace "<<iFDest<<std::endl;
}
GlobalPoint gPointEstNegZ(gPointEst.x(), gPointEst.y(), -fabs(gPointEst.z()));
GlobalPoint gPointEstNorZ(gPointEst.x(), gPointEst.y(), gPointEst.z() );
if ( (vbField_->isZSymmetric() && cVol->inside(gPointEstNegZ))
|| ( !vbField_->isZSymmetric() && cVol->inside(gPointEstNorZ) ) ){
if (debug_){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"The point is inside the volume"<<std::endl;
}
result = SteppingHelixStateInfo::OK;
dist = distToFace[iFDest];
tanDist = tanDistToFace[iFDest];
if (debug_){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Got a point near closest boundary -- face "<<iFDest<<std::endl;
}
break;
}
}
if (result != SteppingHelixStateInfo::OK && expectNewMagVolume){
double sign = dir == alongMomentum ? 1. : -1.;
//check if it's a wrong volume situation
if (nDestSorted-nearParallels > 0 ) result = SteppingHelixStateInfo::WRONG_VOLUME;
else {
//get here if all faces in the corr direction were skipped
Point gPointEst(sv.r3);
double lDist = fabs(distToFace[iDistMin]);
if (lDist > fastSkipDist) lDist = fastSkipDist;
Vector lDelta(sv.p3); lDelta *= sign/curP*lDist;
gPointEst += lDelta;
if (debug_){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Linear est point to shortest dist "<<gPointEst
<<" for iFace "<<iDistMin<<" at distance "<<lDist*sign<<std::endl;
}
GlobalPoint gPointEstNegZ(gPointEst.x(), gPointEst.y(), -fabs(gPointEst.z()));
GlobalPoint gPointEstNorZ(gPointEst.x(), gPointEst.y(), gPointEst.z() );
if ( (vbField_->isZSymmetric() && cVol->inside(gPointEstNegZ))
|| ( !vbField_->isZSymmetric() && cVol->inside(gPointEstNorZ) ) ){
if (debug_){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"The point is inside the volume"<<std::endl;
}
}else {
result = SteppingHelixStateInfo::WRONG_VOLUME;
}
}
if (result == SteppingHelixStateInfo::WRONG_VOLUME){
dist = sign*0.05;
tanDist = dist*1.01;
if( debug_){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Wrong volume located: return small dist, tandist"<<std::endl;
}
}
}
if (result == SteppingHelixStateInfo::INACC){
if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"All faces are too far"<<std::endl;
} else if (result == SteppingHelixStateInfo::WRONG_VOLUME){
if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Appear to be in a wrong volume"<<std::endl;
} else if (result != SteppingHelixStateInfo::OK){
if (debug_) LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Something else went wrong"<<std::endl;
}
return result;
}
| SteppingHelixPropagator::Result SteppingHelixPropagator::refToMatVolume | ( | const SteppingHelixPropagator::StateInfo & | sv, |
| PropagationDirection | dir, | ||
| double & | dist, | ||
| double & | tanDist, | ||
| double | fastSkipDist = 1e12 |
||
| ) | const [protected] |
Definition at line 2204 of file SteppingHelixPropagator.cc.
References alongMomentum, anyDirection, CONE_DT, debug_, SteppingHelixStateInfo::INACC, LogTrace, metname, SteppingHelixStateInfo::NOT_IMPLEMENTED, SteppingHelixStateInfo::OK, SteppingHelixStateInfo::p3, pi, PLANE_DT, SteppingHelixStateInfo::r3, RADIUS_DT, RADIUS_P, refToDest(), query::result, SteppingHelixStateInfo::rzLims, funct::tan(), and UNDEFINED_DT.
Referenced by propagate().
{
static const std::string metname = "SteppingHelixPropagator";
Result result = SteppingHelixStateInfo::NOT_IMPLEMENTED;
double parLim[6] = {sv.rzLims.rMin, sv.rzLims.rMax,
sv.rzLims.zMin, sv.rzLims.zMax,
sv.rzLims.th1, sv.rzLims.th2 };
double distToFace[4] = {0,0,0,0};
double tanDistToFace[4] = {0,0,0,0};
PropagationDirection refDirectionToFace[4] = {anyDirection, anyDirection, anyDirection, anyDirection};
Result resultToFace[4] = {result, result, result, result};
int iFDest = -1;
double curP = sv.p3.mag();
for (unsigned int iFace = 0; iFace < 4; iFace++){
if (debug_){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Start with mat face "<<iFace<<std::endl;
}
double pars[6] = {0,0,0,0,0,0};
DestType dType = UNDEFINED_DT;
if (iFace > 1){
if (fabs(parLim[iFace+2])< 1e-6){//plane
if (sv.r3.z() < 0){
pars[0] = 0; pars[1] = 0; pars[2] = -parLim[iFace];
pars[3] = 0; pars[4] = 0; pars[5] = 1;
} else {
pars[0] = 0; pars[1] = 0; pars[2] = parLim[iFace];
pars[3] = 0; pars[4] = 0; pars[5] = 1;
}
dType = PLANE_DT;
} else {
if (sv.r3.z() > 0){
pars[0] = 0; pars[1] = 0;
pars[2] = parLim[iFace];
pars[3] = Geom::pi()/2. - parLim[iFace+2];
} else {
pars[0] = 0; pars[1] = 0;
pars[2] = - parLim[iFace];
pars[3] = Geom::pi()/2. + parLim[iFace+2];
}
dType = CONE_DT;
}
} else {
pars[RADIUS_P] = parLim[iFace];
dType = RADIUS_DT;
}
resultToFace[iFace] =
refToDest(dType, sv, pars,
distToFace[iFace], tanDistToFace[iFace], refDirectionToFace[iFace], fastSkipDist);
if (resultToFace[iFace] != SteppingHelixStateInfo::OK){
if (resultToFace[iFace] == SteppingHelixStateInfo::INACC) result = SteppingHelixStateInfo::INACC;
continue;
}
if (refDirectionToFace[iFace] == dir || fabs(distToFace[iFace]) < 2e-2*fabs(tanDistToFace[iFace]) ){
double sign = dir == alongMomentum ? 1. : -1.;
Point gPointEst(sv.r3);
Vector lDelta(sv.p3); lDelta *= sign*fabs(distToFace[iFace])/curP;
gPointEst += lDelta;
if (debug_){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Linear est point "<<gPointEst
<<std::endl;
}
double lZ = fabs(gPointEst.z());
double lR = gPointEst.perp();
double tan4 = parLim[4] == 0 ? 0 : tan(parLim[4]);
double tan5 = parLim[5] == 0 ? 0 : tan(parLim[5]);
if ( (lZ - parLim[2]) > lR*tan4
&& (lZ - parLim[3]) < lR*tan5
&& lR > parLim[0] && lR < parLim[1]){
if (debug_){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"The point is inside the volume"<<std::endl;
}
//OK, guessed a point still inside the volume
if (iFDest == -1){
iFDest = iFace;
} else {
if (fabs(tanDistToFace[iFDest]) > fabs(tanDistToFace[iFace])){
iFDest = iFace;
}
}
} else {
if (debug_){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"The point is NOT inside the volume"<<std::endl;
}
}
}
}
if (iFDest != -1){
result = SteppingHelixStateInfo::OK;
dist = distToFace[iFDest];
tanDist = tanDistToFace[iFDest];
if (debug_){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Got a point near closest boundary -- face "<<iFDest<<std::endl;
}
} else {
if (debug_){
LogTrace(metname)<<std::setprecision(17)<<std::setw(20)<<std::scientific<<"Failed to find a dest point inside the volume"<<std::endl;
}
}
return result;
}
| void SteppingHelixPropagator::setDebug | ( | bool | debug | ) | [inline] |
Switch debug printouts (to cout) .. very verbose.
Definition at line 152 of file SteppingHelixPropagator.h.
Referenced by SteppingHelixPropagatorESProducer::produce().
| void SteppingHelixPropagator::setEndcapShiftsInZPosNeg | ( | double | valPos, |
| double | valNeg | ||
| ) | [inline] |
set shifts in Z for endcap pieces (includes EE, HE, ME, YE)
Definition at line 195 of file SteppingHelixPropagator.h.
References ecShiftNeg_, and ecShiftPos_.
Referenced by SteppingHelixPropagatorESProducer::produce().
{
ecShiftPos_ = valPos; ecShiftNeg_ = valNeg;
}
| void SteppingHelixPropagator::setIState | ( | const SteppingHelixStateInfo & | sStart | ) | const [protected] |
(Internals) Init starting point
Definition at line 312 of file SteppingHelixPropagator.cc.
References cIndex_(), SteppingHelixStateInfo::covCurv, SteppingHelixStateInfo::hasErrorPropagated_, SteppingHelixStateInfo::isComplete, loadState(), noErrorPropagation_, nPoints_, SteppingHelixStateInfo::p3, Propagator::propagationDirection(), SteppingHelixStateInfo::q, SteppingHelixStateInfo::r3, and svBuf_.
Referenced by propagate(), and propagateWithPath().
{
nPoints_ = 0;
svBuf_[cIndex_(nPoints_)] = sStart; //do this anyways to have a fresh start
if (sStart.isComplete ) {
svBuf_[cIndex_(nPoints_)] = sStart;
nPoints_++;
} else {
loadState(svBuf_[cIndex_(nPoints_)], sStart.p3, sStart.r3, sStart.q,
propagationDirection(), sStart.covCurv);
nPoints_++;
}
svBuf_[cIndex_(0)].hasErrorPropagated_ = sStart.hasErrorPropagated_ & !noErrorPropagation_;
}
| void SteppingHelixPropagator::setMaterialMode | ( | bool | noMaterial | ) | [inline] |
Switch for material effects mode: no material effects if true.
Definition at line 155 of file SteppingHelixPropagator.h.
References noMaterialMode_.
Referenced by MuonSimHitProducer::beginRun(), TrackDetectorAssociator::init(), HTrackAssociator::init(), SteppingHelixPropagatorESProducer::produce(), AlCaHOCalibProducer::produce(), and PropagateToCal::propagate().
{ noMaterialMode_ = noMaterial;}
| void SteppingHelixPropagator::setNoErrorPropagation | ( | bool | noErrorPropagation | ) | [inline] |
Force no error propagation.
Definition at line 158 of file SteppingHelixPropagator.h.
References noErrorPropagation_.
Referenced by SteppingHelixPropagatorESProducer::produce(), and PropagateToCal::propagate().
{ 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 881 of file SteppingHelixPropagator.cc.
References SteppingHelixPropagator::Basis::lX, SteppingHelixPropagator::Basis::lY, SteppingHelixPropagator::Basis::lZ, and metsig::tau.
Referenced by makeAtomStep().
| void SteppingHelixPropagator::setReturnTangentPlane | ( | bool | val | ) | [inline] |
flag to return tangent plane for non-plane input
Definition at line 173 of file SteppingHelixPropagator.h.
References returnTangentPlane_.
Referenced by SteppingHelixPropagatorESProducer::produce().
{ returnTangentPlane_ = val;}
| void SteppingHelixPropagator::setSendLogWarning | ( | bool | val | ) | [inline] |
flag to send LogWarning on failures
Definition at line 176 of file SteppingHelixPropagator.h.
References sendLogWarning_.
Referenced by SteppingHelixPropagatorESProducer::produce().
{ sendLogWarning_ = val;}
| void SteppingHelixPropagator::setUseInTeslaFromMagField | ( | bool | val | ) | [inline] |
force getting field value from MagneticField, not the geometric one
Definition at line 192 of file SteppingHelixPropagator.h.
References useInTeslaFromMagField_.
Referenced by SteppingHelixPropagatorESProducer::produce().
{ 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 180 of file SteppingHelixPropagator.h.
References useIsYokeFlag_.
Referenced by SteppingHelixPropagatorESProducer::produce().
{ useIsYokeFlag_ = val;}
| void SteppingHelixPropagator::setUseMagVolumes | ( | bool | val | ) | [inline] |
Switch to using MagneticField Volumes .. as in VolumeBasedMagneticField.
Definition at line 167 of file SteppingHelixPropagator.h.
References useMagVolumes_.
Referenced by SteppingHelixPropagatorESProducer::produce().
{ useMagVolumes_ = val;}
| void SteppingHelixPropagator::setUseMatVolumes | ( | bool | val | ) | [inline] |
Switch to using Material Volumes .. internally defined for now.
Definition at line 170 of file SteppingHelixPropagator.h.
References useMatVolumes_.
Referenced by SteppingHelixPropagatorESProducer::produce().
{ 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 184 of file SteppingHelixPropagator.h.
References useTuningForL2Speed_.
Referenced by SteppingHelixPropagatorESProducer::produce().
{ 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 189 of file SteppingHelixPropagator.h.
References vbField_.
Referenced by SteppingHelixPropagatorESProducer::produce().
{ vbField_ = val;}
bool SteppingHelixPropagator::applyRadX0Correction_ [private] |
Definition at line 285 of file SteppingHelixPropagator.h.
Referenced by applyRadX0Correction(), makeAtomStep(), and SteppingHelixPropagator().
AlgebraicMatrix55 SteppingHelixPropagator::covCurvRot_ [mutable, private] |
Definition at line 276 of file SteppingHelixPropagator.h.
Referenced by SteppingHelixPropagator().
AlgebraicMatrix55 SteppingHelixPropagator::dCCurvTransform_ [mutable, private] |
Definition at line 277 of file SteppingHelixPropagator.h.
Referenced by makeAtomStep(), and SteppingHelixPropagator().
bool SteppingHelixPropagator::debug_ [private] |
Definition at line 282 of file SteppingHelixPropagator.h.
Referenced by getNextState(), isYokeVolume(), loadState(), makeAtomStep(), propagate(), refToDest(), refToMagVolume(), refToMatVolume(), setDebug(), and SteppingHelixPropagator().
double SteppingHelixPropagator::defaultStep_ [private] |
Definition at line 294 of file SteppingHelixPropagator.h.
Referenced by propagate(), and SteppingHelixPropagator().
double SteppingHelixPropagator::ecShiftNeg_ [private] |
Definition at line 297 of file SteppingHelixPropagator.h.
Referenced by getDeDx(), setEndcapShiftsInZPosNeg(), and SteppingHelixPropagator().
double SteppingHelixPropagator::ecShiftPos_ [private] |
Definition at line 296 of file SteppingHelixPropagator.h.
Referenced by getDeDx(), setEndcapShiftsInZPosNeg(), and SteppingHelixPropagator().
const MagneticField* SteppingHelixPropagator::field_ [private] |
Definition at line 279 of file SteppingHelixPropagator.h.
Referenced by getNextState(), loadState(), magneticField(), makeAtomStep(), propagate(), and SteppingHelixPropagator().
Definition at line 274 of file SteppingHelixPropagator.h.
Referenced by propagate().
const int SteppingHelixPropagator::MAX_POINTS = 7 [static, private] |
Definition at line 270 of file SteppingHelixPropagator.h.
Referenced by cIndex_(), and SteppingHelixPropagator().
const int SteppingHelixPropagator::MAX_STEPS = 10000 [static, private] |
Definition at line 269 of file SteppingHelixPropagator.h.
Referenced by propagate().
bool SteppingHelixPropagator::noErrorPropagation_ [private] |
Definition at line 284 of file SteppingHelixPropagator.h.
Referenced by setIState(), setNoErrorPropagation(), and SteppingHelixPropagator().
bool SteppingHelixPropagator::noMaterialMode_ [private] |
Definition at line 283 of file SteppingHelixPropagator.h.
Referenced by getDeDx(), setMaterialMode(), and SteppingHelixPropagator().
int SteppingHelixPropagator::nPoints_ [mutable, private] |
Definition at line 271 of file SteppingHelixPropagator.h.
Referenced by propagate(), and setIState().
bool SteppingHelixPropagator::returnTangentPlane_ [private] |
Definition at line 290 of file SteppingHelixPropagator.h.
Referenced by propagateWithPath(), setReturnTangentPlane(), and SteppingHelixPropagator().
bool SteppingHelixPropagator::sendLogWarning_ [private] |
Definition at line 291 of file SteppingHelixPropagator.h.
Referenced by propagate(), propagateWithPath(), setSendLogWarning(), and SteppingHelixPropagator().
StateInfo SteppingHelixPropagator::svBuf_[MAX_POINTS+1] [mutable, private] |
Definition at line 272 of file SteppingHelixPropagator.h.
Referenced by propagate(), propagateWithPath(), setIState(), and SteppingHelixPropagator().
const AlgebraicSymMatrix55 SteppingHelixPropagator::unit55_ [private] |
Definition at line 281 of file SteppingHelixPropagator.h.
Referenced by makeAtomStep(), and SteppingHelixPropagator().
bool SteppingHelixPropagator::useInTeslaFromMagField_ [private] |
Definition at line 289 of file SteppingHelixPropagator.h.
Referenced by getNextState(), loadState(), makeAtomStep(), setUseInTeslaFromMagField(), and SteppingHelixPropagator().
bool SteppingHelixPropagator::useIsYokeFlag_ [private] |
Definition at line 287 of file SteppingHelixPropagator.h.
Referenced by getDeDx(), getNextState(), loadState(), propagate(), setUseIsYokeFlag(), and SteppingHelixPropagator().
bool SteppingHelixPropagator::useMagVolumes_ [private] |
Definition at line 286 of file SteppingHelixPropagator.h.
Referenced by getDeDx(), getNextState(), loadState(), makeAtomStep(), propagate(), setUseMagVolumes(), and SteppingHelixPropagator().
bool SteppingHelixPropagator::useMatVolumes_ [private] |
Definition at line 288 of file SteppingHelixPropagator.h.
Referenced by propagate(), setUseMatVolumes(), and SteppingHelixPropagator().
bool SteppingHelixPropagator::useTuningForL2Speed_ [private] |
Definition at line 292 of file SteppingHelixPropagator.h.
Referenced by propagate(), setUseTuningForL2Speed(), and SteppingHelixPropagator().
const VolumeBasedMagneticField* SteppingHelixPropagator::vbField_ [private] |
Definition at line 280 of file SteppingHelixPropagator.h.
Referenced by getNextState(), loadState(), makeAtomStep(), refToMagVolume(), setVBFPointer(), and SteppingHelixPropagator().
1.7.3