AlignableCDCRecoHit Class Reference

This class is used to transfer CDC information to the track fit and Millepede. More...

#include <AlignableCDCRecoHit.h>

Inheritance diagram for AlignableCDCRecoHit:

Public Member Functions
	~AlignableCDCRecoHit () override
	Destructor.
AlignableCDCRecoHit *	clone () const override
	Creating a copy of this hit.
virtual std::pair< std::vector< int >, TMatrixD >	globalDerivatives (const genfit::StateOnPlane *sop) override
	Labels and derivatives of residuals (local measurement coordinates) w.r.t.
virtual TMatrixD	localDerivatives (const genfit::StateOnPlane *sop) override
	Derivatives for (local) fit parameters.
	CDCRecoHit ()
	Inherit constructors.
	CDCRecoHit (const CDCHit *cdcHit, const genfit::TrackCandHit *trackCandHit)
	Inherit constructors.
WireID	getWireID () const
	Getter for WireID object.
genfit::SharedPlanePtr	constructPlane (const genfit::StateOnPlane &state) const override
	Methods that actually interface to Genfit.
std::vector< genfit::MeasurementOnPlane * >	constructMeasurementsOnPlane (const genfit::StateOnPlane &state) const override
	build MeasurementsOnPlane
virtual const genfit::HMatrixU *	constructHMatrix (const genfit::AbsTrackRep *) const override
	construct error matrix
std::vector< double >	timeDerivativesMeasurementsOnPlane (const genfit::StateOnPlane &state) const
	Get the time derivative of the MesuredStateOnPlane (derived from the track fit).
bool	getFlyByDistanceVector (B2Vector3D &pointingVector, B2Vector3D &trackDir, const genfit::AbsTrackRep *rep=nullptr, bool usePlaneFromFit=false)
	Get the vector pointing from the wire to the fitted trajectory as well as the direction of the track in the fitted point.
void	setLeftRightResolution (int lr)
	select how to resolve the left/right ambiguity: -1: negative (left) side on vector (wire direction) x (track direction) 0: mirrors enter with same weight, DAF will decide.
bool	isLeftRightMeasurement () const override
	CDC RecoHits always have left-right ambiguity.
int	getLeftRightResolution () const override
	Getter for left/right passage flag.
const CDCHit *	getCDCHit () const
	get the pointer to the CDCHit object that was used to create this CDCRecoHit object.

Static Public Member Functions
static void	setTranslators (CDC::ADCCountTranslatorBase *const adcCountTranslator, CDC::CDCGeometryTranslatorBase *const cdcGeometryTranslator, CDC::TDCCountTranslatorBase *const tdcCountTranslator, bool useTrackTime=false, bool cosmics=false)
	Setter for the Translators.

Static Public Attributes
static bool	s_enableTrackT0LocalDerivative = true
	Static enabling(true) or disabling(false) addition of local derivative for track T0.
static bool	s_enableWireSaggingGlobalDerivative = false
	Static enabling(true) or disabling(false) addition of global derivative for wire sagging coefficient (per wire)
static bool	s_enableWireByWireAlignmentGlobalDerivatives = false
	Static enabling(true) or disabling(false) addition of global derivatives for wire-by-wire alignment.

Protected Member Functions
	ClassDefOverride (CDCRecoHit, 11)
	ROOT Macro.

Protected Attributes
WireID	m_wireID
	Wire Identifier.
const CDCHit *	m_cdcHit
	Pointer to the CDCHit used to created this CDCRecoHit.
unsigned short	m_tdcCount
	TDC Count as out of CDCHit.
unsigned short	m_adcCount
	ADC Count as out of CDCHit.
signed char	m_leftRight
	Flag showing left/right passage.

Static Protected Attributes
static std::unique_ptr< CDC::ADCCountTranslatorBase >	s_adcCountTranslator = 0
	Object for ADC Count translation.
static std::unique_ptr< CDC::CDCGeometryTranslatorBase >	s_cdcGeometryTranslator = 0
	Object for geometry translation.
static std::unique_ptr< CDC::TDCCountTranslatorBase >	s_tdcCountTranslator = 0
	Object for getting drift-length and -resolution.
static bool	s_useTrackTime = false
	Whether to use the track time or not when building the measurementOnPlane.
static bool	s_cosmics = false
	Switch to use cosmic events, or physics events from IP.

Private Member Functions
	ClassDefOverride (AlignableCDCRecoHit, 2)
	ROOT Macro.

Detailed Description

This class is used to transfer CDC information to the track fit and Millepede.

Definition at line 21 of file AlignableCDCRecoHit.h.

Constructor & Destructor Documentation

~AlignableCDCRecoHit()

~AlignableCDCRecoHit ( )

inlineoverride

Destructor.

Definition at line 35 of file AlignableCDCRecoHit.h.

{}

Member Function Documentation

CDCRecoHit() [1/2]

CDCRecoHit

(

)

Inherit constructors.

Definition at line 36 of file CDCRecoHit.cc.

  : genfit::AbsMeasurement(1),
    m_wireID(WireID()), m_cdcHit(nullptr), m_tdcCount(0), m_adcCount(0), m_leftRight(0)
{
}

