Add longitudinal position residual

Examples/DOCAWireHitUpdater.hh

@@ -7,25 +7,23 @@ namespace KinKal {
     public:
       DOCAWireHitUpdater(double mindoca,double maxdoca ) : mindoca_(mindoca), maxdoca_(maxdoca) {}
       // define the state given the (presumably unbiased) distance of closest approach
-      WireHitState wireHitState(double doca) const;
+      void updateWireHitState(WireHitState &state, double doca) const;
       double minDOCA() const { return mindoca_; }
       double maxDOCA() const { return maxdoca_; }
     private:
       double mindoca_; // minimum DOCA value to sign LR ambiguity
       double maxdoca_; // maximum DOCA to still use a hit
   };
 
-  WireHitState DOCAWireHitUpdater::wireHitState(double doca) const {
-    WireHitState state;
+  void DOCAWireHitUpdater::updateWireHitState(WireHitState &state, double doca) const {
     if(fabs(doca) > maxdoca_ ) {
-      state = WireHitState::inactive; // disable the hit if it's an outlier
+      state.state_ = WireHitState::inactive; // disable the hit if it's an outlier
     } else if(fabs(doca) > mindoca_ ) {
-      state = doca > 0.0 ? WireHitState::right : WireHitState::left;
+      state.state_ = doca > 0.0 ? WireHitState::right : WireHitState::left;
     } else {
       // too close to the wire: don't try to disambiguate LR sign
-      state = WireHitState::null;
+      state.state_ = WireHitState::null;
     }
-    return state;
   }
 }
 #endif

Examples/ScintHit.hh

