Alignment updates

offline/packages/TrackerMillepedeAlignment/HelicalFitter.h

@@ -76,8 +76,15 @@ class HelicalFitter : public SubsysReco, public PHParameterInterface
   void set_vertex_param_fixed(unsigned int param){ fixed_vertex_params.insert(param);}
   void set_straight_line_fit(bool flag) {straight_line_fit = flag; }
   void set_eta_cut(double eta_cut) {m_eta_cut = eta_cut;}
+  void set_pt_cut(double pt_cut) {m_pt_min = pt_cut;}
+  void set_acts_mode(bool acts_mode) {m_acts_mode = acts_mode;}
+  void set_fixed_vtx(bool fixed_vtx) {m_fixed_vtx = fixed_vtx;}
+  void set_fixed_vtx_x(float fixed_vtx_x) {m_fixed_vtx_x = fixed_vtx_x;}
+  void set_fixed_vtx_y(float fixed_vtx_y) {m_fixed_vtx_y = fixed_vtx_y;}
+  //Fixes MVTX half in order to calculate projected residuals from east to west half or vice versa.
   //-1 is regular operation, 0 is east fixed, 1 is west fixed
-  void set_do_mvtx_half(int half) {do_mvtx_half = half; }
+  void set_do_mvtx_half(int half) {m_do_mvtx_half = half; }
+  void set_is_cosmics() {is_cosmics=true;}
   void set_fitted_subsystems(bool si, bool tpc, bool full)
   {
     fitsilicon = si;
@@ -122,7 +129,7 @@ class HelicalFitter : public SubsysReco, public PHParameterInterface
   Acts::Vector3 getPCALinePoint(const Acts::Vector3& global, const Acts::Vector3& tangent, const Acts::Vector3& posref);
   Acts::Vector3 get_line_plane_intersection(const Acts::Vector3& PCA, const Acts::Vector3& tangent,
                                             const Acts::Vector3& sensor_center, const Acts::Vector3& sensor_normal);
-  std::pair<Acts::Vector3, Acts::Vector3> get_helix_tangent(const std::vector<float>& fitpars, Acts::Vector3 global);
+  std::pair<Acts::Vector3, Acts::Vector3> get_helix_tangent(const std::vector<float>& fitpars, Acts::Vector3 global, bool is_cosmics);
   Acts::Vector3 get_helix_surface_intersection(const Surface& surf, std::vector<float>& fitpars, Acts::Vector3 global);
   Acts::Vector3 get_helix_surface_intersection(const Surface& surf, std::vector<float>& fitpars, Acts::Vector3 global, Acts::Vector3& pca, Acts::Vector3& tangent);
 
@@ -200,10 +207,17 @@ class HelicalFitter : public SubsysReco, public PHParameterInterface
   bool fitsilicon{true};
   bool fittpc{false};
   bool fitfulltrack{false};
+  bool m_acts_mode{false};
+  bool m_fixed_vtx{false};
+
+  float m_fixed_vtx_x{0.0};
+  float m_fixed_vtx_y{0.0};
+
 
   float dca_cut{0.19};  // cm
 
   float m_eta_cut{99999.};
+  float m_pt_min{0.5};
 
   SvtxVertexMap* m_vertexmap{nullptr};
   SvtxTrackMap* m_trackmap{nullptr};
@@ -224,7 +238,8 @@ class HelicalFitter : public SubsysReco, public PHParameterInterface
   bool use_event_vertex{false};
   bool use_intt_zfit{false};
   bool straight_line_fit = false;
-  int do_mvtx_half = -1;
+  int m_do_mvtx_half = -1;
+  bool is_cosmics {false};
 
   int event{0};
 

offline/packages/TrackerMillepedeAlignment/MakeMilleFiles.cc

@@ -5,6 +5,7 @@
 /// Tracking includes
 
 #include <trackbase/MvtxDefs.h>
+#include <trackbase/InttDefs.h>
 #include <trackbase/TpcDefs.h>  // for side
 #include <trackbase/TrackFitUtils.h>
 #include <trackbase/TrkrCluster.h>
@@ -96,9 +97,11 @@ int MakeMilleFiles::InitRun(PHCompositeNode* topNode)
     m_file = TFile::Open(m_tfile_name.c_str(), "RECREATE");
     m_ntuple = new TNtuple (
       "ntp", "ntp",
-      "dXdR:dXdX0:dXdY0:dXdZs:dXdZ0:"
+      "layer:stave:chip:residualX:residualY:clusphi:xglob:yglob:zglob:"
+      "errX:errY:l0:l1:phi:theta:qoverp:time:"
+      "dXdR:dXdZ0:dXdphi:dXdtheta:dXdqoverp:dXdt:"
       "dXdalpha:dXdbeta:dXdgamma:dXdx:dXdy:dXdz:"
-      "dYdR:dYdX0:dYdY0:dYdZs:dYdZ0:"
+      "dYdR:dYdZ0:dYdphi:dYdtheta:dYdqoverp:dYdt:"
       "dYdalpha:dYdbeta:dYdgamma:dYdx:dYdy:dYdz"
     );
     m_ntuple->SetDirectory(m_file);
@@ -165,6 +168,14 @@ int MakeMilleFiles::process_event(PHCompositeNode* /*topNode*/)
 
     //! Make any desired track cuts here
     //! Maybe set a lower pT limit - low pT tracks are not very sensitive to alignment
+    if(track->get_pt() < m_minPt)
+    {
+      if (Verbosity() > 0)
+      {
+        std::cout << "Skipping track with pT " << track->get_pt() << " < " << m_minPt << std::endl;
+      }
+      continue;
+    }
     addTrackToMilleFile(statevec);
 
     //! Only take tracks that have 2 mm within event vertex
@@ -478,16 +489,36 @@ void MakeMilleFiles::addTrackToMilleFile(SvtxAlignmentStateMap::StateVec& statev
     TrkrCluster* cluster = _cluster_map->findCluster(ckey);
     const unsigned int layer = TrkrDefs::getLayer(ckey);
     const unsigned int trkrid = TrkrDefs::getTrkrId(ckey);
-    const SvtxAlignmentState::ResidualVector residual = state->get_residual();
+     uint8_t stave = -1;
+      uint8_t chip = -1;
+      if (trkrid == TrkrDefs::mvtxId)
+      {
+        // need stave to get clamshell
+         stave = MvtxDefs::getStaveId(ckey);
+         chip = MvtxDefs::getChipId(ckey);
+      }
+      else if(trkrid == TrkrDefs::inttId)
+      {
+        stave = InttDefs::getLadderZId(ckey);
+        chip = InttDefs::getLadderPhiId(ckey);
+      }
+    if(m_ignore_tpc && trkrid == TrkrDefs::tpcId)
+      continue;
+    const SvtxAlignmentState::ResidualVector residual = state->get_residual();// / Acts::UnitConstants::cm;
+    //auto acts_pars = state->parameters();
+    //std::cout<<"acts_par(0): "<<acts_pars(0)<<std::endl;
     const Acts::Vector3 global = _tGeometry->getGlobalPosition(ckey, cluster);
 
     // need standard deviation of measurements
     SvtxAlignmentState::ResidualVector clus_sigma = SvtxAlignmentState::ResidualVector::Zero();
 
-    double clusRadius = sqrt(global[0] * global[0] + global[1] * global[1]);
-    auto para_errors = _ClusErrPara.get_clusterv5_modified_error(cluster, clusRadius, ckey);
-    double phierror = sqrt(para_errors.first);
-    double zerror = sqrt(para_errors.second);
+    //double clusRadius = sqrt(global[0] * global[0] + global[1] * global[1]);
+    double clusphi = atan2(global[1] , global[0]);
+    //auto para_errors = _ClusErrPara.get_clusterv5_modified_error(cluster, clusRadius, ckey);
+    //double phierror = sqrt(para_errors.first);
+    //double zerror = sqrt(para_errors.second);
+    double phierror = cluster->getRPhiError();
+    double zerror = cluster->getZError();
     clus_sigma(1) = zerror * Acts::UnitConstants::cm;
     clus_sigma(0) = phierror * Acts::UnitConstants::cm;
 
@@ -502,11 +533,11 @@ void MakeMilleFiles::addTrackToMilleFile(SvtxAlignmentStateMap::StateVec& statev
     int glbl_label[SvtxAlignmentState::NGL];
     if (layer < 3)
     {
-      AlignmentDefs::getMvtxGlobalLabels(surf, glbl_label, mvtx_group);
+      AlignmentDefs::getMvtxGlobalLabels(surf, ckey, glbl_label, mvtx_group);
     }
     else if (layer > 2 && layer < 7)
     {
-      AlignmentDefs::getInttGlobalLabels(surf, glbl_label, intt_group);
+      AlignmentDefs::getInttGlobalLabels(surf, ckey, glbl_label, intt_group);
     }
     else if (layer < 55)
     {
@@ -524,15 +555,20 @@ void MakeMilleFiles::addTrackToMilleFile(SvtxAlignmentStateMap::StateVec& statev
 
     float glbl_derivative[SvtxAlignmentState::NRES][SvtxAlignmentState::NGL]{};
     float lcl_derivative[SvtxAlignmentState::NRES][SvtxAlignmentState::NLOC]{};
-
+    float lcl_derivative_psuedo[SvtxAlignmentState::NLOC][SvtxAlignmentState::NLOC]{};
+    float lcl_meas_psuedo[SvtxAlignmentState::NLOC]{};
+    float glbl_derivative_dummy_empty[SvtxAlignmentState::NGL]{};
+    int glbl_label_dummy_empty[SvtxAlignmentState::NGL]{};
+    SvtxAlignmentState::ActsTrackParamsVector lcl_trackpars = state->get_acts_track_params();
+
     /// For N residual local coordinates x, z
     for (int i = 0; i < SvtxAlignmentState::NRES; ++i)
     {
+
       // Add the measurement separately for each coordinate direction to Mille
       for (int j = 0; j < SvtxAlignmentState::NGL; ++j)
       {
         glbl_derivative[i][j] = state->get_global_derivative_matrix()(i, j);
-
         if (is_layer_fixed(layer) ||
             is_layer_param_fixed(layer, j, fixed_layer_gparams))
         {
@@ -542,7 +578,7 @@ void MakeMilleFiles::addTrackToMilleFile(SvtxAlignmentStateMap::StateVec& statev
         if (trkrid == TrkrDefs::mvtxId)
         {
           // need stave to get clamshell
-          auto stave = MvtxDefs::getStaveId(ckey);
+           auto stave = MvtxDefs::getStaveId(ckey);
           auto clamshell = AlignmentDefs::getMvtxClamshell(layer, stave);
           if (is_layer_param_fixed(layer, j, fixed_layer_gparams) ||
               is_mvtx_layer_fixed(layer, clamshell))
@@ -552,6 +588,7 @@ void MakeMilleFiles::addTrackToMilleFile(SvtxAlignmentStateMap::StateVec& statev
         }
         else if (trkrid == TrkrDefs::inttId)
         {
+
           if (is_layer_param_fixed(layer, j, fixed_layer_gparams) ||
               is_layer_fixed(layer))
           {
@@ -619,10 +656,24 @@ void MakeMilleFiles::addTrackToMilleFile(SvtxAlignmentStateMap::StateVec& statev
       }
     }
 
+    // For 6 local track parameters, the associated psuedo mille measurements that complete the local derivative
+    for (int i = 0; i < SvtxAlignmentState::NLOC; ++i)
+    {
+      for (int j = 0; j < SvtxAlignmentState::NLOC; ++j)
+      {
+        lcl_derivative_psuedo[i][j] = state->get_local_derivative_psuedo_matrix()(i, j);
+      }
+      lcl_meas_psuedo[i] = state->get_local_psuedo_measurement_err()(i);
+
+      _mille->mille(SvtxAlignmentState::NLOC, lcl_derivative_psuedo[i], SvtxAlignmentState::NGL, glbl_derivative_dummy_empty, glbl_label_dummy_empty, 0.0, lcl_meas_psuedo[i]);
+    }
+
     float ntp_data[] = {
-      lcl_derivative[0][0], lcl_derivative[0][1], lcl_derivative[0][2], lcl_derivative[0][3], lcl_derivative[0][4],
+      (float) layer, (float) stave, (float) chip, (float) residual(0), (float) residual(1), (float) clusphi, (float) global[0], (float) global[1], (float) global[2],
+      (float) clus_sigma(0), (float) clus_sigma(1), (float) lcl_trackpars(0), (float) lcl_trackpars(1), (float) lcl_trackpars(2), (float) lcl_trackpars(3), (float) lcl_trackpars(4), (float) lcl_trackpars(5),
+      lcl_derivative[0][0], lcl_derivative[0][1], lcl_derivative[0][2], lcl_derivative[0][3], lcl_derivative[0][4], lcl_derivative[0][5],
       glbl_derivative[0][0], glbl_derivative[0][1], glbl_derivative[0][2], glbl_derivative[0][3], glbl_derivative[0][4], glbl_derivative[0][5],
-      lcl_derivative[1][0], lcl_derivative[1][1], lcl_derivative[1][2], lcl_derivative[1][3], lcl_derivative[1][4],
+      lcl_derivative[1][0], lcl_derivative[1][1], lcl_derivative[1][2], lcl_derivative[1][3], lcl_derivative[1][4], lcl_derivative[1][5],
       glbl_derivative[1][0], glbl_derivative[1][1], glbl_derivative[1][2], glbl_derivative[1][3], glbl_derivative[1][4], glbl_derivative[1][5],
     };
 
@@ -707,3 +758,5 @@ bool MakeMilleFiles::is_tpc_sector_fixed(unsigned int layer, unsigned int sector
 
   return ret;
 }
+
+

offline/packages/TrackerMillepedeAlignment/MakeMilleFiles.h

@@ -69,6 +69,10 @@ class MakeMilleFiles : public SubsysReco
   void set_layer_gparam_fixed(unsigned int layer, unsigned int param);
   void set_layer_lparam_fixed(unsigned int layer, unsigned int param);
 
+  void set_pt_cut(float pt) { m_minPt = pt; }
+
+  void set_ignore_tpc() {m_ignore_tpc = true;}
+
   void set_layers_fixed(unsigned int minlayer, unsigned int maxlayer);
   void set_error_inflation_factor(unsigned int layer, float factor)
   {
@@ -118,6 +122,7 @@ class MakeMilleFiles : public SubsysReco
 
   bool m_useEventVertex = false;
   bool _binary = true;
+  float m_minPt = 0.0;
 
   Acts::Vector2 m_vtxSigma = {0.1, 0.1};
 
@@ -133,6 +138,8 @@ class MakeMilleFiles : public SubsysReco
   std::set<std::pair<unsigned int, unsigned int>> fixed_mvtx_layers;
   std::set<std::pair<unsigned int, unsigned int>> fixed_layer_gparams, fixed_layer_lparams;
 
+  bool m_ignore_tpc = false;
+
   std::string m_constraintFileName = "mp2con.txt";
   std::ofstream m_constraintFile;
 

offline/packages/TrackingDiagnostics/TrackResiduals.cc

@@ -538,6 +538,12 @@ void TrackResiduals::fillVertexTree(PHCompositeNode* topNode)
         {
           continue;
         }
+        m_pcax_vtx_trk.push_back(track->get_x());
+        m_pcay_vtx_trk.push_back(track->get_y());
+        m_pcaz_vtx_trk.push_back(track->get_z());
+        m_px_vtx_trk.push_back(track->get_px());
+        m_py_vtx_trk.push_back(track->get_py());
+        m_pz_vtx_trk.push_back(track->get_pz());
         for (const auto& ckey : get_cluster_keys(track))
         {
           TrkrCluster* cluster = clustermap->findCluster(ckey);
@@ -600,7 +606,7 @@ void TrackResiduals::circleFitClusters(
 
   auto xyparams = TrackFitUtils::line_fit(xypoints);
   auto yzLineParams = TrackFitUtils::line_fit(yzpoints);
-  auto fitpars = TrackFitUtils::fitClusters(global_vec, keys, false);
+  auto fitpars = TrackFitUtils::fitClusters(global_vec, keys, false,false,false,true);
   // auto fitpars = TrackFitUtils::fitClusters(global_vec, keys, !m_linefitTPCOnly);
   m_xyint = std::get<1>(xyparams);
   m_xyslope = std::get<0>(xyparams);
@@ -1783,6 +1789,12 @@ void TrackResiduals::createBranches()
   m_vertextree->Branch("gz", &m_clusgz);
   m_vertextree->Branch("gr", &m_clusgr);
   m_vertextree->Branch("mbdcharge", &m_totalmbd, "m_totalmbd/F");
+  m_vertextree->Branch("pcax_vtx_trk", &m_pcax_vtx_trk);
+  m_vertextree->Branch("pcay_vtx_trk", &m_pcay_vtx_trk);
+  m_vertextree->Branch("pcaz_vtx_trk", &m_pcaz_vtx_trk);
+  m_vertextree->Branch("px_vtx_trk", &m_px_vtx_trk);
+  m_vertextree->Branch("py_vtx_trk", &m_py_vtx_trk);
+  m_vertextree->Branch("pz_vtx_trk", &m_pz_vtx_trk);
 
   m_hittree = new TTree("hittree", "A tree with all hits");
   m_hittree->Branch("run", &m_runnumber, "m_runnumber/I");

offline/packages/TrackingDiagnostics/TrackResiduals.h

@@ -243,6 +243,12 @@ class TrackResiduals : public SubsysReco
 
   int m_ntracks = std::numeric_limits<int>::quiet_NaN();
   int m_nvertices = std::numeric_limits<int>::quiet_NaN();
+  std::vector<float> m_pcax_vtx_trk;
+  std::vector<float> m_pcay_vtx_trk;
+  std::vector<float> m_pcaz_vtx_trk;
+  std::vector<float> m_px_vtx_trk;
+  std::vector<float> m_py_vtx_trk;
+  std::vector<float> m_pz_vtx_trk;
 
   //! cluster tree info
   float m_sclusgr = std::numeric_limits<float>::quiet_NaN();

offline/packages/TrackingDiagnostics/TrackSeedTrackMapConverter.cc

@@ -324,7 +324,7 @@ int TrackSeedTrackMapConverter::process_event(PHCompositeNode* /*unused*/)
       float Y0 = trackSeed->get_Y0();
       float Z0 = trackSeed->get_Z0();
       float slope = trackSeed->get_slope();
-
+      
       std::vector<float> fitpars = {R, X0, Y0, slope, Z0};
 
       svtxtrack->set_crossing(trackSeed->get_crossing());
@@ -447,7 +447,8 @@ int TrackSeedTrackMapConverter::process_event(PHCompositeNode* /*unused*/)
       std::cout << "Inserting svtxtrack into map " << std::endl;
       svtxtrack->identify();
     }
-
+    //if(m_cosmics)
+
     m_trackMap->insert(svtxtrack.get());
   }