CDCRecoHit() [2/2]

CDCRecoHit	(	const CDCHit *	cdcHit,
		const genfit::TrackCandHit *	trackCandHit )

Inherit constructors.

Definition at line 40 of file CDCRecoHit.cc.

  : genfit::AbsMeasurement(1), m_cdcHit(cdcHit), m_leftRight(0)
{
  if (s_adcCountTranslator == 0 || s_cdcGeometryTranslator == 0 || s_tdcCountTranslator == 0) {
    B2FATAL("Can't produce CDCRecoHits without setting of the translators.");
  }
 
  // get information from cdcHit into local variables.
  m_wireID = cdcHit->getID();
  m_tdcCount = cdcHit->getTDCCount();
  m_adcCount = cdcHit->getADCCount();
 
  // set l-r info
  const genfit::WireTrackCandHit* aTrackCandHitPtr =  dynamic_cast<const genfit::WireTrackCandHit*>(trackCandHit);
  if (aTrackCandHitPtr) {
    signed char lrInfo = aTrackCandHitPtr->getLeftRightResolution();
    B2DEBUG(250, "l/r: " << int(lrInfo));
    setLeftRightResolution(lrInfo);
  }
}

clone()

AlignableCDCRecoHit * clone ( ) const

inlineoverride

Creating a copy of this hit.

Definition at line 38 of file AlignableCDCRecoHit.h.

    {
      return new AlignableCDCRecoHit(*this);
    }

constructHMatrix()

const genfit::HMatrixU * constructHMatrix ( const genfit::AbsTrackRep * rep ) const

overridevirtualinherited

construct error matrix

Definition at line 252 of file CDCRecoHit.cc.

{
  if (!dynamic_cast<const genfit::RKTrackRep*>(rep)) {
    B2FATAL("CDCRecoHit can only handle state vectors of type RKTrackRep!");
  }
 
  return new genfit::HMatrixU();
}

constructMeasurementsOnPlane()

std::vector< genfit::MeasurementOnPlane * > constructMeasurementsOnPlane ( const genfit::StateOnPlane & state ) const

overrideinherited

build MeasurementsOnPlane

Definition at line 141 of file CDCRecoHit.cc.

{
  double z = state.getPos().Z();
  const B2Vector3D& p = state.getMom();
  // Calculate alpha and theta.  A description was given in
  // https://indico.mpp.mpg.de/getFile.py/access?contribId=5&sessionId=3&resId=0&materialId=slides&confId=3195
 
//Comment out the following 2 lines since they introduce dependence on cdclib (or circular dependence betw. cdc_objects and cdclib).
//  double alpha = CDCGeometryPar::Instance().getAlpha(state.getPlane()->getO(), p);
//  double theta = CDCGeometryPar::Instance().getTheta(p);
 
//N.B. The following 8 lines are tentative to avoid the circular dependence mentioned above ! The definitions of alpha and theta should be identical to those defined in CDCGeometryPar.
  const double wx = state.getPlane()->getO().X();
  const double wy = state.getPlane()->getO().Y();
  const double px = p.X();
  const double py = p.Y();
  const double cross = wx * py - wy * px;
  const double dot   = wx * px + wy * py;
  double alpha = atan2(cross, dot);
  double theta = atan2(p.Perp(), p.Z());
  /*
  double alpha0 =  CDCGeometryPar::Instance().getAlpha(state.getPlane()->getO(), p);
  double theta0 =  CDCGeometryPar::Instance().getTheta(p);
  if (alpha != alpha0) {
    std::cout <<"alpha,alpha0= " << alpha <<" "<< alpha0 << std::endl;
    exit(-1);
  }
  if (theta != theta0) {;
    std::cout <<"theta,theta0= " << theta <<" "<< theta0 << std::endl;
    exit(-2);
  }
  */
 
  double trackTime = s_useTrackTime ? state.getTime() : 0;
  //temp4cosmics
  //  std::cout <<"phi,trackTime= " << atan2(py,px) <<" "<< trackTime << std::endl;
  if (s_cosmics) {
    if (atan2(py, px) > 0.) {
      //    if (atan2(wy,wx) > 0.) {
      trackTime *= -1.;
    }
  }
 
  // The meaning of the left / right flag (called
  // 'ambiguityDiscriminator' in TDCCounTranslatorBase) is inferred
  // from CDCGeometryPar::getNewLeftRightRaw().
  double mL = s_tdcCountTranslator->getDriftLength(m_tdcCount, m_wireID, trackTime,
                                                   false, //left
                                                   z, alpha, theta, m_adcCount);
  double mR = s_tdcCountTranslator->getDriftLength(m_tdcCount, m_wireID, trackTime,
                                                   true, //right
                                                   z, alpha, theta, m_adcCount);
  double VL = s_tdcCountTranslator->getDriftLengthResolution(mL, m_wireID,
                                                             false, //left
                                                             z, alpha, theta);
  double VR = s_tdcCountTranslator->getDriftLengthResolution(mR, m_wireID,
                                                             true, //right
                                                             z, alpha, theta);
 
  // static to avoid constructing these over and over.
  static TVectorD m(1);
  static TMatrixDSym cov(1);
 
  m(0) = mR;
  cov(0, 0) = VR;
  auto mopR = new genfit::MeasurementOnPlane(m, cov, state.getPlane(), state.getRep(),
                                             constructHMatrix(state.getRep()));
  m(0) = -mL; // Convert from unsigned drift length to signed coordinate.
  cov(0, 0) = VL;
  auto mopL = new genfit::MeasurementOnPlane(m, cov, state.getPlane(), state.getRep(),
                                             constructHMatrix(state.getRep()));
 
  // set left/right weights
  if (m_leftRight < 0) {
    mopL->setWeight(1);
    mopR->setWeight(0);
  } else if (m_leftRight > 0) {
    mopL->setWeight(0);
    mopR->setWeight(1);
  } else {
    // In absence of L/R information set equal weight for mirror hits.
    // We have this weight decrease as the drift distance increases.
    // Since the average will always be close to the wire, the bias
    // introduced by the averaged hit would increase as the drift
    // radius increases.  Reducing the initial weight effectively
    // counteracts this.  This way the DAF can figure the ambiguities
    // out quickly.
    //
    // Wire spacing in the radial direction according to TDR is 10 mm
    // in the innermost superlayer, 18 mm otherwise.  Use these values
    // times a safety margin of 1.5 as upper bound for the drift
    // radii.  The max distance between the mirror hits is twice the
    // maximal drift radius.
    double rMax = 1.5 * (m_wireID.getISuperLayer() == 0 ? 1. : 1.8);
    double weight = 0.5 * square(std::max(0., 1 - (mR + mL) / 2 / rMax));
    mopL->setWeight(weight);
    mopR->setWeight(weight);
  }
 
  // Ignore hits with negative drift times.  For these, the
  // TDCCountTranslator returns a negative drift length.
  if (mL < 0. || mR < 0.) {
    B2DEBUG(150, "Ignoring hit with negative drift time.");
    mopL->setWeight(0);
    mopR->setWeight(0);
  }
 
  return {mopL, mopR};
}