@@ -44,7 +44,7 @@ namespace KinKal {
 
   template <class KTRAJ> ScintHit<KTRAJ>::ScintHit(PCA const& pca, double tvar, double wvar) :
     saxis_(pca.sensorTraj()), tvar_(tvar), wvar_(wvar),
-    tpca_(pca.localTraj(),saxis_,pca.precision(),pca.tpData(),pca.dDdP(),pca.dTdP())
+    tpca_(pca.localTraj(),saxis_,pca.precision(),pca.tpData(),pca.dDdP(),pca.dTdP(),pca.dLdP())
   {}
 
   template <class KTRAJ> Residual const& ScintHit<KTRAJ>::refResidual(unsigned ires) const {

Examples/SimpleWireHit.hh

@@ -22,11 +22,11 @@ namespace KinKal {
       using CA = ClosestApproach<KTRAJ,SensorLine>;
       using KTRAJPTR = std::shared_ptr<KTRAJ>;
       using PTRAJ = ParticleTrajectory<KTRAJ>;
-      enum Dimension { dresid=0, tresid=1};  // residual dimensions
+      enum Dimension { dresid=0, tresid=1, lresid=2};  // residual dimensions
 
       SimpleWireHit(BFieldMap const& bfield, PCA const& pca, WireHitState const& whstate, double mindoca,
-          double driftspeed, double tvar, double tot, double totvar, double rcell,int id);
-      unsigned nResid() const override { return 2; } // 2 residuals
+          double driftspeed, double tvar, double tot, double totvar, double longval, double longvar, double rcell,int id);
+      unsigned nResid() const override { return 3; } // 2 residuals
       double time() const override { return ca_.particleToca(); }
       Residual const& refResidual(unsigned ires=dresid) const override;
       void updateReference(PTRAJ const& ptraj) override;
@@ -56,11 +56,12 @@ namespace KinKal {
                   // is the effective signal propagation velocity along the wire, and the time range describes the active wire length
                   // (when multiplied by the propagation velocity).
       CA ca_; // reference time and position of closest approach to the wire; this is generally biased by the hit
-      std::array<Residual,2> rresid_; // residuals WRT most recent reference
+      std::array<Residual,3> rresid_; // residuals WRT most recent reference
       double mindoca_; // effective minimum DOCA used when assigning LR ambiguity, used to define null hit properties
       double dvel_; // constant drift speed
       double tvar_; // constant time variance
       double tot_, totvar_; // TimeOverThreshold and variance
+      double longval_, longvar_; // longitudinal value and variance
       double rcell_; // straw radius
       int id_; // id
       void updateResiduals();
@@ -69,15 +70,20 @@ namespace KinKal {
   //trivial 'updater' that sets the wire hit state to null
   class NullWireHitUpdater {
     public:
-      WireHitState wireHitState() const { return WireHitState(WireHitState::null); }
+      void updateWireHitState(WireHitState &state) const { state.state_ = WireHitState::null; }
+  };
+  //trivial 'updater' that enables longitudinal measurement
+  class LongWireHitUpdater {
+    public:
+      void updateWireHitState(WireHitState &state) const {state.lstate_ = WireHitState::longactive; }
   };
 
   template <class KTRAJ> SimpleWireHit<KTRAJ>::SimpleWireHit(BFieldMap const& bfield, PCA const& pca, WireHitState const& whstate,
-      double mindoca, double driftspeed, double tvar, double tot, double totvar, double rcell, int id) :
+      double mindoca, double driftspeed, double tvar, double tot, double totvar, double longval, double longvar, double rcell, int id) :
     bfield_(bfield),
     whstate_(whstate), wire_(pca.sensorTraj()),
-    ca_(pca.localTraj(),wire_,pca.precision(),pca.tpData(),pca.dDdP(),pca.dTdP()), // must be explicit to get the right sensor traj reference
-    mindoca_(mindoca), dvel_(driftspeed), tvar_(tvar), tot_(tot), totvar_(totvar), rcell_(rcell), id_(id) {
+    ca_(pca.localTraj(),wire_,pca.precision(),pca.tpData(),pca.dDdP(),pca.dTdP(),pca.dLdP()), // must be explicit to get the right sensor traj reference
+    mindoca_(mindoca), dvel_(driftspeed), tvar_(tvar), tot_(tot), totvar_(totvar), longval_(longval), longvar_(longvar), rcell_(rcell), id_(id) {
     }
 
   template <class KTRAJ> void SimpleWireHit<KTRAJ>::updateReference(PTRAJ const& ptraj) {
@@ -93,22 +99,29 @@ namespace KinKal {
     if(first){
       // look for an updater; if found, use it to update the state
       auto nwhu = miconfig.findUpdater<NullWireHitUpdater>();
+      auto twhu = miconfig.findUpdater<LongWireHitUpdater>();
       auto dwhu = miconfig.findUpdater<DOCAWireHitUpdater>();
       if(nwhu != 0 && dwhu != 0)throw std::invalid_argument(">1 SimpleWireHit updater specified");
+      if(twhu != 0){
+        twhu->updateWireHitState(whstate_);
+      }
       if(nwhu != 0){
         mindoca_ = cellRadius();
-        whstate_ = nwhu->wireHitState();
+        nwhu->updateWireHitState(whstate_);
         // set the residuals based on this state
       } else if(dwhu != 0){
         // update minDoca (for null ambiguity error estimate)
         mindoca_ = std::min(dwhu->minDOCA(),cellRadius());
         // compute the unbiased closest approach.  This is brute-force
         // a more clever solution is to linearly correct the residuals for the change in parameters
         auto uca = this->unbiasedClosestApproach();
-        whstate_ = uca.usable() ? dwhu->wireHitState(uca.doca()) : WireHitState(WireHitState::inactive);
+        if (uca.usable())
+          dwhu->updateWireHitState(whstate_,uca.doca());
+        else
+          whstate_ = WireHitState(WireHitState::inactive);
       }
     }
-    rresid_[tresid] = rresid_[dresid] = Residual();
+    rresid_[tresid] = rresid_[dresid] = rresid_[lresid] = Residual();
     if(whstate_.active()){
       rresid_[tresid] = Residual(ca_.deltaT() - tot_, totvar_,0.0,true,ca_.dTdP()); // always constrain to TOT; this stabilizes the fit
       if(whstate_.useDrift()){
@@ -121,13 +134,18 @@ namespace KinKal {
         double nulldvar = dvel_*dvel_*(ca_.deltaT()*ca_.deltaT()+0.8);
         rresid_[dresid] = Residual(ca_.doca(),nulldvar,0.0,true,ca_.dDdP());
       }
+      if(whstate_.useLong()){
+        double calong = (ca_.sensorPoca().Vect() - wire_.middle()).Dot(ca_.sensorDirection());
+        double lr = calong - longval_;
+        rresid_[lresid] = Residual(lr,longvar_,0.0,true,ca_.dLdP());
+      }
     }
     // now update the weight
     this->updateWeight(miconfig);
   }
 
   template <class KTRAJ> Residual const& SimpleWireHit<KTRAJ>::refResidual(unsigned ires) const {
-    if(ires >tresid)throw std::invalid_argument("Invalid residual");
+    if(ires >lresid)throw std::invalid_argument("Invalid residual");
     return rresid_[ires];
   }
 
@@ -155,11 +173,21 @@ namespace KinKal {
         ost << "null";
         break;
     }
+    switch(whstate_.lstate_) {
+      case WireHitState::longinactive:
+        ost << " longinactive";
+        break;
+      case WireHitState::longactive:
+        ost << " longactive";
+        break;
+    }
     if(detail > 0){
       if(rresid_[tresid].active())
         ost << " Active Time Residual " << rresid_[tresid];
       if(rresid_[dresid].active())
         ost << " Active Distance Residual " << rresid_[dresid];
+      if(rresid_[lresid].active())
+        ost << " Active Long Residual " << rresid_[lresid];
       ost << std::endl;
     }
     if(detail > 1) {

Examples/StrawXing.hh

@@ -49,7 +49,7 @@ namespace KinKal {
   template <class KTRAJ> StrawXing<KTRAJ>::StrawXing(PCA const& pca, StrawMaterial const& smat) :
     axis_(pca.sensorTraj()),
     smat_(smat),
-    tpca_(pca.localTraj(),axis_,pca.precision(),pca.tpData(),pca.dDdP(),pca.dTdP()),
+    tpca_(pca.localTraj(),axis_,pca.precision(),pca.tpData(),pca.dDdP(),pca.dTdP(),pca.dLdP()),
     toff_(smat.wireRadius()/pca.particleTraj().speed(pca.particleToca())), // locate the effect to 1 side of the wire to avoid overlap with hits
     varscale_(1.0)
   {}

Examples/WireHitStructs.hh

@@ -13,11 +13,16 @@ namespace KinKal {
   // struct describing wire hit internal state
   struct WireHitState {
     enum State { inactive=-2, left=-1, null=0, right=1};  // state description
+    enum LState { longinactive=0, longactive=1}; // use longitudinal measurement 
     State state_; // left-right ambiguity
+    LState lstate_; // time division
     bool useDrift() const { return state_ == left || state_ == right; }
+    bool useLong() const { return lstate_ == longactive; }
     bool active() const { return state_ != inactive; }
     bool operator == (WireHitState::State state) const { return state_ == state; }
     bool operator != (WireHitState::State state) const { return state_ != state; }
+    bool operator == (WireHitState::LState state) const { return lstate_ == state; }
+    bool operator != (WireHitState::LState state) const { return lstate_ != state; }
     double lrSign() const {
       switch (state_) {
         case left:
@@ -28,7 +33,7 @@ namespace KinKal {
           return 0.0;
       }
     }
-    WireHitState(State state = inactive) : state_(state) {}
+    WireHitState(State state = inactive, LState lstate = longinactive) : state_(state), lstate_(lstate) {}
   };
 }
 #endif

Tests/FitTest.hh

@@ -66,7 +66,7 @@ using namespace KinKal;
 using namespace std;
 // avoid confusion with root
 void print_usage() {
-  printf("Usage: FitTest  --momentum f --simparticle i --fitparticle i--charge i --nhits i --hres f --seed i -ambigdoca f --nevents i --simmat i--fitmat i --ttree i --Bz f --dBx f --dBy f --dBz f--Bgrad f --tolerance f --TFilesuffix c --PrintBad i --PrintDetail i --ScintHit i --invert i --Schedule a --ssmear i --constrainpar i --inefficiency f --extend s --TimeBuffer f --matvarscale i\n");
+  printf("Usage: FitTest  --momentum f --simparticle i --fitparticle i--charge i --nhits i --hres f --seed i -ambigdoca f --nevents i --simmat i--fitmat i --ttree i --Bz f --dBx f --dBy f --dBz f--Bgrad f --tolerance f --TFilesuffix c --PrintBad i --PrintDetail i --ScintHit i --invert i --Schedule a --ssmear i --constrainpar i --inefficiency f --extend s --TimeBuffer f --matvarscale i --uselong i\n");
 }
 
 // utility function to compute transverse distance between 2 similar trajectories.  Also
@@ -147,7 +147,7 @@ int testState(KinKal::Track<KTRAJ> const& kktrk,DVEC sigmas) {
   return retval;
 }
 
-int makeConfig(string const& cfile, KinKal::Config& config,bool mvarscale=true) {
+int makeConfig(string const& cfile, KinKal::Config& config,bool mvarscale=true, bool uselong=false) {
   string fullfile;
   if(strncmp(cfile.c_str(),"/",1) == 0) {
     fullfile = string(cfile);
@@ -181,6 +181,8 @@ int makeConfig(string const& cfile, KinKal::Config& config,bool mvarscale=true)
         ss >> temp >> utype;
         MetaIterConfig miconfig(temp);
         miconfig.addUpdater(StrawXingConfig(0.3,5.0,10.0,mvarscale)); // hardcoded values, should come from outside, FIXME
+        if(uselong)
+          miconfig.addUpdater(std::any(LongWireHitUpdater()));
         if(utype == 0 ){
           cout << "NullWireHitUpdater for iteration " << nmiter << endl;
           miconfig.addUpdater(std::any(NullWireHitUpdater()));
@@ -255,6 +257,7 @@ int FitTest(int argc, char *argv[],KinKal::DVEC const& sigmas) {
   bool fitmat(true);
   bool extend(false);
   bool mvarscale(true);
+  bool uselong(false);
   string exfile;
   BFieldMap *BF(0);
   double Bgrad(0.0), dBx(0.0), dBy(0.0), dBz(0.0), Bz(1.0);
@@ -302,6 +305,7 @@ int FitTest(int argc, char *argv[],KinKal::DVEC const& sigmas) {
     {"extend",     required_argument, 0, 'X'  },
     {"TimeBuffer",     required_argument, 0, 'W'  },
     {"MatVarScale",     required_argument, 0, 'v'  },
+    {"uselong", required_argument, 0, 'l' },
     {NULL, 0,0,0}
   };
 
@@ -368,6 +372,8 @@ int FitTest(int argc, char *argv[],KinKal::DVEC const& sigmas) {
       case 'X' : exfile = optarg;
                  extend = true;
                  break;
+      case 'l' : uselong = atoi(optarg);
+                 break;
       default: print_usage();
                exit(EXIT_FAILURE);
     }
@@ -392,12 +398,12 @@ int FitTest(int argc, char *argv[],KinKal::DVEC const& sigmas) {
   toy.setTolerance(tol/10.0); // finer precision on sim
   // setup fit configuration
   Config config;
-  makeConfig(sfile,config,mvarscale);
+  makeConfig(sfile,config,mvarscale,uselong);
   cout << "Main fit " << config << endl;
   // read the schedule from the file
   Config exconfig;
   if(extend){
-    makeConfig(exfile,exconfig,mvarscale);
+    makeConfig(exfile,exconfig,mvarscale,uselong);
     cout << "Extension " << exconfig << endl;
   }
   // generate hits

Tests/ToyMC.hh

@@ -45,7 +45,7 @@ namespace KKTest {
         mom_(mom), icharge_(icharge),
         tr_(iseed), nhits_(nhits), simmat_(simmat), scinthit_(scinthit), ambigdoca_(ambigdoca), simmass_(simmass),
         sprop_(0.8*CLHEP::c_light), sdrift_(0.065),
-        zrange_(zrange), rstraw_(2.5), rwire_(0.025), wthick_(0.015), wlen_(1000.0), sigt_(3.0), sigtot_(7.0), ineff_(0.05),
+        zrange_(zrange), rstraw_(2.5), rwire_(0.025), wthick_(0.015), wlen_(1000.0), sigt_(3.0), sigtot_(7.0), siglong_(40.0), ineff_(0.05),
         scitsig_(0.5), shPosSig_(10.0), shmax_(40.0), cz_(zrange_+50), clen_(200.0), cprop_(0.8*CLHEP::c_light),
         osig_(10.0), ctmin_(0.5), ctmax_(0.8), tol_(1e-5), tprec_(1e-8), t0off_(700.0),
         smat_(matdb_,rstraw_, wthick_, 3*wthick_, rwire_),
@@ -93,6 +93,7 @@ namespace KKTest {
       double rwire_, wthick_, wlen_; // wire radius, thickness, length
       double sigt_; // drift time resolution in ns
       double sigtot_; // TOT drift time resolution (ns)
+      double siglong_; // Longitudinal pos resolution (mm)
       double ineff_; // hit inefficiency
                      // time hit parameters
       double scitsig_, shPosSig_, shmax_, cz_, clen_, cprop_;
@@ -155,7 +156,9 @@ namespace KKTest {
         double mindoca = std::min(ambigdoca_,rstraw_);
         // true TOT is the drift time
         double tot = tr_.Gaus(tp.deltaT(),sigtot_);
-        thits.push_back(std::make_shared<WIREHIT>(bfield_, tp, whstate, mindoca, sdrift_, sigt_*sigt_, tot, sigtot_*sigtot_,  rstraw_, ihit));
+        double longpos = (tp.sensorPoca().Vect()-tline.middle()).Dot(tp.sensorDirection());
+        double longval = tr_.Gaus(longpos,siglong_);
+        thits.push_back(std::make_shared<WIREHIT>(bfield_, tp, whstate, mindoca, sdrift_, sigt_*sigt_, tot, sigtot_*sigtot_, longval, siglong_*siglong_, rstraw_, ihit));
       }
       // compute material effects and change trajectory accordingly
       auto xing = std::make_shared<STRAWXING>(tp,smat_);

Trajectory/CentralHelix.cc

@@ -93,9 +93,11 @@ namespace KinKal {
     setBNom(trot,bnom);
   }
 
-  CentralHelix::CentralHelix(Parameters const &pdata, double mass, int charge, double bnom, TimeRange const& range) : trange_(range),  pars_(pdata), mass_(mass), bnom_(VEC3(0.0,0.0,bnom)){
+  CentralHelix::CentralHelix(Parameters const &pdata, double mass, int charge, double bnom, TimeRange const& range) : CentralHelix(pdata, mass, charge, VEC3(0.0,0.0,bnom), range) {}
+  CentralHelix::CentralHelix(Parameters const &pdata, double mass, int charge, VEC3 const& bnom, TimeRange const& range) : trange_(range),  pars_(pdata), mass_(mass), bnom_(bnom){
     //FIXME for now just ignore sign
     abscharge_ = abs(charge);
+    setTransforms();
   }
 
   CentralHelix::CentralHelix(Parameters const &pdata, CentralHelix const& other) : CentralHelix(other) {
@@ -442,7 +444,7 @@ namespace KinKal {
     // work in local coordinate system to avoid additional matrix mulitplications
     auto xvec = localPosition(time);
     auto mvec = localMomentum(time);
-    VEC3 BxdB =VEC3(0.0,0.0,1.0).Cross(dB)/bnomR();
+    VEC3 BxdB =VEC3(0.0,0.0,1.0).Cross(dBloc)/bnomR();
     VEC3 dx = xvec.Cross(BxdB);
     VEC3 dm = mvec.Cross(BxdB);
     // convert these to a full state vector change

Trajectory/CentralHelix.hh

@@ -46,6 +46,7 @@ namespace KinKal {
       CentralHelix(VEC4 const& pos, MOM4 const& mom, int charge, VEC3 const& bnom, TimeRange const& range=TimeRange());
       CentralHelix(VEC4 const& pos, MOM4 const& mom, int charge, double bnom, TimeRange const& range=TimeRange());
       // construct from explicit parametric and kinematic info
+      CentralHelix(Parameters const &pdata, double mass, int charge, VEC3 const& bnom, TimeRange const& range);
       CentralHelix(Parameters const &pdata, double mass, int abscharge, double bnom, TimeRange const& range);
       // copy payload and adjust for a different BFieldMap and range
       CentralHelix(CentralHelix const& other, VEC3 const& bnom, double trot);