constructPlane()

genfit::SharedPlanePtr constructPlane ( const genfit::StateOnPlane & state ) const

overrideinherited

Methods that actually interface to Genfit.

Definition at line 81 of file CDCRecoHit.cc.

{
  // We find the plane in two steps: first we neglect wire sag to get
  // a good estimate of the z coordinate of the crossing.  Then we use
  // this z coordinate to find the point of closest approach to the
  // sagging wire and its local direction.
 
  // Don't use clone: we don't want to extrapolate covariances if
  // state is a genfit::MeasuredStateOnPlane.
  genfit::StateOnPlane st(state);
 
  const B2Vector3D& noSagWire1(s_cdcGeometryTranslator->getWireBackwardPosition(m_wireID));
  const B2Vector3D& noSagWire2(s_cdcGeometryTranslator->getWireForwardPosition(m_wireID));
 
  // unit vector along the wire
  B2Vector3D noSagWireDirection = noSagWire2 - noSagWire1;
  noSagWireDirection.SetMag(1.);
 
  // point of closest approach
  const genfit::AbsTrackRep* rep = state.getRep();
  rep->extrapolateToLine(st, noSagWire1, noSagWireDirection);
  const B2Vector3D& noSagPoca = rep->getPos(st);
 
  double zPOCA = (noSagWire1.Z()
                  + noSagWireDirection.Dot(noSagPoca - noSagWire1) * noSagWireDirection.Z());
 
  // Now re-extrapolate taking Z of trajectory into account.
  const B2Vector3D& wire1(s_cdcGeometryTranslator->getWireBackwardPosition(m_wireID, zPOCA));
  const B2Vector3D& wire2(s_cdcGeometryTranslator->getWireForwardPosition(m_wireID, zPOCA));
 
  // unit vector along the wire (will become V of plane)
  B2Vector3D wireDirection = wire2 - wire1;
  wireDirection.SetMag(1.);
 
  // point of closest approach
  rep->extrapolateToLine(st, wire1, wireDirection);
  const B2Vector3D& poca = rep->getPos(st);
  B2Vector3D dirInPoca = rep->getMom(st);
  dirInPoca.SetMag(1.);
  const B2Vector3D& pocaOnWire = wire1 + wireDirection.Dot(poca - wire1) * wireDirection;
  //temp
  //  std::cout << (noSagWire1 + noSagWireDirection.Dot(noSagPoca - noSagWire1) * noSagWireDirection).Y() <<" "<< pocaOnWire.Y() <<" " << (noSagWire1 + noSagWireDirection.Dot(noSagPoca - noSagWire1) * noSagWireDirection - pocaOnWire).Y() << std::endl;
  //  std::cout << (noSagWire1 + noSagWireDirection.Dot(noSagPoca - noSagWire1) * noSagWireDirection).Perp() <<" "<< pocaOnWire.Perp() << std::endl;
  //  std::cout << (noSagWire1 + noSagWireDirection.Dot(noSagPoca - noSagWire1) * noSagWireDirection).Z() <<" "<< pocaOnWire.Z() << std::endl;
 
  // check if direction is parallel to wire
  if (fabs(wireDirection.Angle(dirInPoca)) < 0.01) {
    genfit::Exception exc("CDCRecoHit::constructPlane(): Cannot construct detector plane, direction is parallel to wire", __LINE__,
                          __FILE__);
    throw exc;
  }
 
  // construct orthogonal (unit) vector
  const B2Vector3D& U = wireDirection.Cross(dirInPoca);
 
  genfit::SharedPlanePtr pl = genfit::SharedPlanePtr(new genfit::DetPlane(pocaOnWire, U, wireDirection));
  //pl->Print();
  return pl;
}

getCDCHit()

const CDCHit * getCDCHit ( ) const

inlineinherited

get the pointer to the CDCHit object that was used to create this CDCRecoHit object.

Can be NULL if CDCRecoHit was not created with the CDCRecoHit(const CDCHit* cdcHit) constructor

Definition at line 112 of file CDCRecoHit.h.

    {
      return m_cdcHit;
    }

getFlyByDistanceVector()

bool getFlyByDistanceVector ( B2Vector3D & pointingVector, B2Vector3D & trackDir, const genfit::AbsTrackRep * rep = nullptr, bool usePlaneFromFit = false )

inherited

Get the vector pointing from the wire to the fitted trajectory as well as the direction of the track in the fitted point.

Uses the cardinal TrackRep if rep == NULL. Returns false if the track is not fitted in this point for the requested TrackRep. If the hit does not belong to a track, an error is issued and false returned. If usePlaneFromFit is used, constructPlane is not called to evaluate the closest point on the wire to the track. Instead the origin of the plane of the fitted state is used (which will be the same point if the wire is not bent).

Definition at line 369 of file CDCRecoHit.cc.

{
  const genfit::TrackPoint* tp = this->getTrackPoint();
  if (!tp) {
    B2ERROR("No genfit::TrackPoint for CDCRecoHit.");
    return false;
  }
  const genfit::AbsFitterInfo* fi = tp->getFitterInfo(rep);
  if (!fi) {
    B2DEBUG(200, "No genfit::AbsFitterInfo for this CDCRecoHit.");
    return false;
  }
 
  const genfit::MeasuredStateOnPlane& mop = fi->getFittedState();
  B2Vector3D fittedPoca = mop.getPos();
  // constructPlane places the coordinate center in the POCA to the
  // wire.  Using this is the default behavior.  If this should be too
  // slow, as it has to re-evaluate the POCA, alternatively the user
  // can set usePlaneFromFit which uses the plane determined by the
  // track fit.
  B2Vector3D pocaOnWire = (usePlaneFromFit
                           ? mop.getPlane()->getO()
                           : this->constructPlane(mop)->getO());
 
  // The vector from the wire to the track.
  pointingVector = fittedPoca - pocaOnWire;
 
  trackDir = mop.getMom();
  trackDir.SetMag(1.);
  return true;
}

getLeftRightResolution()

int getLeftRightResolution ( ) const

inlineoverrideinherited

Getter for left/right passage flag.

Definition at line 106 of file CDCRecoHit.h.

{ return m_leftRight; }

getWireID()

WireID getWireID ( ) const

inlineinherited

Getter for WireID object.

Definition at line 49 of file CDCRecoHit.h.

    {
      return m_wireID;
    }

globalDerivatives()

std::pair< std::vector< int >, TMatrixD > globalDerivatives ( const genfit::StateOnPlane * sop )

overridevirtual

Labels and derivatives of residuals (local measurement coordinates) w.r.t.

alignment/calibration parameters Matrix "G" of derivatives valid for given prediction of track state:

G(i, j) = d_residual_i/d_parameter_j

For 2D measurement (u,v):

G = ( du/da du/db du/dc ... ) ( dv/da dv/db dv/dc ... )

for calibration parameters a, b, c.

For 1D measurement:

G = ( 0 0 0 ... ) ( dv/da dv/db dv/dc ... ) for V-strip,

G = ( du/da du/db du/dc ... ) ( 0 0 0 ... ) for U-strip,

Measurements with more dimensions (slopes, curvature) should provide full 4-5Dx(n params) matrix (state as (q/p, u', v', u, v) or (u', v', u, v))

Parameters

sop	Predicted state of the track as linearization point around which derivatives of alignment/calibration parameters shall be computed

Returns

pair<vector<int>, TMatrixD> With matrix with number of rows = dimension of residual, number of columns = number of parameters. number of columns must match vector<int>.size().

Definition at line 28 of file AlignableCDCRecoHit.cc.

{
  GlobalDerivatives globals;
  unsigned short LR = (int(m_leftRight) > 0.) ? 1 : 0;
 
  const B2Vector3D& mom = sop->getMom();
  const B2Vector3D& wirePositon = sop->getPlane()->getO();
  const unsigned short layer = getWireID().getICLayer();
 
  const CDCGeometryPar& cdcgeo = CDCGeometryPar::Instance();
  double alpha = cdcgeo.getAlpha(wirePositon, mom);
  double theta = cdcgeo.getTheta(mom);
  const TVectorD& stateOnPlane = sop->getState();
  const double driftLengthEstimate = std::abs(stateOnPlane(3));
  const double driftTime = cdcgeo.getDriftTime(driftLengthEstimate, layer, LR, alpha, theta);
  const double driftVelocity = cdcgeo.getDriftV(driftTime, layer, LR, alpha, theta);
 
  // CDC Calibration -------------------------------------------------
 
  // Time zero calibration (per wire)
  if (driftTime > 20 && driftTime < 200 && fabs(driftVelocity) < 1.0e-2) {
    globals.add(
      GlobalLabel::construct<CDCTimeZeros>(getWireID(), 0),
      driftVelocity * double(int(m_leftRight))
    );
  }
 
  // Time walk calibration (per board)
  if (driftTime > 20 && driftTime < 200 && fabs(driftVelocity) < 1.0e-2 && m_adcCount < 400 && m_adcCount > 2) {
    DBObjPtr<CDCTimeWalks> tws;
    unsigned short board = CDCGeometryPar::Instance().getBoardID(getWireID());
    const std::vector<float> param  = tws->getTimeWalkParams(board);
    globals.add(
      GlobalLabel::construct<CDCTimeWalks>(board, 0),
      -driftVelocity * exp(-param[1]* m_adcCount) * double(int(m_leftRight))
    );
    globals.add(
      GlobalLabel::construct<CDCTimeWalks>(board, 1),
      driftVelocity * m_adcCount * param[0]*exp(-param[1]* m_adcCount) * double(int(m_leftRight))
    );
  }
 
  // Space time relations calibration
  if (driftTime > 20 && driftTime < 500 && fabs(driftVelocity) < 1.0e-2) {
    // TODO: ugly to need to ask XTRelations for something here...
    // Can't I get this CDCGeometryPar or sth. like this?
 
    theta = cdcgeo.getOutgoingTheta(alpha, theta);
    alpha = cdcgeo.getOutgoingAlpha(alpha);
    DBObjPtr<CDCXtRelations> xts;
    auto xtid = xts->getXtID(getWireID().getICLayer(), LR, (float)alpha, (float)theta);
    const auto& par = xts->getXtParams(xtid);
    auto boundary = par.at(6);
    if (driftTime < boundary) {
      globals.add(
        GlobalLabel::construct<CDCXtRelations>(xtid, 0),
        ROOT::Math::Chebyshev5(driftTime, 1, 0, 0, 0, 0, 0) * double(int(m_leftRight))
      );
      globals.add(
        GlobalLabel::construct<CDCXtRelations>(xtid, 1),
        ROOT::Math::Chebyshev5(driftTime, 0, 1, 0, 0, 0, 0) * double(int(m_leftRight))
      );
      globals.add(
        GlobalLabel::construct<CDCXtRelations>(xtid, 2),
        ROOT::Math::Chebyshev5(driftTime, 0, 0, 1, 0, 0, 0) * double(int(m_leftRight))
      );
      globals.add(
        GlobalLabel::construct<CDCXtRelations>(xtid, 3),
        ROOT::Math::Chebyshev5(driftTime, 0, 0, 0, 1, 0, 0) * double(int(m_leftRight))
      );
      globals.add(
        GlobalLabel::construct<CDCXtRelations>(xtid, 4),
        ROOT::Math::Chebyshev5(driftTime, 0, 0, 0, 0, 1, 0) * double(int(m_leftRight))
      );
      globals.add(
        GlobalLabel::construct<CDCXtRelations>(xtid, 5),
        ROOT::Math::Chebyshev5(driftTime, 0, 0, 0, 0, 0, 1) * double(int(m_leftRight))
      );
    } else {
      globals.add(
        GlobalLabel::construct<CDCXtRelations>(xtid, 7),
        double(int(m_leftRight))
      );
      globals.add(
        GlobalLabel::construct<CDCXtRelations>(xtid, 8),
        (driftTime - boundary) * double(int(m_leftRight))
      );
    }
  }
 
 
  // CDC Alignment ---------------------------------------------------
 
  auto tdir = sop->getDir();
  auto ndir = sop->getPlane()->getNormal();
  auto udir = sop->getPlane()->getU();
  auto vdir = sop->getPlane()->getV();
  auto pos  = sop->getPos();
 
  auto tn = tdir[0] * ndir[0] + tdir[1] * ndir[1] + tdir[2] * ndir [2];
  // d residual / d measurement
  auto drdm = TMatrixD(3, 3);
  drdm.UnitMatrix();
  for (int i = 0; i < 3; ++i) {
    for (int j = 0; j < 3; ++j) {
      drdm(i, j) -= tdir[i] * ndir[j] / tn;
    }
  }
  // d measurement / d global rigid body param
  auto dmdg = TMatrixD(3, 6);
  dmdg.Zero();
  dmdg[0][0] = 1.; dmdg[0][4] = -pos[2]; dmdg[0][5] =  pos[1];
  dmdg[1][1] = 1.; dmdg[1][3] =  pos[2]; dmdg[1][5] = -pos[0];
  dmdg[2][2] = 1.; dmdg[2][3] = -pos[1]; dmdg[2][4] =  pos[0];
  // d local residual / d global residual
  auto drldrg = TMatrixD(3, 3);
  for (int i = 0; i < 3; ++i) {
    drldrg(0, i) = udir[i];
    drldrg(1, i) = vdir[i];
    drldrg(2, i) = ndir[i];
  }
  // d local residual / d global rigid body param
  auto drldg = drldrg * (drdm * dmdg);
 
  // wire ends
  const double zWireM = s_cdcGeometryTranslator->getWireBackwardPosition(getWireID(), CDCGeometryPar::c_Aligned)[2];
  const double zWireP = s_cdcGeometryTranslator->getWireForwardPosition(getWireID(), CDCGeometryPar::c_Aligned)[2];
  // relative Z position [0..1]
  const double zRel = std::max(0., std::min(1., (pos[2] - zWireM) / (zWireP - zWireM)));
 
  // Layer alignment
  // wire 511 = no wire (0 is reserved for existing wires) - this has to be compatible with code in CDCGeometryPar::setWirPosAlignParams
  auto layerID = WireID(getWireID().getICLayer(), 511);
 
  // Alignment of layer X (bwd)
  globals.add(
    GlobalLabel::construct<CDCAlignment>(layerID, CDCAlignment::layerX),
    drldg(0, 0)
  );
 
  // Alignment of layer Y (bwd)
  globals.add(
    GlobalLabel::construct<CDCAlignment>(layerID, CDCAlignment::layerY),
    drldg(0, 1)
  );
 
  // Alignment of layer rotation (gamma) (bwd)
  globals.add(
    GlobalLabel::construct<CDCAlignment>(layerID, CDCAlignment::layerPhi),
    drldg(0, 5)
  );
 
  // Difference between wire ends (end plates)
  // Alignment of layer dX, dX = forward - backward endplate
  globals.add(
    GlobalLabel::construct<CDCAlignment>(layerID, CDCAlignment::layerDx),
    drldg(0, 0) * zRel
  );
 
  // Alignment of layer dY, dY = forward - backward endplate
  globals.add(
    GlobalLabel::construct<CDCAlignment>(layerID, CDCAlignment::layerDy),
    drldg(0, 1) * zRel
  );
 
  // Alignment of layer rotation difference d(gamma or phi), dPhi = forward - backward endplate
  globals.add(
    GlobalLabel::construct<CDCAlignment>(layerID, CDCAlignment::layerDPhi),
    drldg(0, 5) * zRel
  );
 
  //WARNING: experimental (disabled by default)
  // Wire-by-wire alignment
  if (s_enableWireByWireAlignmentGlobalDerivatives) {
    // How much shift (in X or Y) on BWD wire-end will change the residual at estimated track crossing
    // with the wire (at relative z-position on wire = zRel)
    double zRelM = fabs(1. - zRel);
    // Same as above but for FWD wire-end (residual at zRel = zRel * delta(X or Y at FWD wire-end)
    double zRelP = fabs(zRel - 0.);
 
    // Alignment of wires X in global coords at BWD wire-end
    globals.add(
      GlobalLabel::construct<CDCAlignment>(getWireID().getEWire(), CDCAlignment::wireBwdX),
      drldg(0, 0) * zRelM
    );
    // Alignment of wires X in global coords at FWD wire-end
    globals.add(
      GlobalLabel::construct<CDCAlignment>(getWireID().getEWire(), CDCAlignment::wireFwdX),
      drldg(0, 0) * zRelP
    );
 
    // Alignment of wires Y in global coords at BWD wire-end
    globals.add(
      GlobalLabel::construct<CDCAlignment>(getWireID().getEWire(), CDCAlignment::wireBwdY),
      drldg(0, 1) * zRelM
    );
    // Alignment of wires Y in global coords at FWD wire-end
    globals.add(
      GlobalLabel::construct<CDCAlignment>(getWireID().getEWire(), CDCAlignment::wireFwdY),
      drldg(0, 1) * zRelP
    );
  }
 
  //WARNING: experimental (disabled by default)
  //TODO: missing some factors (we do not align directly the wire-sag coefficient, but
  // wire tension ... coef = pi * ro * r * r / 8 / tension
  //TODO: need to get these numbers from CDCGeometryPar!
  if (s_enableWireSaggingGlobalDerivative) {
    globals.add(
      GlobalLabel::construct<CDCAlignment>(getWireID().getEWire(), CDCAlignment::wireTension),
      drldg(0, 1) * 4.0 * zRel * (1.0 - zRel)
    );
  }
 
 
 
  return globals;
}

isLeftRightMeasurement()

bool isLeftRightMeasurement ( ) const

inlineoverrideinherited

CDC RecoHits always have left-right ambiguity.

Definition at line 103 of file CDCRecoHit.h.

{ return true; }

localDerivatives()

TMatrixD localDerivatives ( const genfit::StateOnPlane * sop )

overridevirtual

Derivatives for (local) fit parameters.

Parameters

sop	State on virtual plane to calculate derivatives

Returns

TMatrixD of local derivatives, number of columns = number of params, number of rows = 2 (or measurement dimension if > 2)

Definition at line 248 of file AlignableCDCRecoHit.cc.

{
  if (!s_enableTrackT0LocalDerivative)
    return TMatrixD();
 
  unsigned short LR = (int(m_leftRight) > 0.) ? 1 : 0;
 
  const B2Vector3D& mom = sop->getMom();
  const B2Vector3D& wirePositon = sop->getPlane()->getO();
  const unsigned short layer = getWireID().getICLayer();
 
  const CDCGeometryPar& cdcgeo = CDCGeometryPar::Instance();
  const double alpha = cdcgeo.getAlpha(wirePositon, mom);
  const double theta = cdcgeo.getTheta(mom);
  const TVectorD& stateOnPlane = sop->getState();
  const double driftLengthEstimate = std::abs(stateOnPlane(3));
  const double driftTime = cdcgeo.getDriftTime(driftLengthEstimate, layer, LR, alpha, theta);
  const double driftVelocity = cdcgeo.getDriftV(driftTime, layer, LR, alpha, theta);
 
  TMatrixD locals(2, 1);
  if (driftTime > 20 && driftTime < 200 && fabs(driftVelocity) < 1.0e-2) {
    locals(0, 0) = - double(int(m_leftRight)) * driftVelocity;
    locals(1, 0) = 0.; // insensitive coordinate along wire
  }
 
  return locals;
}

setLeftRightResolution()

void setLeftRightResolution ( int lr )

inlineinherited

select how to resolve the left/right ambiguity: -1: negative (left) side on vector (wire direction) x (track direction) 0: mirrors enter with same weight, DAF will decide.

1: positive (right) side on vector (wire direction) x (track direction) where the wire direction is pointing towards +z except for small corrections such as stereo angle, sagging

Definition at line 100 of file CDCRecoHit.h.

{ m_leftRight = lr; }

setTranslators()

void setTranslators ( CDC::ADCCountTranslatorBase *const adcCountTranslator, CDC::CDCGeometryTranslatorBase *const cdcGeometryTranslator, CDC::TDCCountTranslatorBase *const tdcCountTranslator, bool useTrackTime = false, bool cosmics = false )

staticinherited

Setter for the Translators.

Definition at line 34 of file CDCRecoHit.cc.

{
  s_adcCountTranslator.reset(adcCountTranslator);
  s_cdcGeometryTranslator.reset(cdcGeometryTranslator);
  s_tdcCountTranslator.reset(tdcCountTranslator);
  s_useTrackTime = useTrackTime;
  //temp4cosmics
  s_cosmics = cosmics;
}

timeDerivativesMeasurementsOnPlane()

std::vector< double > timeDerivativesMeasurementsOnPlane ( const genfit::StateOnPlane & state ) const

inherited

Get the time derivative of the MesuredStateOnPlane (derived from the track fit).

Definition at line 262 of file CDCRecoHit.cc.

{
  double z = state.getPos().Z();
  const B2Vector3D& p = state.getMom();
  // Calculate alpha and theta.  A description was given by in
  // https://indico.mpp.mpg.de/getFile.py/access?contribId=5&sessionId=3&resId=0&materialId=slides&confId=3195
 
//Comment out the following 2 lines since they introduce dependence on cdclib (or circular dependence betw. cdc_objects and cdclib).
//  double alpha = CDCGeometryPar::Instance().getAlpha(state.getPlane()->getO(), p);
//  double theta = CDCGeometryPar::Instance().getTheta(p);
 
//N.B. The following 8 lines are tentative to avoid the circular dependence mentioned above ! The definitions of alpha and theta should be identical to those defined in CDCGeometryPar.
  const double wx = state.getPlane()->getO().X();
  const double wy = state.getPlane()->getO().Y();
  const double px = p.X();
  const double py = p.Y();
  const double cross = wx * py - wy * px;
  const double dot   = wx * px + wy * py;
  double alpha = atan2(cross, dot);
  double theta = atan2(p.Perp(), p.Z());
  /*
  double alpha0 =  CDCGeometryPar::Instance().getAlpha(state.getPlane()->getO(), p);
  double theta0 =  CDCGeometryPar::Instance().getTheta(p);
  if (alpha != alpha0) {
    std::cout <<"alpha,alpha0= " << alpha <<" "<< alpha0 << std::endl;
    exit(-1);
  }
  if (theta != theta0) {;
    std::cout <<"theta,theta0= " << theta <<" "<< theta0 << std::endl;
    exit(-2);
  }
  */
 
  double trackTime = s_useTrackTime ? state.getTime() : 0;
 
  // The meaning of the left / right flag (called
  // 'ambiguityDiscriminator' in TDCCounTranslatorBase) is inferred
  // from CDCGeometryPar::getNewLeftRightRaw().
  auto fL = [&](const double & t) -> double {
    return s_tdcCountTranslator->getDriftLength(m_tdcCount, m_wireID, t,
                                                false, //left
                                                z, alpha, theta, m_adcCount); };
  auto fR = [&](const double & t) -> double {
    return s_tdcCountTranslator->getDriftLength(m_tdcCount, m_wireID, t,
                                                true, //right
                                                z, alpha, theta, m_adcCount); };
 
  // Calculate derivative for all left and right mirror hit.
  //
  // It's a polynomial, but let's not meddle with the innard of the
  // CDC code for now.
  //
  // The algorithm follows the one in TF1::Derivative() :
  //   df(x) = (4 D(h/2) - D(h)) / 3
  // with D(h) = (f(x + h) - f(x - h)) / (2 h).
  double rightShort[2], rightFull[2];
  double leftShort[2], leftFull[2];
  const double defaultStepT = 1e-3 * trackTime;
  double stepT;
  {
    double temp = trackTime + defaultStepT / 2;
    // Find the actual size of the step, which will differ from
    // defaultStepX due to roundoff.  This is the step-size we will
    // use for this direction.  Idea taken from Numerical Recipes,
    // 3rd ed., section 5.7.
    //
    // Note that if a number is exactly representable, it's double
    // will also be exact.  Outside denormals, this also holds for
    // halving.  Unless the exponent changes (which it only will in
    // the vicinity of zero) adding or subtracing doesn't make a
    // difference.
    //
    // We determine the roundoff error for the half-step.  If this
    // is exactly representable, the full step will also be.
    stepT = 2 * (temp - trackTime);
 
    rightShort[0] = fL(trackTime + .5 * stepT);
    rightShort[1] = fR(trackTime + .5 * stepT);
  }
  {
    leftShort[0] = fL(trackTime - .5 * stepT);
    leftShort[1] = fR(trackTime - .5 * stepT);
  }
  {
    rightFull[0] = fL(trackTime + stepT);
    rightFull[1] = fR(trackTime + stepT);
  }
  {
    leftFull[0] = fL(trackTime - stepT);
    leftFull[1] = fR(trackTime - stepT);
  }
 
  // Calculate the derivatives for the individual components of
  // the track parameters.
  double derivFull[2];
  double derivShort[2];
  for (size_t j = 0; j < 2; ++j) {
    derivFull[j] = (rightFull[j] - leftFull[j]) / (2.*stepT);
    derivShort[j] = (rightShort[j] - leftShort[j]) / stepT;
  }
  //std::cout << rightShort[0] << " " << derivShort[0] << " " << trackTime << std::endl;
  return { +(4.*derivShort[0] - derivFull[0]) / 3.,
           -(4.*derivShort[1] - derivFull[1]) / 3.};
}

Member Data Documentation

m_adcCount

unsigned short m_adcCount

protectedinherited

ADC Count as out of CDCHit.

Definition at line 150 of file CDCRecoHit.h.

m_cdcHit

const CDCHit* m_cdcHit

protectedinherited

Pointer to the CDCHit used to created this CDCRecoHit.

Definition at line 144 of file CDCRecoHit.h.

m_leftRight

signed char m_leftRight

protectedinherited

Flag showing left/right passage.

Definition at line 153 of file CDCRecoHit.h.

m_tdcCount

unsigned short m_tdcCount

protectedinherited

TDC Count as out of CDCHit.

Definition at line 147 of file CDCRecoHit.h.

m_wireID

WireID m_wireID

protectedinherited

Wire Identifier.

Definition at line 141 of file CDCRecoHit.h.

s_adcCountTranslator

std::unique_ptr< ADCCountTranslatorBase > s_adcCountTranslator = 0

staticprotectedinherited

Object for ADC Count translation.

Definition at line 120 of file CDCRecoHit.h.

s_cdcGeometryTranslator

std::unique_ptr< CDCGeometryTranslatorBase > s_cdcGeometryTranslator = 0

staticprotectedinherited

Object for geometry translation.

Definition at line 123 of file CDCRecoHit.h.

s_cosmics

bool s_cosmics = false

staticprotectedinherited

Switch to use cosmic events, or physics events from IP.

Default value is 0, which means "physics event" mode.

Definition at line 136 of file CDCRecoHit.h.

s_enableTrackT0LocalDerivative

bool s_enableTrackT0LocalDerivative = true

static

Static enabling(true) or disabling(false) addition of local derivative for track T0.

Definition at line 25 of file AlignableCDCRecoHit.h.

s_enableWireByWireAlignmentGlobalDerivatives

bool s_enableWireByWireAlignmentGlobalDerivatives = false

static

Static enabling(true) or disabling(false) addition of global derivatives for wire-by-wire alignment.

Definition at line 29 of file AlignableCDCRecoHit.h.

s_enableWireSaggingGlobalDerivative

bool s_enableWireSaggingGlobalDerivative = false

static

Static enabling(true) or disabling(false) addition of global derivative for wire sagging coefficient (per wire)

Definition at line 27 of file AlignableCDCRecoHit.h.

s_tdcCountTranslator

std::unique_ptr< TDCCountTranslatorBase > s_tdcCountTranslator = 0

staticprotectedinherited

Object for getting drift-length and -resolution.

Definition at line 126 of file CDCRecoHit.h.

s_useTrackTime

bool s_useTrackTime = false

staticprotectedinherited

Whether to use the track time or not when building the measurementOnPlane.

This needs to be in sync with the TDCCountTranslator.

Definition at line 131 of file CDCRecoHit.h.

---

The documentation for this class was generated from the following files:

• alignment/reconstruction/include/AlignableCDCRecoHit.h
• alignment/reconstruction/src/AlignableCDCRecoHit.cc