diff --git a/offline/packages/TrackerMillepedeAlignment/HelicalFitter.cc b/offline/packages/TrackerMillepedeAlignment/HelicalFitter.cc index 8331a8b491..76d82c2130 100644 --- a/offline/packages/TrackerMillepedeAlignment/HelicalFitter.cc +++ b/offline/packages/TrackerMillepedeAlignment/HelicalFitter.cc @@ -90,7 +90,7 @@ HelicalFitter::HelicalFitter(const std::string& name) vertexPosition(1) = 0; vtx_sigma(0) = 0.005; - vtx_sigma(1) = 0.005; + vtx_sigma(1) = 0.01; } //____________________________________________________________________________.. @@ -131,20 +131,20 @@ int HelicalFitter::InitRun(PHCompositeNode* topNode) fout = new TFile(ntuple_outfilename.c_str(), "recreate"); if (straight_line_fit) { - ntp = new TNtuple("ntp", "HF ntuple", "event:trkid:layer:nsilicon:crosshalfmvtx:ntpc:nclus:trkrid:sector:side:subsurf:phi:glbl0:glbl1:glbl2:glbl3:glbl4:glbl5:sensx:sensy:sensz:normx:normy:normz:sensxideal:sensyideal:senszideal:normxideal:normyideal:normzideal:xglobideal:yglobideal:zglobideal:XYs:Y0:Zs:Z0:xglob:yglob:zglob:xfit:yfit:zfit:xfitMvtxHalf:yfitMvtxHalf:zfitMvtxHalf:pcax:pcay:pcaz:tangx:tangy:tangz:X:Y:fitX:fitY:fitXMvtxHalf:fitYMvtxHalf:dXdXYs:dXdY0:dXdZs:dXdZ0:dXdalpha:dXdbeta:dXdgamma:dXdx:dXdy:dXdz:dYdXYs:dYdY0:dYdZs:dYdZ0:dYdalpha:dYdbeta:dYdgamma:dYdx:dYdy:dYdz"); + ntp = new TNtuple("ntp", "HF ntuple", "event:trkid:layer:nsilicon:quality:crosshalfmvtx:ntpc:nclus:trkrid:sector:side:subsurf:phi:glbl0:glbl1:glbl2:glbl3:glbl4:glbl5:sensx:sensy:sensz:normx:normy:normz:sensxideal:sensyideal:senszideal:normxideal:normyideal:normzideal:xglobideal:yglobideal:zglobideal:XYs:Y0:Zs:Z0:xglob:yglob:zglob:xfit:yfit:zfit:xfitMvtxHalf:yfitMvtxHalf:zfitMvtxHalf:pcax:pcay:pcaz:tangx:tangy:tangz:X:Y:fitX:fitY:fitXMvtxHalf:fitYMvtxHalf:dXdXYs:dXdY0:dXdZs:dXdZ0:dXdalpha:dXdbeta:dXdgamma:dXdx:dXdy:dXdz:dYdXYs:dYdY0:dYdZs:dYdZ0:dYdalpha:dYdbeta:dYdgamma:dYdx:dYdy:dYdz"); } else { - ntp = new TNtuple("ntp", "HF ntuple", "event:trkid:layer:nsilicon:ntpc:nclus:trkrid:sector:side:subsurf:phi:glbl0:glbl1:glbl2:glbl3:glbl4:glbl5:sensx:sensy:sensz:normx:normy:normz:sensxideal:sensyideal:senszideal:normxideal:normyideal:normzideal:xglobideal:yglobideal:zglobideal:R:X0:Y0:Zs:Z0:xglob:yglob:zglob:xfit:yfit:zfit:pcax:pcay:pcaz:tangx:tangy:tangz:X:Y:fitX:fitY:dXdR:dXdX0:dXdY0:dXdZs:dXdZ0:dXdalpha:dXdbeta:dXdgamma:dXdx:dXdy:dXdz:dYdR:dYdX0:dYdY0:dYdZs:dYdZ0:dYdalpha:dYdbeta:dYdgamma:dYdx:dYdy:dYdz"); + ntp = new TNtuple("ntp", "HF ntuple", "event:trkid:layer:nsilicon:crosshalfmvtx:ntpc:nclus:trkrid:quality:charge:crossing:sector:side:subsurf:phi:glbl0:glbl1:glbl2:glbl3:glbl4:glbl5:sensx:sensy:sensz:normx:normy:normz:sensxideal:sensyideal:senszideal:normxideal:normyideal:normzideal:xglobideal:yglobideal:zglobideal:R:X0:Y0:Zs:Z0:xglob:yglob:zglob:xfit:yfit:zfit:pcax:pcay:pcaz:tangx:tangy:tangz:X:Y:errX:errY:fitX:fitY:dXdR:dXdX0:dXdY0:dXdZs:dXdZ0:dXdalpha:dXdbeta:dXdgamma:dXdx:dXdy:dXdz:dYdR:dYdX0:dYdY0:dYdZs:dYdZ0:dYdalpha:dYdbeta:dYdgamma:dYdx:dYdy:dYdz"); } if (straight_line_fit) { - track_ntp = new TNtuple("track_ntp", "HF track ntuple", "track_id:residual_x:residual_y:residualxsigma:residualysigma:dXdXYs:dXdY0:dXdZs:dXdZ0:dXdx:dXdy:dXdz:dYdXYs:dYdY0:dYdZs:dYdZ0:dYdx:dYdy:dYdz:track_xvtx:track_yvtx:track_zvtx:event_xvtx:event_yvtx:event_zvtx:track_phi:perigee_phi:track_eta"); + track_ntp = new TNtuple("track_ntp", "HF track ntuple", "track_id:quality:residual_x:residual_y:residualxsigma:residualysigma:dXdXYs:dXdY0:dXdZs:dXdZ0:dXdx:dXdy:dXdz:dYdXYs:dYdY0:dYdZs:dYdZ0:dYdx:dYdy:dYdz:track_xvtx:track_yvtx:track_zvtx:event_xvtx:event_yvtx:event_zvtx:track_phi:perigee_phi:track_eta"); } else { - track_ntp = new TNtuple("track_ntp", "HF track ntuple", "track_id:residual_x:residual_y:residualxsigma:residualysigma:dXdR:dXdX0:dXdY0:dXdZs:dXdZ0:dXdx:dXdy:dXdz:dYdR:dYdX0:dYdY0:dYdZs:dYdZ0:dYdx:dYdy:dYdz:track_xvtx:track_yvtx:track_zvtx:event_xvtx:event_yvtx:event_zvtx:track_phi:perigee_phi"); + track_ntp = new TNtuple("track_ntp", "HF track ntuple", "track_id:quality:residual_x:residual_y:residualxsigma:residualysigma:dXdR:dXdX0:dXdY0:dXdZs:dXdZ0:dXdx:dXdy:dXdz:dYdR:dYdX0:dYdY0:dYdZs:dYdZ0:dYdx:dYdy:dYdz:track_xvtx:track_yvtx:track_zvtx:event_xvtx:event_yvtx:event_zvtx:quality:track_phi:perigee_phi:track_eta:track_p:track_pt"); } } @@ -223,14 +223,15 @@ int HelicalFitter::process_event(PHCompositeNode* /*unused*/) unsigned int nsilicon = 0; unsigned int ntpc = 0; unsigned int nclus = 0; + int trkid_correction = 0; bool h2h_flag = false; bool mvtx_east_only = false; bool mvtx_west_only = false; - if (do_mvtx_half == 0) + if (m_do_mvtx_half == 0) { mvtx_east_only = true; } - if (do_mvtx_half == 1) + if (m_do_mvtx_half == 1) { mvtx_west_only = true; } @@ -241,6 +242,8 @@ int HelicalFitter::process_event(PHCompositeNode* /*unused*/) std::vector cumulative_vertex; std::vector cumulative_someseed; std::vector cumulative_newTrack; + std::vector cumulative_trackid; + std::vector cumulative_event_vertex; if (fittpc && _track_map_tpc != nullptr) { @@ -263,6 +266,7 @@ int HelicalFitter::process_event(PHCompositeNode* /*unused*/) } if (!tracklet) { + trkid_correction++; continue; } @@ -273,6 +277,7 @@ int HelicalFitter::process_event(PHCompositeNode* /*unused*/) getTrackletClusterList(tracklet, cluskey_vec); if (cluskey_vec.size() < 3) { + trkid_correction++; continue; } int nintt = 0; @@ -283,6 +288,8 @@ int HelicalFitter::process_event(PHCompositeNode* /*unused*/) nintt++; } } + if(nintt<2) + continue; // store cluster global positions in a vector global_vec and cluskey_vec @@ -314,7 +321,7 @@ int HelicalFitter::process_event(PHCompositeNode* /*unused*/) } else { - if (fitsilicon && nintt < 2) + if (fitsilicon && global_vec.size() < 4) { continue; // discard incomplete seeds } @@ -323,8 +330,12 @@ int HelicalFitter::process_event(PHCompositeNode* /*unused*/) continue; } - fitpars = TrackFitUtils::fitClusters(global_vec, cluskey_vec); // do helical fit - + fitpars = TrackFitUtils::fitClusters(global_vec, cluskey_vec,false,false,false,is_cosmics); // do helical fit + fitpars_mvtx_half = TrackFitUtils::fitClusters(global_vec, cluskey_vec, use_intt_zfit, mvtx_east_only, mvtx_west_only, is_cosmics); + if (fitpars_mvtx_half.size() < 3) + { + fitpars_mvtx_half = fitpars; + } if (fitpars.size() == 0) { continue; // discard this track, not enough clusters to fit @@ -406,7 +417,8 @@ int HelicalFitter::process_event(PHCompositeNode* /*unused*/) } else { - fitpars = TrackFitUtils::fitClusters(global_vec, cluskey_vec, use_intt_zfit); // do helical fit + fitpars = TrackFitUtils::fitClusters(global_vec, cluskey_vec, use_intt_zfit,false,false,is_cosmics); // do helical fit + fitpars_mvtx_half = TrackFitUtils::fitClusters(global_vec, cluskey_vec, use_intt_zfit, mvtx_east_only, mvtx_west_only, is_cosmics); if (fitpars.size() == 0) { @@ -519,10 +531,60 @@ int HelicalFitter::process_event(PHCompositeNode* /*unused*/) } continue; } - if (std::abs(newTrack.get_eta()) > m_eta_cut) + //std::cout<<"tracklet eta2: "<get_eta()< m_eta_cut || newTrack.get_pt() < m_pt_min) { continue; } + + Acts::Vector3 event_vtx(0.0,0.0,0.0); + bool passed_vtx_flag = false; + unsigned int abs_cross = newTrack.get_crossing(); + if (newTrack.get_crossing() < 0) + abs_cross = INT_MAX; + if(use_event_vertex) + { + for (const auto& [vtxkey, vertex] : *m_vertexmap) + { + for (auto trackiter = vertex->begin_tracks(); trackiter != vertex->end_tracks(); ++trackiter) + { + if (!m_acts_mode && (trackid-trkid_correction) == (*trackiter) && vertex->size_tracks()>3)//&&(newTrack.get_crossing())==short(vertex->get_beam_crossing()) )//&& vtxtrack->get_crossing()== newTrack.get_crossing()) + { + event_vtx(0) = vertex->get_x(); + event_vtx(1) = vertex->get_y(); + event_vtx(2) = vertex->get_z(); + passed_vtx_flag = true; + if (Verbosity() > 0) + { + std::cout << "vtx crossing: "<get_beam_crossing()<<"setting event_vertex for trackid " << trackid << " to vtxid " << vtxkey<< std::endl; + std::cout<< "crossing:" << newTrack.get_crossing()<<"trk vtx "<get_beam_crossing()) == short(abs_cross) && vertex->size_tracks()>2) + { + event_vtx(0) = vertex->get_x(); + event_vtx(1) = vertex->get_y(); + event_vtx(2) = vertex->get_z(); + passed_vtx_flag = true; + } + } + } + } + else + { + passed_vtx_flag=true; + } + + if(m_fixed_vtx) + { + event_vtx(0) = m_fixed_vtx_x; + event_vtx(1) = m_fixed_vtx_y; + } + if(!passed_vtx_flag) + { + continue; + } + cumulative_global_vec.push_back(global_vec); cumulative_cluskey_vec.push_back(cluskey_vec); cumulative_vertex.push_back(track_vtx); @@ -530,6 +592,9 @@ int HelicalFitter::process_event(PHCompositeNode* /*unused*/) cumulative_fitpars_mvtx_half_vec.push_back(fitpars_mvtx_half); cumulative_someseed.push_back(someseed); cumulative_newTrack.push_back(newTrack); + cumulative_trackid.push_back(trackid-trkid_correction); + cumulative_event_vertex.push_back(event_vtx); + } // terminate loop over tracks @@ -551,6 +616,8 @@ int HelicalFitter::process_event(PHCompositeNode* /*unused*/) for (unsigned int trackid = 0; trackid < accepted_tracks; ++trackid) { + + const auto& global_vec = cumulative_global_vec[trackid]; const auto& cluskey_vec = cumulative_cluskey_vec[trackid]; auto fitpars = cumulative_fitpars_vec[trackid]; @@ -558,8 +625,10 @@ int HelicalFitter::process_event(PHCompositeNode* /*unused*/) const auto& someseed = cumulative_someseed[trackid]; auto newTrack = cumulative_newTrack[trackid]; SvtxAlignmentStateMap::StateVec statevec; + int trackid_test=cumulative_trackid[trackid]; + auto event_vtx = cumulative_event_vertex[trackid]; + float pull_cumulative = 0.; - // get the residuals and derivatives for all clusters for (unsigned int ivec = 0; ivec < global_vec.size(); ++ivec) { auto global = global_vec[ivec]; @@ -590,6 +659,7 @@ int HelicalFitter::process_event(PHCompositeNode* /*unused*/) else { fitpoint = get_helix_surface_intersection(surf, fitpars, global, helix_pca, helix_tangent); + fitpoint_mvtx_half = get_helix_surface_intersection(surf, fitpars_mvtx_half, global, helix_pca, helix_tangent); } // fitpoint is the point where the helical fit intersects the plane of the surface @@ -609,7 +679,114 @@ int HelicalFitter::process_event(PHCompositeNode* /*unused*/) } Acts::Vector2 residual(xloc - fitpoint_local(0), zloc - fitpoint_local(1)); + unsigned int const layer = TrkrDefs::getLayer(cluskey_vec[ivec]); + + SvtxTrackState_v1 svtxstate(fitpoint.norm()); + svtxstate.set_x(fitpoint(0)); + svtxstate.set_y(fitpoint(1)); + svtxstate.set_z(fitpoint(2)); + std::pair tangent; + if (straight_line_fit) + { + tangent = get_line_tangent(fitpars, global); + } + else + { + tangent = get_helix_tangent(fitpars, global, is_cosmics); + } + svtxstate.set_px(someseed.get_p() * tangent.second.x()); + svtxstate.set_py(someseed.get_p() * tangent.second.y()); + svtxstate.set_pz(someseed.get_p() * tangent.second.z()); + newTrack.insert_state(&svtxstate); + + if (Verbosity() > 1) + { + Acts::Vector3 loc_check = surf->localToGlobalTransform(_tGeometry->geometry().getGeoContext()).inverse() * (global * Acts::UnitConstants::cm); + loc_check /= Acts::UnitConstants::cm; + std::cout << " layer " << layer << std::endl + << " cluster global " << global(0) << " " << global(1) << " " << global(2) << std::endl + << " fitpoint " << fitpoint(0) << " " << fitpoint(1) << " " << fitpoint(2) << std::endl + << " fitpoint_local " << fitpoint_local(0) << " " << fitpoint_local(1) << " " << fitpoint_local(2) << std::endl + << " cluster local x " << cluster->getLocalX() << " cluster local y " << cluster->getLocalY() << std::endl + << " cluster global to local x " << loc_check(0) << " local y " << loc_check(1) << " local z " << loc_check(2) << std::endl + << " cluster local residual x " << residual(0) << " cluster local residual y " << residual(1) << std::endl; + } + + if (Verbosity() > 1) + { + Acts::Transform3 transform = surf->localToGlobalTransform(_tGeometry->geometry().getGeoContext()); + std::cout << "Transform is:" << std::endl; + std::cout << transform.matrix() << std::endl; + Acts::Vector3 loc_check = surf->localToGlobalTransform(_tGeometry->geometry().getGeoContext()).inverse() * (global * Acts::UnitConstants::cm); + loc_check /= Acts::UnitConstants::cm; + unsigned int const sector = TpcDefs::getSectorId(cluskey_vec[ivec]); + unsigned int const side = TpcDefs::getSide(cluskey_vec[ivec]); + std::cout << " layer " << layer << " sector " << sector << " side " << side << " subsurf " << cluster->getSubSurfKey() << std::endl + << " cluster global " << global(0) << " " << global(1) << " " << global(2) << std::endl + << " fitpoint " << fitpoint(0) << " " << fitpoint(1) << " " << fitpoint(2) << std::endl + << " fitpoint_local " << fitpoint_local(0) << " " << fitpoint_local(1) << " " << fitpoint_local(2) << std::endl + << " cluster local x " << cluster->getLocalX() << " cluster local y " << cluster->getLocalY() << std::endl + << " cluster global to local x " << loc_check(0) << " local y " << loc_check(1) << " local z " << loc_check(2) << std::endl + << " cluster local residual x " << residual(0) << " cluster local residual y " << residual(1) << std::endl; + } + + // need standard deviation of measurements + Acts::Vector2 clus_sigma = getClusterError(cluster, cluskey, global); + pull_cumulative += ((residual(0)*residual(0))/(clus_sigma(0)*clus_sigma(0)))+((residual(1)*residual(1))/(clus_sigma(1)*clus_sigma(1))); + } + + // get the residuals and derivatives for all clusters + for (unsigned int ivec = 0; ivec < global_vec.size(); ++ivec) + { + auto global = global_vec[ivec]; + auto cluskey = cluskey_vec[ivec]; + auto cluster = _cluster_map->findCluster(cluskey); + + if (!cluster) + { + continue; + } + + unsigned int const trkrid = TrkrDefs::getTrkrId(cluskey); + + // What we need now is to find the point on the surface at which the helix would intersect + // If we have that point, we can transform the fit back to local coords + // we have fitpars for the helix, and the cluster key - from which we get the surface + + Surface const surf = _tGeometry->maps().getSurface(cluskey, cluster); + Acts::Vector3 helix_pca(0, 0, 0); + Acts::Vector3 helix_tangent(0, 0, 0); + Acts::Vector3 fitpoint; + Acts::Vector3 fitpoint_mvtx_half; + if (straight_line_fit) + { + fitpoint = get_line_surface_intersection(surf, fitpars); + fitpoint_mvtx_half = get_line_surface_intersection(surf, fitpars_mvtx_half); + } + else + { + fitpoint = get_helix_surface_intersection(surf, fitpars, global, helix_pca, helix_tangent); + fitpoint_mvtx_half = get_helix_surface_intersection(surf, fitpars_mvtx_half, global, helix_pca, helix_tangent); + } + + // fitpoint is the point where the helical fit intersects the plane of the surface + // Now transform the helix fitpoint to local coordinates to compare with cluster local coordinates + Acts::Vector3 fitpoint_local = surf->localToGlobalTransform(_tGeometry->geometry().getGeoContext()).inverse() * (fitpoint * Acts::UnitConstants::cm); + Acts::Vector3 fitpoint_mvtx_half_local = surf->localToGlobalTransform(_tGeometry->geometry().getGeoContext()).inverse() * (fitpoint_mvtx_half * Acts::UnitConstants::cm); + + fitpoint_local /= Acts::UnitConstants::cm; + fitpoint_mvtx_half_local /= Acts::UnitConstants::cm; + + auto xloc = cluster->getLocalX(); // in cm + auto zloc = cluster->getLocalY(); + + if (trkrid == TrkrDefs::tpcId) + { + zloc = convertTimeToZ(cluskey, cluster); + } + + Acts::Vector2 residual(xloc - fitpoint_local(0), zloc - fitpoint_local(1)); unsigned int const layer = TrkrDefs::getLayer(cluskey_vec[ivec]); float const phi = atan2(global(1), global(0)); @@ -624,7 +801,7 @@ int HelicalFitter::process_event(PHCompositeNode* /*unused*/) } else { - tangent = get_helix_tangent(fitpars, global); + tangent = get_helix_tangent(fitpars, global, is_cosmics); } svtxstate.set_px(someseed.get_p() * tangent.second.x()); @@ -800,9 +977,9 @@ int HelicalFitter::process_event(PHCompositeNode* /*unused*/) } if (straight_line_fit) { - float ntp_data[78] = { - (float) event, (float) trackid, - (float) layer, (float) nsilicon, (float) h2h_flag, (float) ntpc, (float) nclus, (float) trkrid, (float) sector, (float) side, + float ntp_data[79] = { + (float) event, (float) trackid_test, + (float) layer, (float) nsilicon, (float) pull_cumulative, (float) h2h_flag, (float) ntpc, (float) nclus, (float) trkrid, (float) sector, (float) side, (float) subsurf, phi, (float) glbl_label[0], (float) glbl_label[1], (float) glbl_label[2], (float) glbl_label[3], (float) glbl_label[4], (float) glbl_label[5], (float) sensorCenter(0), (float) sensorCenter(1), (float) sensorCenter(2), @@ -822,13 +999,13 @@ int HelicalFitter::process_event(PHCompositeNode* /*unused*/) lcl_derivativeY[0], lcl_derivativeY[1], lcl_derivativeY[2], lcl_derivativeY[3], glbl_derivativeY[0], glbl_derivativeY[1], glbl_derivativeY[2], glbl_derivativeY[3], glbl_derivativeY[4], glbl_derivativeY[5]}; - ntp->Fill(ntp_data); + ntp->Fill(ntp_data); } else { - float ntp_data[75] = { - (float) event, (float) trackid, - (float) layer, (float) nsilicon, (float) ntpc, (float) nclus, (float) trkrid, (float) sector, (float) side, + float ntp_data[81] = { + (float) event, (float) trackid_test, + (float) layer, (float) nsilicon, (float) h2h_flag, (float) ntpc, (float) nclus, (float) trkrid, (float) pull_cumulative, (float) newTrack.get_charge(), (float) newTrack.get_crossing(), (float) sector, (float) side, (float) subsurf, phi, (float) glbl_label[0], (float) glbl_label[1], (float) glbl_label[2], (float) glbl_label[3], (float) glbl_label[4], (float) glbl_label[5], (float) sensorCenter(0), (float) sensorCenter(1), (float) sensorCenter(2), @@ -841,7 +1018,7 @@ int HelicalFitter::process_event(PHCompositeNode* /*unused*/) (float) fitpoint(0), (float) fitpoint(1), (float) fitpoint(2), (float) tangent.first.x(), (float) tangent.first.y(), (float) tangent.first.z(), (float) tangent.second.x(), (float) tangent.second.y(), (float) tangent.second.z(), - xloc, zloc, (float) fitpoint_local(0), (float) fitpoint_local(1), + xloc, zloc, (float) clus_sigma(0), (float) clus_sigma(1), (float) fitpoint_local(0), (float) fitpoint_local(1), lcl_derivativeX[0], lcl_derivativeX[1], lcl_derivativeX[2], lcl_derivativeX[3], lcl_derivativeX[4], glbl_derivativeX[0], glbl_derivativeX[1], glbl_derivativeX[2], glbl_derivativeX[3], glbl_derivativeX[4], glbl_derivativeX[5], lcl_derivativeY[0], lcl_derivativeY[1], lcl_derivativeY[2], lcl_derivativeY[3], lcl_derivativeY[4], @@ -860,7 +1037,10 @@ int HelicalFitter::process_event(PHCompositeNode* /*unused*/) } } - if (!isnan(residual(0)) && clus_sigma(0) < 1.0) // discards crazy clusters + bool pull_cumulative_pass = true; + if (pull_cumulative>2000) + pull_cumulative_pass = false; + if (!isnan(residual(0)) && clus_sigma(0) < 1.0&&pull_cumulative_pass) // discards crazy clusters { if (arr_has_nan(lcl_derivativeX)) { @@ -872,10 +1052,15 @@ int HelicalFitter::process_event(PHCompositeNode* /*unused*/) std::cerr << "glbl_derivativeX is NaN" << std::endl; continue; } - _mille->mille(AlignmentDefs::NLC, lcl_derivativeX, AlignmentDefs::NGL, glbl_derivativeX, glbl_label, residual(0), errinf * clus_sigma(0)); + if(Verbosity() > 2){ + std::cout << "xyz layer " << layer << " buffers:" << std::endl; + AlignmentDefs::printBuffers(0, residual, clus_sigma, lcl_derivativeX, glbl_derivativeX, glbl_label); + + } + _mille->mille(AlignmentDefs::NLC, lcl_derivativeX, AlignmentDefs::NGL, glbl_derivativeX, glbl_label, residual(0), errinf * clus_sigma(0)); } - if (!isnan(residual(1)) && clus_sigma(1) < 1.0) + if (!isnan(residual(1)) && clus_sigma(1) < 1.0&&pull_cumulative_pass) { if (arr_has_nan(lcl_derivativeY)) { @@ -891,47 +1076,17 @@ int HelicalFitter::process_event(PHCompositeNode* /*unused*/) } } - m_alignmentmap->insertWithKey(trackid, statevec); - m_trackmap->insertWithKey(&newTrack, trackid); + m_alignmentmap->insertWithKey(trackid_test, statevec); + m_trackmap->insertWithKey(&newTrack, trackid_test); // if cosmics, end here, if collision track, continue with vtx // skip the common vertex requirement for this track unless there are 3 tracks in the event - if (accepted_tracks < 3) + if (accepted_tracks < 1) { _mille->end(); continue; } // calculate vertex residual with perigee surface - //------------------------------------------------------- - - Acts::Vector3 event_vtx(averageVertex(0), averageVertex(1), averageVertex(2)); - - if (m_vertexmap) - { - for (const auto& [vtxkey, vertex] : *m_vertexmap) - { - for (auto trackiter = vertex->begin_tracks(); trackiter != vertex->end_tracks(); ++trackiter) - { - SvtxTrack* vtxtrack = m_trackmap->get(*trackiter); - if (vtxtrack) - { - unsigned int const vtxtrackid = vtxtrack->get_id(); - if (trackid == vtxtrackid) - { - event_vtx(0) = vertex->get_x(); - event_vtx(1) = vertex->get_y(); - event_vtx(2) = vertex->get_z(); - if (Verbosity() > 0) - { - std::cout << " setting event_vertex for trackid " << trackid << " to vtxid " << vtxkey - << " vtx " << event_vtx(0) << " " << event_vtx(1) << " " << event_vtx(2) << std::endl; - } - } - } - } - } - } - // The residual for the vtx case is (event vtx - track vtx) // that is -dca float dca3dxy = 0; @@ -949,7 +1104,7 @@ int HelicalFitter::process_event(PHCompositeNode* /*unused*/) // These are local coordinate residuals in the perigee surface Acts::Vector2 vtx_residual(-dca3dxy, -dca3dz); - + float lclvtx_derivativeX[AlignmentDefs::NLC]; float lclvtx_derivativeY[AlignmentDefs::NLC]; if (straight_line_fit) @@ -1004,8 +1159,10 @@ int HelicalFitter::process_event(PHCompositeNode* /*unused*/) std::cout << i << ", "; } } - - if (!isnan(vtx_residual(0))) + bool pull_cumulative_pass = true; + if (pull_cumulative<2000) + pull_cumulative_pass = false; + if (!isnan(vtx_residual(0))&&pull_cumulative_pass) { if (arr_has_nan(lclvtx_derivativeX)) { @@ -1019,7 +1176,7 @@ int HelicalFitter::process_event(PHCompositeNode* /*unused*/) } _mille->mille(AlignmentDefs::NLC, lclvtx_derivativeX, AlignmentDefs::NGLVTX, glblvtx_derivativeX, AlignmentDefs::glbl_vtx_label, vtx_residual(0), vtx_sigma(0)); } - if (!isnan(vtx_residual(1))) + if (!isnan(vtx_residual(1))&&pull_cumulative_pass) { if (arr_has_nan(lclvtx_derivativeY)) { @@ -1042,9 +1199,11 @@ int HelicalFitter::process_event(PHCompositeNode* /*unused*/) float const perigee_phi = atan2(r(1), r(0)); float const track_phi = atan2(newTrack.get_py(), newTrack.get_px()); float const track_eta = atanh(newTrack.get_pz() / newTrack.get_p()); + float const track_p = newTrack.get_p(); + float const track_pt = newTrack.get_pt(); if (straight_line_fit) { - float ntp_data[28] = {(float) trackid, (float) vtx_residual(0), (float) vtx_residual(1), (float) vtx_sigma(0), (float) vtx_sigma(1), + float ntp_data[29] = {(float) trackid_test, (float) pull_cumulative, (float) vtx_residual(0), (float) vtx_residual(1), (float) vtx_sigma(0), (float) vtx_sigma(1), lclvtx_derivativeX[0], lclvtx_derivativeX[1], lclvtx_derivativeX[2], lclvtx_derivativeX[3], glblvtx_derivativeX[0], glblvtx_derivativeX[1], glblvtx_derivativeX[2], lclvtx_derivativeY[0], lclvtx_derivativeY[1], lclvtx_derivativeY[2], lclvtx_derivativeY[3], @@ -1056,13 +1215,13 @@ int HelicalFitter::process_event(PHCompositeNode* /*unused*/) } else { - float ntp_data[29] = {(float) trackid, (float) vtx_residual(0), (float) vtx_residual(1), (float) vtx_sigma(0), (float) vtx_sigma(1), + float ntp_data[34] = {(float) trackid_test, (float) pull_cumulative, (float) vtx_residual(0), (float) vtx_residual(1), (float) vtx_sigma(0), (float) vtx_sigma(1), lclvtx_derivativeX[0], lclvtx_derivativeX[1], lclvtx_derivativeX[2], lclvtx_derivativeX[3], lclvtx_derivativeX[4], glblvtx_derivativeX[0], glblvtx_derivativeX[1], glblvtx_derivativeX[2], lclvtx_derivativeY[0], lclvtx_derivativeY[1], lclvtx_derivativeY[2], lclvtx_derivativeY[3], lclvtx_derivativeY[4], glblvtx_derivativeY[0], glblvtx_derivativeY[1], glblvtx_derivativeY[2], newTrack.get_x(), newTrack.get_y(), newTrack.get_z(), - (float) event_vtx(0), (float) event_vtx(1), (float) event_vtx(2), track_phi, perigee_phi}; + (float) event_vtx(0), (float) event_vtx(1), (float) event_vtx(2), pull_cumulative, track_phi, perigee_phi, track_eta, track_p, track_pt}; track_ntp->Fill(ntp_data); } @@ -1100,7 +1259,7 @@ Acts::Vector3 HelicalFitter::get_helix_surface_intersection(const Surface& surf, // there are analytic solutions for a line-plane intersection. // to use this, need to get the vector tangent to the helix near the measurement and a point on it. - std::pair const line = get_helix_tangent(fitpars, std::move(global)); + std::pair const line = get_helix_tangent(fitpars, std::move(global), is_cosmics); Acts::Vector3 const pca = line.first; Acts::Vector3 const tangent = line.second; @@ -1157,7 +1316,7 @@ Acts::Vector3 HelicalFitter::get_helix_surface_intersection(const Surface& surf, // there are analytic solutions for a line-plane intersection. // to use this, need to get the vector tangent to the helix near the measurement and a point on it. - std::pair const line = get_helix_tangent(fitpars, std::move(global)); + std::pair const line = get_helix_tangent(fitpars, std::move(global), is_cosmics); pca = line.first; tangent = line.second; @@ -1276,20 +1435,12 @@ Acts::Vector3 HelicalFitter::get_line_plane_intersection(const Acts::Vector3& PC float const d = (sensor_center - PCA).dot(sensor_normal) / tangent.dot(sensor_normal); Acts::Vector3 intersection = PCA + d * tangent; - /* - std::cout << " sensor center " << sensor_center(0) << " " << sensor_center(1) << " " << sensor_center(2) << std::endl; - std::cout << " intersection " << intersection(0) << " " << intersection(1) << " " << intersection(2) << std::endl; - std::cout << " PCA " << PCA(0) << " " << PCA(1) << " " << PCA(2) << std::endl; - std::cout << " tangent " << tangent(0) << " " << tangent(1) << " " << tangent(2) << std::endl; - std::cout << " d " << d << std::endl; - */ - return intersection; } -std::pair HelicalFitter::get_helix_tangent(const std::vector& fitpars, Acts::Vector3 global) +std::pair HelicalFitter::get_helix_tangent(const std::vector& fitpars, Acts::Vector3 global, bool is_cosmics) { - auto pair = TrackFitUtils::get_helix_tangent(fitpars, global); + auto pair = TrackFitUtils::get_helix_tangent(fitpars, global, is_cosmics); /* save for posterity purposes if(Verbosity() > 2) @@ -1501,8 +1652,8 @@ void HelicalFitter::getTrackletClusterList(TrackSeed* tracklet, std::vector 2 && layer < 7) + //// drop INTT clusters for now -- TEMPORARY! + // if (layer > 2 && layer < 5) //{ // continue; //} @@ -1525,7 +1676,7 @@ Acts::Vector2 HelicalFitter::getClusterError(TrkrCluster* cluster, TrkrDefs::clu auto para_errors = _ClusErrPara.get_clusterv5_modified_error(cluster, clusRadius, cluskey); double const phierror = sqrt(para_errors.first); double const zerror = sqrt(para_errors.second); - clus_sigma(1) = zerror; + clus_sigma(1) = zerror * 2.0; clus_sigma(0) = phierror; return clus_sigma; @@ -1555,7 +1706,7 @@ void HelicalFitter::getLocalDerivativesXY(const Surface& surf, const Acts::Vecto { std::cout << "Call get_helix_tangent for best fit fitpars" << std::endl; } - std::pair const tangent = get_helix_tangent(fitpars, global); + std::pair const tangent = get_helix_tangent(fitpars, global, is_cosmics); Acts::Vector3 projX(0, 0, 0), projY(0, 0, 0); get_projectionXY(surf, tangent, projX, projY); @@ -1808,7 +1959,7 @@ void HelicalFitter::getGlobalDerivativesXY(const Surface& surf, const Acts::Vect } else { - tangent = get_helix_tangent(fitpars, global); + tangent = get_helix_tangent(fitpars, global, is_cosmics); } Acts::Vector3 projX(0, 0, 0), projY(0, 0, 0); diff --git a/offline/packages/TrackerMillepedeAlignment/HelicalFitter.h b/offline/packages/TrackerMillepedeAlignment/HelicalFitter.h index 12fe31ccb5..ae7c6187c8 100644 --- a/offline/packages/TrackerMillepedeAlignment/HelicalFitter.h +++ b/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 get_helix_tangent(const std::vector& fitpars, Acts::Vector3 global); + std::pair get_helix_tangent(const std::vector& fitpars, Acts::Vector3 global, bool is_cosmics); Acts::Vector3 get_helix_surface_intersection(const Surface& surf, std::vector& fitpars, Acts::Vector3 global); Acts::Vector3 get_helix_surface_intersection(const Surface& surf, std::vector& 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}; diff --git a/offline/packages/TrackerMillepedeAlignment/MakeMilleFiles.cc b/offline/packages/TrackerMillepedeAlignment/MakeMilleFiles.cc index cb50fd55f0..e93b9cb23a 100644 --- a/offline/packages/TrackerMillepedeAlignment/MakeMilleFiles.cc +++ b/offline/packages/TrackerMillepedeAlignment/MakeMilleFiles.cc @@ -5,6 +5,7 @@ /// Tracking includes #include +#include #include // for side #include #include @@ -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): "<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; } + + diff --git a/offline/packages/TrackerMillepedeAlignment/MakeMilleFiles.h b/offline/packages/TrackerMillepedeAlignment/MakeMilleFiles.h index 0733e22bb6..4b507141bb 100644 --- a/offline/packages/TrackerMillepedeAlignment/MakeMilleFiles.h +++ b/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> fixed_mvtx_layers; std::set> fixed_layer_gparams, fixed_layer_lparams; + bool m_ignore_tpc = false; + std::string m_constraintFileName = "mp2con.txt"; std::ofstream m_constraintFile; diff --git a/offline/packages/TrackingDiagnostics/TrackResiduals.cc b/offline/packages/TrackingDiagnostics/TrackResiduals.cc index 7ca8c6323e..3953538da3 100644 --- a/offline/packages/TrackingDiagnostics/TrackResiduals.cc +++ b/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"); diff --git a/offline/packages/TrackingDiagnostics/TrackResiduals.h b/offline/packages/TrackingDiagnostics/TrackResiduals.h index c14085c2b9..c0749a76a4 100644 --- a/offline/packages/TrackingDiagnostics/TrackResiduals.h +++ b/offline/packages/TrackingDiagnostics/TrackResiduals.h @@ -243,6 +243,12 @@ class TrackResiduals : public SubsysReco int m_ntracks = std::numeric_limits::quiet_NaN(); int m_nvertices = std::numeric_limits::quiet_NaN(); + std::vector m_pcax_vtx_trk; + std::vector m_pcay_vtx_trk; + std::vector m_pcaz_vtx_trk; + std::vector m_px_vtx_trk; + std::vector m_py_vtx_trk; + std::vector m_pz_vtx_trk; //! cluster tree info float m_sclusgr = std::numeric_limits::quiet_NaN(); diff --git a/offline/packages/TrackingDiagnostics/TrackSeedTrackMapConverter.cc b/offline/packages/TrackingDiagnostics/TrackSeedTrackMapConverter.cc index 969e746c82..df9e7cb27f 100644 --- a/offline/packages/TrackingDiagnostics/TrackSeedTrackMapConverter.cc +++ b/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 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()); } diff --git a/offline/packages/trackbase/TrackFitUtils.cc b/offline/packages/trackbase/TrackFitUtils.cc index e6ecbb16ab..8773b61651 100644 --- a/offline/packages/trackbase/TrackFitUtils.cc +++ b/offline/packages/trackbase/TrackFitUtils.cc @@ -37,7 +37,7 @@ namespace } } // namespace -std::pair TrackFitUtils::get_helix_tangent(const std::vector& fitpars, Acts::Vector3& global) +std::pair TrackFitUtils::get_helix_tangent(const std::vector& fitpars, Acts::Vector3& global, bool is_cosmics) { // no analytic solution for the coordinates of the closest approach of a helix to a point // Instead, we get the PCA in x and y to the circle, and the PCA in z to the z vs R line at the R of the PCA @@ -52,7 +52,12 @@ std::pair TrackFitUtils::get_helix_tangent(const s // The radius of the PCA determines the z position: float const pca_circle_radius = pca_circle.norm(); // radius of the PCA of the circle to the point - float const pca_z = pca_circle_radius * zslope + z0; + float ztmp; + if(is_cosmics) + ztmp = pca_circle(0)* zslope + z0; + else + ztmp = pca_circle_radius * zslope + z0; + float const pca_z = ztmp; Acts::Vector3 const pca(pca_circle(0), pca_circle(1), pca_z); // now we want a second point on the helix so we can get a local straight line approximation to the track @@ -62,7 +67,12 @@ std::pair TrackFitUtils::get_helix_tangent(const s float const d_angle = 0.005; float const newx = radius * std::cos(angle_pca + d_angle) + x0; float const newy = radius * std::sin(angle_pca + d_angle) + y0; - float const newz = std::sqrt(newx * newx + newy * newy) * zslope + z0; + float ztmp2; + if(is_cosmics) + ztmp2 = newx * zslope + z0; + else + ztmp2 = std::sqrt(newx * newx + newy * newy) * zslope + z0; + float const newz = ztmp2; Acts::Vector3 const second_point_pca(newx, newy, newz); // pca and second_point_pca define a straight line approximation to the track @@ -557,7 +567,7 @@ unsigned int TrackFitUtils::addClusters(std::vector& fitpars, //_________________________________________________________________________________ Acts::Vector3 TrackFitUtils::get_helix_pca(std::vector& fitpars, - const Acts::Vector3& global) + const Acts::Vector3& global, bool is_cosmics) { // no analytic solution for the coordinates of the closest approach of a helix to a point // Instead, we get the PCA in x and y to the circle, and the PCA in z to the z vs R line at the R of the PCA @@ -572,7 +582,12 @@ Acts::Vector3 TrackFitUtils::get_helix_pca(std::vector& fitpars, // The radius of the PCA determines the z position: float const pca_circle_radius = pca_circle.norm(); - float const pca_z = pca_circle_radius * zslope + z0; + float ztmp; + if(is_cosmics) + ztmp = pca_circle(0)* zslope + z0; + else + ztmp = pca_circle_radius * zslope + z0; + float const pca_z = ztmp; Acts::Vector3 const pca(pca_circle(0), pca_circle(1), pca_z); // now we want a second point on the helix so we can get a local straight line approximation to the track @@ -580,7 +595,12 @@ Acts::Vector3 TrackFitUtils::get_helix_pca(std::vector& fitpars, float const projection = 0.25; // cm Acts::Vector3 const second_point = pca + projection * pca / pca.norm(); Acts::Vector2 second_point_pca_circle = get_circle_point_pca(radius, x0, y0, second_point); - float const second_point_pca_z = second_point_pca_circle.norm() * zslope + z0; + float ztmp2; + if(is_cosmics) + ztmp2 = second_point_pca_circle(0)* zslope + z0; + else + ztmp2 = second_point_pca_circle.norm() * zslope + z0; + float const second_point_pca_z = ztmp2; Acts::Vector3 const second_point_pca(second_point_pca_circle(0), second_point_pca_circle(1), second_point_pca_z); // pca and second_point_pca define a straight line approximation to the track @@ -610,7 +630,7 @@ Acts::Vector2 TrackFitUtils::get_circle_point_pca(float radius, float x0, float //_________________________________________________________________________________ std::vector TrackFitUtils::fitClusters(std::vector& global_vec, const std::vector &cluskey_vec, - bool use_intt) + bool use_intt, bool mvtx_east, bool mvtx_west, bool is_cosmics) { std::vector fitpars; @@ -619,15 +639,32 @@ std::vector TrackFitUtils::fitClusters(std::vector& global { return fitpars; } - std::tuple circle_fit_pars = TrackFitUtils::circle_fit_by_taubin(global_vec); + //std::tuple circle_fit_pars = TrackFitUtils::circle_fit_by_taubin(global_vec); + std::tuple circle_fit_pars ; + bool cross_mvtx_half = false; + if ((mvtx_east || mvtx_west)) + { + cross_mvtx_half = TrackFitUtils::isTrackCrossMvtxHalf(cluskey_vec); + } + else + { + circle_fit_pars = TrackFitUtils::circle_fit_by_taubin(global_vec); + } // It is problematic that the large errors on the INTT strip z values are not allowed for - drop the INTT from the z line fit std::vector global_vec_noINTT; for (unsigned int ivec = 0; ivec < global_vec.size(); ++ivec) { unsigned int const trkrid = TrkrDefs::getTrkrId(cluskey_vec[ivec]); + if ((mvtx_east || mvtx_west) && trkrid == TrkrDefs::mvtxId) + { + if (cross_mvtx_half && TrackFitUtils::includeMvtxHit(cluskey_vec[ivec], mvtx_east, mvtx_west)) + { + global_vec_noINTT.push_back(global_vec[ivec]); + } + } - if (trkrid != TrkrDefs::inttId && cluskey_vec[ivec] != 0) + else if (trkrid != TrkrDefs::inttId and cluskey_vec[ivec] != 0) { global_vec_noINTT.push_back(global_vec[ivec]); } @@ -641,7 +678,11 @@ std::vector TrackFitUtils::fitClusters(std::vector& global { return fitpars; } - std::tuple line_fit_pars = TrackFitUtils::line_fit(global_vec_noINTT); + std::tuple line_fit_pars; + if(is_cosmics) + line_fit_pars = TrackFitUtils::line_fit_xz(global_vec_noINTT); + else + line_fit_pars = TrackFitUtils::line_fit(global_vec_noINTT); fitpars.push_back(std::get<0>(circle_fit_pars)); fitpars.push_back(std::get<1>(circle_fit_pars)); diff --git a/offline/packages/trackbase/TrackFitUtils.h b/offline/packages/trackbase/TrackFitUtils.h index ae9dc22481..24aa14bc34 100644 --- a/offline/packages/trackbase/TrackFitUtils.h +++ b/offline/packages/trackbase/TrackFitUtils.h @@ -119,15 +119,16 @@ namespace TrackFitUtils const unsigned int& startLayer, const unsigned int& endLayer); - std::pair get_helix_tangent(const std::vector& fitpars, Acts::Vector3& global); + std::pair get_helix_tangent(const std::vector& fitpars, Acts::Vector3& global, bool is_cosmics=false); - Acts::Vector3 get_helix_pca(std::vector& fitpars, const Acts::Vector3& global); + Acts::Vector3 get_helix_pca(std::vector& fitpars, const Acts::Vector3& global, bool is_cosmics= false); Acts::Vector2 get_circle_point_pca(float radius, float x0, float y0, Acts::Vector3 global); std::vector fitClusters(std::vector& global_vec, const std::vector &cluskey_vec, - bool use_intt = false); + bool use_intt = false, bool mvtx_east = false, bool mvtx_west = false, bool is_cosmics=false); + void getTrackletClusters(ActsGeometry* _tGeometry, TrkrClusterContainer* _cluster_map, std::vector& global_vec, diff --git a/offline/packages/trackbase_historic/SvtxAlignmentState.cc b/offline/packages/trackbase_historic/SvtxAlignmentState.cc index e9c0d2cab0..281c504305 100644 --- a/offline/packages/trackbase_historic/SvtxAlignmentState.cc +++ b/offline/packages/trackbase_historic/SvtxAlignmentState.cc @@ -4,7 +4,10 @@ namespace { SvtxAlignmentState::GlobalMatrix globalMatrix = SvtxAlignmentState::GlobalMatrix::Zero(); SvtxAlignmentState::LocalMatrix localMatrix = SvtxAlignmentState::LocalMatrix::Zero(); + SvtxAlignmentState::LocalMatrixPsuedo localMatrixPsuedo = SvtxAlignmentState::LocalMatrixPsuedo::Zero(); + SvtxAlignmentState::LocalMeasErrPsuedo localMeasErrPsuedo = SvtxAlignmentState::LocalMeasErrPsuedo::Zero(); SvtxAlignmentState::ResidualVector residual = SvtxAlignmentState::ResidualVector::Zero(); + SvtxAlignmentState::ActsTrackParamsVector trackParams = SvtxAlignmentState::ActsTrackParamsVector::Zero(); } // namespace const SvtxAlignmentState::ResidualVector& SvtxAlignmentState::get_residual() const @@ -17,7 +20,22 @@ const SvtxAlignmentState::LocalMatrix& SvtxAlignmentState::get_local_derivative_ return localMatrix; } +const SvtxAlignmentState::LocalMatrixPsuedo& SvtxAlignmentState::get_local_derivative_psuedo_matrix() const +{ + return localMatrixPsuedo; +} + +const SvtxAlignmentState::LocalMeasErrPsuedo& SvtxAlignmentState::get_local_psuedo_measurement_err() const +{ + return localMeasErrPsuedo; +} + const SvtxAlignmentState::GlobalMatrix& SvtxAlignmentState::get_global_derivative_matrix() const { return globalMatrix; } + +const SvtxAlignmentState::ActsTrackParamsVector& SvtxAlignmentState::get_acts_track_params() const +{ + return trackParams; +} diff --git a/offline/packages/trackbase_historic/SvtxAlignmentState.h b/offline/packages/trackbase_historic/SvtxAlignmentState.h index a1b207800c..1adc36a44c 100644 --- a/offline/packages/trackbase_historic/SvtxAlignmentState.h +++ b/offline/packages/trackbase_historic/SvtxAlignmentState.h @@ -20,7 +20,10 @@ class SvtxAlignmentState : public PHObject typedef Eigen::Matrix GlobalMatrix; typedef Eigen::Matrix LocalMatrix; + typedef Eigen::Matrix LocalMatrixPsuedo; + typedef Eigen::Matrix LocalMeasErrPsuedo; typedef Eigen::Matrix ResidualVector; + typedef Eigen::Matrix ActsTrackParamsVector; ~SvtxAlignmentState() override {} @@ -34,13 +37,21 @@ class SvtxAlignmentState : public PHObject virtual void set_residual(const ResidualVector&) {} virtual void set_local_derivative_matrix(const LocalMatrix&) {} + //Local derivative matrix psuedo measurement for ACTS -> Millepede + virtual void set_local_derivative_psuedo_matrix(const LocalMatrixPsuedo&) {} + //setter for eigen vals only called in above func for psuedo matrix + virtual void set_local_psuedo_measurement_err(const LocalMeasErrPsuedo&) {} virtual void set_global_derivative_matrix(const GlobalMatrix&) {} virtual void set_cluster_key(TrkrDefs::cluskey) {} + virtual void set_acts_track_params(const ActsTrackParamsVector&) {} virtual const ResidualVector& get_residual() const; virtual const LocalMatrix& get_local_derivative_matrix() const; + virtual const LocalMatrixPsuedo& get_local_derivative_psuedo_matrix() const; + virtual const LocalMeasErrPsuedo& get_local_psuedo_measurement_err() const; virtual const GlobalMatrix& get_global_derivative_matrix() const; virtual TrkrDefs::cluskey get_cluster_key() const { return UINT_MAX; } + virtual const ActsTrackParamsVector& get_acts_track_params() const; protected: SvtxAlignmentState() {} diff --git a/offline/packages/trackbase_historic/SvtxAlignmentState_v1.cc b/offline/packages/trackbase_historic/SvtxAlignmentState_v1.cc index 28e8d2c5f2..70fc083b4d 100644 --- a/offline/packages/trackbase_historic/SvtxAlignmentState_v1.cc +++ b/offline/packages/trackbase_historic/SvtxAlignmentState_v1.cc @@ -5,6 +5,7 @@ SvtxAlignmentState_v1::SvtxAlignmentState_v1() , m_localDeriv(LocalMatrix::Zero()) , m_globalDeriv(GlobalMatrix::Zero()) , m_cluskey(UINT_MAX) + , m_trackParams(ActsTrackParamsVector::Zero()) { } diff --git a/offline/packages/trackbase_historic/SvtxAlignmentState_v1.h b/offline/packages/trackbase_historic/SvtxAlignmentState_v1.h index 0bb10850a3..2e3ad95ac8 100644 --- a/offline/packages/trackbase_historic/SvtxAlignmentState_v1.h +++ b/offline/packages/trackbase_historic/SvtxAlignmentState_v1.h @@ -35,17 +35,23 @@ class SvtxAlignmentState_v1 : public SvtxAlignmentState { m_cluskey = key; } + void set_acts_track_params(const ActsTrackParamsVector& p) override + { + m_trackParams = p; + } const ResidualVector& get_residual() const override { return m_residual; } const LocalMatrix& get_local_derivative_matrix() const override { return m_localDeriv; } const GlobalMatrix& get_global_derivative_matrix() const override { return m_globalDeriv; } TrkrDefs::cluskey get_cluster_key() const override { return m_cluskey; } + const ActsTrackParamsVector& get_acts_track_params() const override { return m_trackParams; } private: ResidualVector m_residual; LocalMatrix m_localDeriv; GlobalMatrix m_globalDeriv; TrkrDefs::cluskey m_cluskey; + ActsTrackParamsVector m_trackParams; ClassDefOverride(SvtxAlignmentState_v1, 1) }; diff --git a/offline/packages/trackreco/ActsAlignmentStates.cc b/offline/packages/trackreco/ActsAlignmentStates.cc index 070f78e53f..945acc2655 100644 --- a/offline/packages/trackreco/ActsAlignmentStates.cc +++ b/offline/packages/trackreco/ActsAlignmentStates.cc @@ -173,8 +173,14 @@ void ActsAlignmentStates::fillAlignmentStateMap( // Get the derivative of alignment (global) parameters w.r.t. measurement or residual /// The local bound parameters still have access to global phi/theta + const double l0 = state.smoothed()[Acts::eBoundLoc0]; + const double l1 = state.smoothed()[Acts::eBoundLoc1]; const double phi = state.smoothed()[Acts::eBoundPhi]; const double theta = state.smoothed()[Acts::eBoundTheta]; + const double qoverp = state.smoothed()[Acts::eBoundQOverP]; + const double time = state.smoothed()[Acts::eBoundTime]; + + SvtxAlignmentState::ActsTrackParamsVector track_params = {l0,l1,phi,theta,qoverp,time}; Acts::Vector3 tangent = Acts::makeDirectionFromPhiTheta(phi,theta); @@ -206,17 +212,23 @@ void ActsAlignmentStates::fillAlignmentStateMap( //! this is the derivative of the state wrt to Acts track parameters //! e.g. (d_0, z_0, phi, theta, q/p, t) - auto localDeriv = H * state.jacobian(); + //auto localDeriv = H * state.jacobian(); + auto localDeriv = makeLocalDerivatives(H); + SvtxAlignmentState::LocalMeasErrPsuedo localmeaserrpsuedo = SvtxAlignmentState::LocalMeasErrPsuedo::Zero(); + auto localDerivPsuedo = makeLocalDerivativesPsuedo(state, localmeaserrpsuedo); if (m_verbosity > 2) { - std::cout << "local deriv " << std::endl << localDeriv << std::endl; + std::cout << "local deriv " << std::endl << localDeriv << std::endl; } auto svtxstate = std::make_unique(); svtxstate->set_residual(localResidual); svtxstate->set_local_derivative_matrix(localDeriv); + svtxstate->set_local_derivative_psuedo_matrix(localDerivPsuedo); + svtxstate->set_local_psuedo_measurement_err(localmeaserrpsuedo); svtxstate->set_global_derivative_matrix(globDeriv); + svtxstate->set_acts_track_params(track_params); svtxstate->set_cluster_key(ckey); statevec.push_back(svtxstate.release()); @@ -261,6 +273,79 @@ ActsAlignmentStates::makeGlobalDerivatives(const Acts::Vector3& OM, return globalder; } +SvtxAlignmentState::LocalMatrix +ActsAlignmentStates::makeLocalDerivatives(const auto& H) +{ + SvtxAlignmentState::LocalMatrix localder = SvtxAlignmentState::LocalMatrix::Zero(); + //std::cout<< "H is " << std::endl << H << std::endl; + for (int i = 0; i < SvtxAlignmentState::NRES; ++i) + { + for (int j = 0; j < SvtxAlignmentState::NLOC; ++j) + { + localder(i, j) = H(i, j); + } + } + return localder; +} + +SvtxAlignmentState::LocalMatrixPsuedo +ActsAlignmentStates::makeLocalDerivativesPsuedo(const auto& state, SvtxAlignmentState::LocalMeasErrPsuedo& localmeaserrpsuedo) +{ + SvtxAlignmentState::LocalMatrixPsuedo localderpsuedo = SvtxAlignmentState::LocalMatrixPsuedo::Zero(); + localmeaserrpsuedo(0) = 1.; + Acts::ActsDynamicMatrix updatedCovariance{SvtxAlignmentState::NLOC, + SvtxAlignmentState::NLOC}; + Acts::ActsDynamicMatrix updatedProjection{SvtxAlignmentState::NRES, + SvtxAlignmentState::NLOC}; + + const auto H = state.projectorSubspaceHelper().fullProjector().topLeftCorner( + state.calibratedSize(), Acts::eBoundSize); + + for (int internalTPindex = 0; internalTPindex < SvtxAlignmentState::NLOC; ++internalTPindex) + { + updatedProjection.col(internalTPindex) = H.col(internalTPindex); + for (int internalTPindex2 = 0; internalTPindex2 < SvtxAlignmentState::NLOC; ++internalTPindex2) + { + updatedCovariance(internalTPindex, internalTPindex2) = state.covariance()(internalTPindex, internalTPindex2); + } + } + + const Acts::ActsDynamicMatrix weightMatMeasurements = + updatedProjection.transpose() * state.effectiveCalibratedCovariance().inverse() * + updatedProjection; + + const Acts::ActsDynamicMatrix regularisedCov = + regulariseCovariance(updatedCovariance, -1, -1, 1.e-9); + + const Acts::ActsDynamicMatrix correlationTerm = + getInverseComplement(regularisedCov, weightMatMeasurements); + + + Eigen::SelfAdjointEigenSolver eigenSolver( + correlationTerm); + if (eigenSolver.info() != Eigen::Success) { + std::cout << " FAILED to find decompose correlation term" << std::endl; + return localderpsuedo; + } + const Acts::ActsDynamicVector eigenVals = eigenSolver.eigenvalues(); + const Acts::ActsDynamicMatrix eigenVecs = eigenSolver.eigenvectors(); + + /// convert each EV to a pseudo-measurement + for (long iMeas = 0; iMeas < eigenVecs.rows(); + ++iMeas) { // fill the local derivatives from the current eigenvector + // skip negative EV + if (eigenVals(iMeas) <= 0) { + continue; + } + localmeaserrpsuedo(iMeas) = 1./std::sqrt(eigenVals(iMeas)); + for (std::size_t iPar = 0; iPar < SvtxAlignmentState::NLOC; ++iPar) { + localderpsuedo(iPar, iMeas) = eigenVecs(iPar, iMeas); + } + } + + return localderpsuedo; +} + //______________________________________________________ std::pair ActsAlignmentStates::get_projectionXY(const Acts::Surface& surface, const Acts::Vector3& tangent) { @@ -293,3 +378,65 @@ std::pair ActsAlignmentStates::get_projectionXY(co return std::make_pair(projx, projy); } + +const Acts::ActsDynamicMatrix ActsAlignmentStates::regulariseCovariance(const Acts::ActsDynamicMatrix& inputCov, + double conditionCutOff, + double removeHugeLeading, + double stabilisationDiag) +{ + Acts::ActsDynamicMatrix out = + Acts::ActsDynamicMatrix::Zero(inputCov.rows(), inputCov.cols()); + + /// optionally add a tiny diagonal matrix for additional stabilisation + Acts::ActsDynamicMatrix regularisation = + Acts::ActsDynamicMatrix::Zero(inputCov.rows(), inputCov.cols()); + if (stabilisationDiag > 0) { + regularisation = stabilisationDiag * Acts::ActsDynamicMatrix::Identity( + inputCov.rows(), inputCov.cols()); + } + auto eigensolver = Eigen::SelfAdjointEigenSolver( + inputCov + regularisation); + + if (eigensolver.info() != Eigen::Success) { + std::cout << " FAILED to find eigenvec" << std::endl; + return out; + } + auto eigenVals = eigensolver.eigenvalues(); + auto eigenVecs = eigensolver.eigenvectors(); + + std::size_t maxIndex = eigenVals.size() - 1; + double lambdaMax = eigenVals(eigenVals.size() - 1); + // check for a huge leading eigenvalue - this happens when the time coordinate + // is unconstrained + if (removeHugeLeading > 0 && + lambdaMax > removeHugeLeading * eigenVals(maxIndex - 1)) { + --maxIndex; + lambdaMax = eigenVals(maxIndex); + } + double lambdaMin = eigenVals(0); + if (conditionCutOff > 0) { + lambdaMin = conditionCutOff * lambdaMax; + } + // clamp the EV to the permitted interval + for (Eigen::Index i = 0; i < eigenVals.size(); ++i) { + out(i, i) = std::clamp(eigenVals(i), lambdaMin, lambdaMax); + } + // then return the regularised matrix as V D V^T + out = out * eigenVecs.transpose(); + out = eigenVecs * out; + return out; +} + +const Acts::ActsDynamicMatrix ActsAlignmentStates::getInverseComplement(const Acts::ActsDynamicMatrix& target, + const Acts::ActsDynamicMatrix& existing_sol) +{ + Acts::ActsDynamicMatrix Rhs = + Acts::ActsDynamicMatrix::Identity(target.rows(), target.cols()) - + target * existing_sol; + // call decomposition from Eigen (prefer over llt for semi-def. matrices) + auto LDLT = target.ldlt(); + Acts::ActsDynamicMatrix solution = LDLT.solve(Rhs); + // finally, symmetrise to correct for floating point effects + solution = 0.5 * (solution + solution.transpose()); + return solution; +} \ No newline at end of file diff --git a/offline/packages/trackreco/ActsAlignmentStates.h b/offline/packages/trackreco/ActsAlignmentStates.h index 9d891c906d..edd47e6f27 100644 --- a/offline/packages/trackreco/ActsAlignmentStates.h +++ b/offline/packages/trackreco/ActsAlignmentStates.h @@ -55,9 +55,15 @@ class ActsAlignmentStates private: std::pair get_projectionXY(const Acts::Surface& surface, const Acts::Vector3& tangent); - + const Acts::ActsDynamicMatrix regulariseCovariance(const Acts::ActsDynamicMatrix& inputCov, + double conditionCutOff, + double removeHugeLeading, + double stabilisationDiag); + const Acts::ActsDynamicMatrix getInverseComplement(const Acts::ActsDynamicMatrix& target, + const Acts::ActsDynamicMatrix& existing_sol); SvtxAlignmentState::GlobalMatrix makeGlobalDerivatives(const Acts::Vector3& OM, const std::pair& projxy); - + SvtxAlignmentState::LocalMatrix makeLocalDerivatives(const auto& H); + SvtxAlignmentState::LocalMatrixPsuedo makeLocalDerivativesPsuedo(const auto& state, SvtxAlignmentState::LocalMeasErrPsuedo& localmeaserrpsuedo); //! verbosity int m_verbosity = 0; diff --git a/offline/packages/trackreco/MakeSourceLinks.cc b/offline/packages/trackreco/MakeSourceLinks.cc index 58b1b2a3fb..8191e4670a 100644 --- a/offline/packages/trackreco/MakeSourceLinks.cc +++ b/offline/packages/trackreco/MakeSourceLinks.cc @@ -327,7 +327,8 @@ SourceLinkVec MakeSourceLinks::getSourceLinksClusterMover( TrkrClusterContainer* clusterContainer, ActsGeometry* tGeometry, const TpcGlobalPositionWrapper& globalPositionWrapper, - short int crossing) + short int crossing, + bool use_modified_clus_error) { if (m_verbosity > 1) { @@ -525,12 +526,22 @@ SourceLinkVec MakeSourceLinks::getSourceLinksClusterMover( Acts::ActsSquareMatrix<2> cov = Acts::ActsSquareMatrix<2>::Zero(); - double clusRadius = sqrt(global[0] * global[0] + global[1] * global[1]); - auto para_errors = ClusterErrorPara::get_clusterv5_modified_error(cluster, clusRadius, cluskey); - cov(Acts::eBoundLoc0, Acts::eBoundLoc0) = para_errors.first * Acts::UnitConstants::cm2; - cov(Acts::eBoundLoc0, Acts::eBoundLoc1) = 0; - cov(Acts::eBoundLoc1, Acts::eBoundLoc0) = 0; - cov(Acts::eBoundLoc1, Acts::eBoundLoc1) = para_errors.second * Acts::UnitConstants::cm2; + if(use_modified_clus_error) + { + double clusRadius = sqrt(global[0] * global[0] + global[1] * global[1]); + auto para_errors = ClusterErrorPara::get_clusterv5_modified_error(cluster, clusRadius, cluskey); + cov(Acts::eBoundLoc0, Acts::eBoundLoc0) = para_errors.first * Acts::UnitConstants::cm2; + cov(Acts::eBoundLoc0, Acts::eBoundLoc1) = 0; + cov(Acts::eBoundLoc1, Acts::eBoundLoc0) = 0; + cov(Acts::eBoundLoc1, Acts::eBoundLoc1) = para_errors.second * Acts::UnitConstants::cm2; + } + else + { + cov(Acts::eBoundLoc0, Acts::eBoundLoc0) = cluster->getRPhiError()*cluster->getRPhiError() * Acts::UnitConstants::cm2; + cov(Acts::eBoundLoc0, Acts::eBoundLoc1) = 0; + cov(Acts::eBoundLoc1, Acts::eBoundLoc0) = 0; + cov(Acts::eBoundLoc1, Acts::eBoundLoc1) = cluster->getZError()*cluster->getZError() * Acts::UnitConstants::cm2; + } ActsSourceLink::Index index = measurements.size(); diff --git a/offline/packages/trackreco/MakeSourceLinks.h b/offline/packages/trackreco/MakeSourceLinks.h index bf0788a097..3239f7a241 100644 --- a/offline/packages/trackreco/MakeSourceLinks.h +++ b/offline/packages/trackreco/MakeSourceLinks.h @@ -72,7 +72,8 @@ class MakeSourceLinks TrkrClusterContainer* /*clusters*/, ActsGeometry* /*geometry*/, const TpcGlobalPositionWrapper& /*globalpositionWrapper*/, - short int crossing); + short int crossing, + bool use_modified_clus_error = true); private: diff --git a/offline/packages/trackreco/Makefile.am b/offline/packages/trackreco/Makefile.am index 7c34ef2652..9801e6dcba 100644 --- a/offline/packages/trackreco/Makefile.am +++ b/offline/packages/trackreco/Makefile.am @@ -55,6 +55,7 @@ pkginclude_HEADERS = \ PHInitVertexing.h \ PHMicromegasTpcTrackMatching.h \ PHSiliconTpcTrackMatching.h \ + PHSiliconTpcTrackMatchingDummy.h \ PHSiliconCosmicSeeding.h \ PHSimpleVertexFinder.h \ PHTrackPruner.h \ @@ -146,6 +147,7 @@ libtrack_reco_la_SOURCES = \ PHSiliconHelicalPropagator.cc \ PHSiliconSeedMerger.cc \ PHSiliconTpcTrackMatching.cc \ + PHSiliconTpcTrackMatchingDummy.cc \ PHSiliconCosmicSeeding.cc \ PHSimpleKFProp.cc \ PHSimpleVertexFinder.cc \ diff --git a/offline/packages/trackreco/PHActsTrkFitter.cc b/offline/packages/trackreco/PHActsTrkFitter.cc index 3d7c6387a3..2cb0ec8c0c 100644 --- a/offline/packages/trackreco/PHActsTrkFitter.cc +++ b/offline/packages/trackreco/PHActsTrkFitter.cc @@ -618,7 +618,8 @@ void PHActsTrkFitter::loopTracks(Acts::Logging::Level logLevel) m_clusterContainer, m_tGeometry, m_globalPositionWrapper, - this_crossing); + this_crossing, + !m_forceSiOnlyFit); //for alignment tests, if forceSiOnlyFit is true, then we want to turn off the parametrized cluster errors when building source links } // tpc source links @@ -1294,8 +1295,8 @@ void PHActsTrkFitter::updateSvtxTrack( track->set_z(params.position(m_transient_geocontext)(2) / Acts::UnitConstants::cm); auto* seed = track->get_tpc_seed(); - - if(!m_forceSiOnlyFit) + + if(!m_forceSiOnlyFit || !m_siOnlyTpcSeedPt) { track->set_px(params.momentum()(0)); track->set_py(params.momentum()(1)); diff --git a/offline/packages/trackreco/PHActsTrkFitter.h b/offline/packages/trackreco/PHActsTrkFitter.h index 2903a67ebe..51ba9b340c 100644 --- a/offline/packages/trackreco/PHActsTrkFitter.h +++ b/offline/packages/trackreco/PHActsTrkFitter.h @@ -86,6 +86,10 @@ class PHActsTrkFitter : public SubsysReco { m_forceSiOnlyFit = forceSiOnlyFit; } + void forceSiOnlyFitTpcSeedPT() + { + m_siOnlyTpcSeedPt = true; + } /// FOR ALIGNMENT STUDIES ONLY, USE AT OWN RISK. With direct navigation, force a fit with only tpc hits and a full /// matched (si+tpc track seed). This requires a standard track fit to be run first, followed by refit configured with @@ -241,6 +245,7 @@ class PHActsTrkFitter : public SubsysReco /// Acts::DirectedNavigator with a list of sorted silicon+MM surfaces bool m_fitSiliconMMs = false; bool m_forceSiOnlyFit = false; + bool m_siOnlyTpcSeedPt = false; bool m_forceTpcOnlyFit = false; /// requires micromegas present when fitting silicon-MM surfaces diff --git a/offline/packages/trackreco/PHCosmicSeeder.cc b/offline/packages/trackreco/PHCosmicSeeder.cc index 1de07cb203..50e53692d6 100644 --- a/offline/packages/trackreco/PHCosmicSeeder.cc +++ b/offline/packages/trackreco/PHCosmicSeeder.cc @@ -16,6 +16,7 @@ #include #include #include +#include #include #include @@ -85,6 +86,10 @@ int PHCosmicSeeder::InitRun(PHCompositeNode* topNode) //____________________________________________________________________________.. int PHCosmicSeeder::process_event(PHCompositeNode* /*unused*/) { + if (Verbosity() > 1) + { + std::cout<<"processing next event"<getHitSetKeys(m_trackerId)) { @@ -100,16 +105,49 @@ int PHCosmicSeeder::process_event(PHCompositeNode* /*unused*/) clusterPositions.insert(std::make_pair(ckey, global)); } } + if(clusterPositions.size()<3) + return Fun4AllReturnCodes::ABORTEVENT; + if (Verbosity() > 1) + { + std::cout<<"processing next event2"<getHitSetKeys(TrkrDefs::TrkrId::inttId)) + { + auto range = m_clusterContainer->getClusters(hitsetkey); + for (auto citer = range.first; citer != range.second; ++citer) + { + const auto ckey = citer->first; + const auto cluster = citer->second; + if(cluster->getMaxAdc() < m_adcCut){ + continue; + } + const auto global = m_tGeometry->getGlobalPosition(ckey, cluster); + clusterPositions.insert(std::make_pair(ckey, global)); + } + } + } + if (Verbosity() > 1) + { + std::cout<<"processing next event 3: "< 2) { std::cout << "cluster map size is " << clusterPositions.size() << std::endl; } + if(m_trackerId==TrkrDefs::TrkrId::mvtxId && (clusterPositions.size()<4 || clusterPositions.size() > 500)) + { + return Fun4AllReturnCodes::ABORTEVENT; + } auto seeds = makeSeeds(clusterPositions); if (Verbosity() > 1) { std::cout << "Initial seed candidate size is " << seeds.size() << std::endl; } + if (int(seeds.size())>500) + return Fun4AllReturnCodes::ABORTEVENT; std::sort(seeds.begin(), seeds.end(), [](const seed& a, const seed& b) { return a.ckeys.size() > b.ckeys.size(); }); @@ -176,10 +214,11 @@ int PHCosmicSeeder::process_event(PHCompositeNode* /*unused*/) { svtxseed->insert_cluster_key(key); } - // if (m_trackerId == TrkrDefs::TrkrId::mvtxId){ - // svtxseed->circleFitByTaubin(clusterPositions, 0, 3); - // svtxseed->lineFit(clusterPositions, 0, 3); - // } + if (m_trackerId == TrkrDefs::TrkrId::mvtxId){ + TrackSeedHelper::circleFitByTaubin(svtxseed.get(),clusterPositions, 0, 7); + TrackSeedHelper::lineFit(svtxseed.get(),clusterPositions, 0, 7); + svtxseed->set_phi(TrackSeedHelper::get_phi(svtxseed.get(),clusterPositions)); + } m_seedContainer->insert(svtxseed.get()); ++iseed; } @@ -187,6 +226,8 @@ int PHCosmicSeeder::process_event(PHCompositeNode* /*unused*/) { std::cout << "Final n seeds: " << m_seedContainer->size() << std::endl; } + if (m_seedContainer->size()==0) + return Fun4AllReturnCodes::ABORTEVENT; ++m_event; return Fun4AllReturnCodes::EVENT_OK; } @@ -218,7 +259,8 @@ PHCosmicSeeder::SeedVector PHCosmicSeeder::findIntersections(PHCosmicSeeder::See //! so merge and delete for (auto key : seed2.ckeys) { - seed1.ckeys.insert(key); + if(m_trackerId!=TrkrDefs::TrkrId::mvtxId || !seedContainsClusterSensor(seed1,key)) + seed1.ckeys.insert(key); } seedsToDelete.insert(j); } @@ -284,7 +326,10 @@ PHCosmicSeeder::SeedVector PHCosmicSeeder::chainSeeds(PHCosmicSeeder::SeedVector float pdiff_tol = 1.0; if (m_trackerId == TrkrDefs::TrkrId::mvtxId) { - pdiff_tol = 0.25; + if(m_zerofield) + pdiff_tol = 0.25; + else + pdiff_tol = 1.0; } float const pdiff = std::abs((seed1.xyslope - seed2.xyslope) / longestxyslope); float const pdiff2 = std::abs((seed1.xyintercept - seed2.xyintercept) / longestxyint); @@ -292,7 +337,7 @@ PHCosmicSeeder::SeedVector PHCosmicSeeder::chainSeeds(PHCosmicSeeder::SeedVector float const pdiff4 = std::abs((seed1.xzslope - seed2.xzslope) / longestxzslope); if (Verbosity() > 1) { - std::cout << "pdiff1,2,3,4 " << pdiff << ", " << pdiff2 << ", " << pdiff3 << ", " << pdiff4 << std::endl; + std::cout << "tol: "< contains_layer = {false, false, false}; + //std::vector contains_layer = {false, false, false}; + int mvtx_layers = 0; for (auto& key : initialSeeds[i].ckeys) { - contains_layer[TrkrDefs::getLayer(key)] = true; + if(int(TrkrDefs::getLayer(key))<3) + mvtx_layers++; + //std::cout< 1 && m_trackerId == TrkrDefs::TrkrId::mvtxId) - { - PHCosmicSeeder::SeedVector emptyseeds; - return emptyseeds; - } - + //if (returnseeds.size() > 1 && m_trackerId == TrkrDefs::TrkrId::mvtxId) + //{ + // PHCosmicSeeder::SeedVector emptyseeds; + // return emptyseeds; + //} +// return returnseeds; } PHCosmicSeeder::SeedVector PHCosmicSeeder::combineSeeds(PHCosmicSeeder::SeedVector& initialSeeds, @@ -358,26 +406,38 @@ PHCosmicSeeder::SeedVector PHCosmicSeeder::combineSeeds(PHCosmicSeeder::SeedVect recalculateSeedLineParameters(seed1, clusterPositions, false); recalculateSeedLineParameters(seed2, clusterPositions, false); - float slope_tol = 0.1; - float incept_tol = 3.0; + float slope_tol_xy = 0.1; + float incept_tol_xy = 3.0; + float slope_tol_xz = 0.1; + float incept_tol_xz = 3.0; if (m_trackerId == TrkrDefs::TrkrId::mvtxId) { - slope_tol = 0.1; - incept_tol = 0.5; + slope_tol_xy = 0.1; + incept_tol_xy = 0.5; + slope_tol_xz = 0.1; + incept_tol_xz = 0.5; + if(!m_zerofield) + { + slope_tol_xy = 1.0; + incept_tol_xy = 2.0; + slope_tol_xz = 1.0; + incept_tol_xz = 2.0; + } } if (Verbosity() > 4) { std::cout << "xy slope diff: " << seed1.xyslope - seed2.xyslope << ", xy int diff " << seed1.xyintercept - seed2.xyintercept << ", xz slope diff " << seed1.xzslope - seed2.xzslope << ", xz int diff " << seed1.xzintercept - seed2.xzintercept << std::endl; } //! These values are tuned on the cosmic data - if (std::abs(seed1.xyslope - seed2.xyslope) < slope_tol && - std::abs(seed1.xyintercept - seed2.xyintercept) < incept_tol && - std::abs(seed1.xzslope - seed2.xzslope) < slope_tol && - std::abs(seed1.xzintercept - seed2.xzintercept) < incept_tol) + if (std::abs(seed1.xyslope - seed2.xyslope) < slope_tol_xy && + std::abs(seed1.xyintercept - seed2.xyintercept) < incept_tol_xy && + std::abs(seed1.xzslope - seed2.xzslope) < slope_tol_xz && + std::abs(seed1.xzintercept - seed2.xzintercept) < incept_tol_xz) { for (auto& key : seed2.ckeys) { - seed1.ckeys.insert(key); + if(m_trackerId!=TrkrDefs::TrkrId::mvtxId || !seedContainsClusterSensor(seed1,key)) + seed1.ckeys.insert(key); } seedsToDelete.insert(j); } @@ -417,6 +477,10 @@ PHCosmicSeeder::makeSeeds(PHCosmicSeeder::PositionMap& clusterPositions) { continue; } + if(m_trackerId==TrkrDefs::TrkrId::mvtxId && TrkrDefs::getHitSetKeyFromClusKey(key1) == TrkrDefs::getHitSetKeyFromClusKey(key2)) + { + continue; + } // make a cut on clusters to at least be close to each other within a few cm float const dist = (pos2 - pos1).norm(); if (m_trackerId == TrkrDefs::TrkrId::mvtxId && (TrkrDefs::getLayer(key1) == TrkrDefs::getLayer(key2))) @@ -461,14 +525,14 @@ PHCosmicSeeder::makeSeeds(PHCosmicSeeder::PositionMap& clusterPositions) std::cout << "doublet has " << dub.ckeys.size() << " keys " << std::endl; for (auto key : dub.ckeys) { - std::cout << "position is " << clusterPositions.find(key)->second.transpose() << std::endl; + std::cout << "position is " << clusterPositions.find(key)->second.transpose() << "layer is " << int(TrkrDefs::getLayer(key))<< std::endl; } } + auto begin = dub.ckeys.begin(); auto pos1 = clusterPositions.find(*(begin))->second; std::advance(begin, 1); - auto pos2 = clusterPositions.find(*(begin))->second; - + auto pos2 = clusterPositions.find(*(begin))->second; for (auto& [key, pos] : clusterPositions) { //! skip existing keys @@ -482,7 +546,11 @@ PHCosmicSeeder::makeSeeds(PHCosmicSeeder::PositionMap& clusterPositions) float dist12_check = 2.; if (m_trackerId == TrkrDefs::TrkrId::mvtxId) { - dist12_check = 3.5; + dist12_check = 4.0; + if(TrkrDefs::getLayer(key) > 2) + { + dist12_check = 15.0; + } } if (dist1 < dist2) { @@ -504,13 +572,22 @@ PHCosmicSeeder::makeSeeds(PHCosmicSeeder::PositionMap& clusterPositions) float const predz2 = dub.yzslope * pos.y() + dub.yzintercept; if (Verbosity() > 2) { + std::cout << "testing key with layer "<< int(TrkrDefs::getLayer(key))< 0.3 || fabs(predz2 - pos.z()) > 0.3)) + if (m_trackerId == TrkrDefs::TrkrId::mvtxId && TrkrDefs::getLayer(key)<3 && (fabs(predz - pos.z()) > 0.5 || fabs(predz2 - pos.z()) > 0.3)) + { + continue; + } + else if( m_trackerId == TrkrDefs::TrkrId::mvtxId && TrkrDefs::getLayer(key) >= 3 && (fabs(predz - pos.z()) > 1.1 || fabs(predz2 - pos.z()) > 1.1)) { continue; } @@ -541,10 +618,10 @@ PHCosmicSeeder::makeSeeds(PHCosmicSeeder::PositionMap& clusterPositions) } std::cout << "seed xy slope: " << seed_A.xyslope << std::endl; std::cout << std::endl - << "seed pos " << std::endl; + << " hitsetkey and seed pos " << std::endl; for (auto& key : seed_A.ckeys) { - std::cout << clusterPositions.find(key)->second.transpose() << std::endl; + std::cout << int(TrkrDefs::getLayer(key))<<" , "<second.transpose() << std::endl; } std::cout << "Done printing seed info" << std::endl; } @@ -615,6 +692,32 @@ void PHCosmicSeeder::recalculateSeedLineParameters(seed& seed_A, seed_A.xzintercept = avgy - seed_A.xzslope * avgx; } } +bool PHCosmicSeeder::seedContainsClusterSensor(seed& seed_in, TrkrDefs::cluskey key) +{ + for (auto& existing_key : seed_in.ckeys) + { + if(TrkrDefs::getTrkrId(key) == TrkrDefs::TrkrId::mvtxId&&TrkrDefs::getTrkrId(existing_key) == TrkrDefs::TrkrId::mvtxId) + { + if(TrkrDefs::getLayer(existing_key) == TrkrDefs::getLayer(key)&& + MvtxDefs::getStaveId(existing_key) == MvtxDefs::getStaveId(key)&& + MvtxDefs::getChipId(existing_key) == MvtxDefs::getChipId(key)) + { + return true; + } + } + if(TrkrDefs::getTrkrId(key) == TrkrDefs::TrkrId::inttId &&TrkrDefs::getTrkrId(existing_key) == TrkrDefs::TrkrId::inttId) + { + //std::cout<(topNode, "ActsGeometry"); diff --git a/offline/packages/trackreco/PHCosmicSeeder.h b/offline/packages/trackreco/PHCosmicSeeder.h index c9088629b3..1932e5a13f 100644 --- a/offline/packages/trackreco/PHCosmicSeeder.h +++ b/offline/packages/trackreco/PHCosmicSeeder.h @@ -7,6 +7,7 @@ #include #include #include +#include #include #include @@ -43,6 +44,7 @@ class PHCosmicSeeder : public SubsysReco void seedAnalysis() { m_analysis = true; } void trackMapName(const std::string &name) { m_trackMapName = name; } void trackerId(TrkrDefs::TrkrId trackerId) { m_trackerId = trackerId; } + void fieldOn() {m_zerofield=false; } private: int getNodes(PHCompositeNode *topNode); @@ -52,6 +54,7 @@ class PHCosmicSeeder : public SubsysReco SeedVector findIntersections(SeedVector &initialSeeds); SeedVector chainSeeds(SeedVector &initialSeeds, PositionMap &clusterPositions); void recalculateSeedLineParameters(seed &seed, PositionMap &clusters, bool isXY); + bool seedContainsClusterSensor(seed& seed_in, TrkrDefs::cluskey key); float m_xyTolerance = 2.; //! cm // float m_xzTolerance = 2.; //! cm @@ -64,6 +67,7 @@ class PHCosmicSeeder : public SubsysReco TFile *m_outfile = nullptr; TNtuple *m_tup = nullptr; bool m_analysis = false; + bool m_zerofield = true; float m_event = 0; }; diff --git a/offline/packages/trackreco/PHSiliconTpcTrackMatchingDummy.cc b/offline/packages/trackreco/PHSiliconTpcTrackMatchingDummy.cc new file mode 100644 index 0000000000..0d1601894a --- /dev/null +++ b/offline/packages/trackreco/PHSiliconTpcTrackMatchingDummy.cc @@ -0,0 +1,567 @@ +#include "PHSiliconTpcTrackMatchingDummy.h" + +/// Tracking includes +#include +#include +#include +#include +#include +#include +#include // for cluskey, getTrkrId, tpcId + +#include +#include +#include +#include + +#include // for SvtxVertex +#include + +#include // for PHG4Hit +#include // for keytype +#include // for PHG4Particle + +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include // for UINT_MAX +#include // for fabs, sqrt +#include // for operator<<, basic_ostream +#include +#include // for _Rb_tree_const_iterator +#include // for pair + +using namespace std; + +//____________________________________________________________________________.. +PHSiliconTpcTrackMatchingDummy::PHSiliconTpcTrackMatchingDummy(const std::string &name) + : SubsysReco(name) + , PHParameterInterface(name) +{ + InitializeParameters(); +} + +//____________________________________________________________________________.. +PHSiliconTpcTrackMatchingDummy::~PHSiliconTpcTrackMatchingDummy() = default; + +//____________________________________________________________________________.. +int PHSiliconTpcTrackMatchingDummy::InitRun(PHCompositeNode *topNode) +{ + UpdateParametersWithMacro(); + if(_test_windows) + { + _file = new TFile(_file_name.c_str(), "RECREATE"); + _tree = new TNtuple("track_match", "track_match", + "event:sicrossing:siq:siphi:sieta:six:siy:siz:sipx:sipy:sipz:tpcq:tpcphi:tpceta:tpcx:tpcy:tpcz:tpcpx:tpcpy:tpcpz:tpcid:siid"); + } + // put these in the output file + cout << PHWHERE << " Search windows: phi " << _phi_search_win << " eta " + << _eta_search_win << " _pp_mode " << _pp_mode << " _use_intt_crossing " << _use_intt_crossing << endl; + + int ret = GetNodes(topNode); + if (ret != Fun4AllReturnCodes::EVENT_OK) + { + return ret; + } + std::istringstream stringline(m_fieldMap); + stringline >> fieldstrength; + + // initialize the WindowMatchers + window_dx.init_bools("dx", _print_windows || Verbosity() >0); + window_dy.init_bools("dy", _print_windows || Verbosity() >0); + window_dz.init_bools("dz", _print_windows || Verbosity() >0); + window_dphi.init_bools("dphi", _print_windows || Verbosity() >0); + window_deta.init_bools("deta", _print_windows || Verbosity() >0); + + + return ret; +} + +//_____________________________________________________________________ +void PHSiliconTpcTrackMatchingDummy::SetDefaultParameters() +{ + // Data on gasses @20 C and 760 Torr from the following source: + // http://www.slac.stanford.edu/pubs/icfa/summer98/paper3/paper3.pdf + // diffusion and drift velocity for 400kV for NeCF4 50/50 from calculations: + // http://skipper.physics.sunysb.edu/~prakhar/tpc/HTML_Gases/split.html + + return; +} + +std::string PHSiliconTpcTrackMatchingDummy::WindowMatcher::print_fn(const Arr3D& dat) { + std::ostringstream os; + if (dat[1]==0.) { + os << dat[0]; + } else { + os << dat[0] << (dat[1]>0 ? "+" : "") << dat[1] <<"*exp("< fn_exp(posLo, posLo_b0, pt) + && delta < fn_exp(posHi, posHi_b0, pt)); + } + } else { + double pt = (tpc_pt fn_exp(negLo, negLo_b0, pt) + && delta < fn_exp(negHi, negHi_b0, pt)); + } + } +} + +//____________________________________________________________________________.. +int PHSiliconTpcTrackMatchingDummy::process_event(PHCompositeNode * /*unused*/) +{ + if(Verbosity() > 2) + { + std::cout << " Warning: PHSiliconTpcTrackMatching " + << ( _zero_field ? "zero field is ON" : " zero field is OFF") << std::endl; + } + // _track_map contains the TPC seed track stubs + // _track_map_silicon contains the silicon seed track stubs + // _svtx_seed_map contains the combined silicon and tpc track seeds + + // in case these objects are in the input file, we clear the nodes and replace them + _svtx_seed_map->Reset(); + _track_map->Reset(); + + if (Verbosity() > 0) + { + cout << PHWHERE << " TPC track map size " << _track_map->size() << " Silicon track map size " << _track_map_silicon->size() << endl; + } + + if (_track_map_silicon->size() == 0) + { + return Fun4AllReturnCodes::EVENT_OK; + } + + unsigned int si_id = 0; + // loop over the silicon seeds and add the crossing to them + for (unsigned int trackid = 0; trackid != _track_map_silicon->size(); ++trackid) + { + _tracklet_si = _track_map_silicon->get(trackid); + if (!_tracklet_si) + { + continue; + } + auto crossing = _tracklet_si->get_crossing(); + if (Verbosity() > 8) + { + std::cout << " silicon stub: " << trackid << " eta " << _tracklet_si->get_eta() + << " pt " << _tracklet_si->get_pt() << " si z " << TrackSeedHelper::get_z(_tracklet_si) + << " crossing " << crossing << std::endl; + } + + if (Verbosity() > 1) + { + cout << " Si track " << trackid << " crossing " << crossing << endl; + } + auto dummy = std::make_unique(); + dummy->set_qOverR(_tracklet_si->get_qOverR()); + dummy->set_phi(_tracklet_si->get_phi()); + dummy->set_X0(_tracklet_si->get_X0()); + dummy->set_Y0(_tracklet_si->get_Y0()); + dummy->set_Z0(_tracklet_si->get_Z0()); + dummy->set_slope(_tracklet_si->get_slope()); + + auto svtxseed = std::make_unique(); + svtxseed->set_silicon_seed_index(si_id); + svtxseed->set_tpc_seed_index(si_id); + // In pp mode, if a matched track does not have INTT clusters we have to find the crossing geometrically + // Record the geometrically estimated crossing in the track seeds for later use if needed + svtxseed->set_crossing_estimate(crossing); + _track_map->insert(dummy.get()); + _svtx_seed_map->insert(svtxseed.get()); + + si_id++; + if (Verbosity() > 1) + { + std::cout << " combined seed id " << _svtx_seed_map->size() - 1 << " si id " << si_id << " tpc id " << si_id << " crossing estimate " << crossing << std::endl; + } + } + + + if (Verbosity() > 0) + { + std::cout << "final svtx seed map size " << _svtx_seed_map->size() << std::endl; + } + + if (Verbosity() > 1) + { + for (const auto &seed : *_svtx_seed_map) + { + seed->identify(); + std::cout << std::endl; + } + + cout << "PHSiliconTpcTrackMatchingDummy::process_event(PHCompositeNode *topNode) Leaving process_event" << endl; + } + m_event++; + return Fun4AllReturnCodes::EVENT_OK; +} + +double PHSiliconTpcTrackMatchingDummy::getBunchCrossing(unsigned int trid, double z_mismatch) +{ + const double vdrift = _tGeometry->get_drift_velocity(); // cm/ns + const double z_bunch_separation = sphenix_constants::time_between_crossings * vdrift; // cm + + // The sign of z_mismatch will depend on which side of the TPC the tracklet is in + TrackSeed *track = _track_map->get(trid); + + // crossing + double crossings = z_mismatch / z_bunch_separation; + + // Check the TPC side for the first cluster in the track + unsigned int side = 10; + std::set side_set; + for (TrackSeed::ConstClusterKeyIter iter = track->begin_cluster_keys(); + iter != track->end_cluster_keys(); + ++iter) + { + TrkrDefs::cluskey cluster_key = *iter; + unsigned int trkrid = TrkrDefs::getTrkrId(cluster_key); + if (trkrid == TrkrDefs::tpcId) + { + side = TpcDefs::getSide(cluster_key); + side_set.insert(side); + } + } + + if (side == 10) + { + return SHRT_MAX; + } + + if (side_set.size() == 2 && Verbosity() > 1) + { + std::cout << " WARNING: tpc seed " << trid << " changed TPC sides, " + << " final side " << side << std::endl; + } + + // if side = 1 (north, +ve z side), a positive t0 will make the cluster late relative to true z, so it will look like z is less positive + // so a negative z mismatch for side 1 means a positive t0, and positive crossing, so reverse the sign for side 1 + if (side == 1) + { + crossings *= -1.0; + } + + if (Verbosity() > 1) + { + std::cout << " gettrackid " << trid << " side " << side << " z_mismatch " << z_mismatch << " crossings " << crossings << std::endl; + } + + return crossings; +} + +int PHSiliconTpcTrackMatchingDummy::End(PHCompositeNode * /*unused*/) +{ + if(_test_windows) + { + _file->cd(); + _tree->Write(); + _file->Close(); + } + return Fun4AllReturnCodes::EVENT_OK; +} + +int PHSiliconTpcTrackMatchingDummy::GetNodes(PHCompositeNode *topNode) +{ + //--------------------------------- + // Get additional objects off the Node Tree + //--------------------------------- + + _cluster_crossing_map = findNode::getClass(topNode, "TRKR_CLUSTERCROSSINGASSOC"); + if (!_cluster_crossing_map) + { + //cerr << PHWHERE << " ERROR: Can't find TRKR_CLUSTERCROSSINGASSOC " << endl; + // return Fun4AllReturnCodes::ABORTEVENT; + } + + _track_map_silicon = findNode::getClass(topNode, _silicon_track_map_name); + if (!_track_map_silicon) + { + cerr << PHWHERE << " ERROR: Can't find SiliconTrackSeedContainer " << endl; + return Fun4AllReturnCodes::ABORTEVENT; + } + + _track_map = findNode::getClass(topNode, _track_map_name); + if (!_track_map) + { + cerr << PHWHERE << " ERROR: Can't find " << _track_map_name.c_str() << endl; + return Fun4AllReturnCodes::ABORTEVENT; + } + + _svtx_seed_map = findNode::getClass(topNode, "SvtxTrackSeedContainer"); + if (!_svtx_seed_map) + { + std::cout << "Creating node SvtxTrackSeedContainer" << std::endl; + /// Get the DST Node + PHNodeIterator iter(topNode); + PHCompositeNode *dstNode = dynamic_cast(iter.findFirst("PHCompositeNode", "DST")); + + /// Check that it is there + if (!dstNode) + { + std::cerr << "DST Node missing, quitting" << std::endl; + throw std::runtime_error("failed to find DST node in PHActsSourceLinks::createNodes"); + } + + /// Get the tracking subnode + PHNodeIterator dstIter(dstNode); + PHCompositeNode *svtxNode = dynamic_cast(dstIter.findFirst("PHCompositeNode", "SVTX")); + + /// Check that it is there + if (!svtxNode) + { + svtxNode = new PHCompositeNode("SVTX"); + dstNode->addNode(svtxNode); + } + + _svtx_seed_map = new TrackSeedContainer_v1(); + PHIODataNode *node = new PHIODataNode(_svtx_seed_map, "SvtxTrackSeedContainer", "PHObject"); + svtxNode->addNode(node); + } + + _cluster_map = findNode::getClass(topNode, _cluster_map_name); + if (!_cluster_map) + { + std::cout << PHWHERE << " ERROR: Can't find node " <<_cluster_map_name << std::endl; + return Fun4AllReturnCodes::ABORTEVENT; + } + + _tGeometry = findNode::getClass(topNode, "ActsGeometry"); + if (!_tGeometry) + { + std::cout << PHWHERE << "Error, can't find acts tracking geometry" << std::endl; + return Fun4AllReturnCodes::ABORTEVENT; + } + + return Fun4AllReturnCodes::EVENT_OK; +} + +void PHSiliconTpcTrackMatchingDummy::checkZMatches( + std::multimap &tpc_matches, + std::multimap &bad_map) +{ + // for _pp_mode=false, assume zero crossings for all track matches + // for _pp_mode=true, do crossing correction on track position z according to side and vdrift + // z matching criteria follows window_z + // there is a dz threshold cut to avoid window_z blow up at low pT + + float vdrift = _tGeometry->get_drift_velocity(); + + for (auto [tpcid, si_id] : tpc_matches) + { + TrackSeed *tpc_track = _track_map->get(tpcid); + TrackSeed *si_track = _track_map_silicon->get(si_id); + + short int crossing = si_track->get_crossing(); + float tpc_pt, tpc_z, si_z; + int tpc_q; + if (_zero_field) { + auto cluster_list_tpc = getTrackletClusterList(tpc_track); + auto cluster_list_si = getTrackletClusterList(si_track); + + tpc_pt = std::get<3>(TrackFitUtils::zero_field_track_params(_tGeometry, _cluster_map, cluster_list_tpc)); + tpc_z = std::get<4>(TrackFitUtils::zero_field_track_params(_tGeometry, _cluster_map, cluster_list_tpc)).z(); + tpc_q = -100; + + si_z = std::get<4>(TrackFitUtils::zero_field_track_params(_tGeometry, _cluster_map, cluster_list_si)).z(); + } else { + tpc_pt = fabs(1. / _tracklet_tpc->get_qOverR()) * (0.3 / 100.) * fieldstrength; + tpc_z = TrackSeedHelper::get_z(tpc_track); + tpc_q = _tracklet_tpc->get_charge(); + si_z = TrackSeedHelper::get_z(si_track); + } + + // get TPC side from one of the TPC clusters + std::vector temp_clusters = getTrackletClusterList(tpc_track); + if(temp_clusters.size() == 0) { continue; } + unsigned int this_side = TpcDefs::getSide(temp_clusters[0]); + + bool is_posQ = (tpc_q>0.); + + float z_mismatch = tpc_z - si_z; + float tpc_z_corrected = _clusterCrossingCorrection.correctZ(tpc_z, this_side, crossing); + float z_mismatch_corrected = tpc_z_corrected - si_z; + + bool z_match = false; + if (_pp_mode) + { + if (crossing == SHRT_MAX) + { + if (Verbosity() > 2) + { + std::cout << " drop si_track " << si_id << " with eta " << si_track->get_eta() << " and z " << TrackSeedHelper::get_z(si_track) << " because crossing is undefined " << std::endl; + } + continue; + } + + if (window_dz.in_window(is_posQ, tpc_pt, tpc_z_corrected, si_z) && (fabs(z_mismatch_corrected) < _crossing_deltaz_max)) + { + z_match = true; + } + else if (fabs(z_mismatch_corrected) < _crossing_deltaz_min) + { + z_match = true; + } + } + else + { + if (window_dz.in_window(is_posQ, tpc_pt, tpc_z, si_z) && (fabs(z_mismatch) < _crossing_deltaz_max)) + { + z_match = true; + } + else if (fabs(z_mismatch) < _crossing_deltaz_min) + { + z_match = true; + } + } + + if (z_match) + { + if (Verbosity() > 1) + { + std::cout << " Success: crossing " << crossing << " tpcid " << tpcid << " si id " << si_id + << " tpc z " << tpc_z << " si z " << si_z << " z_mismatch " << z_mismatch << "tpc z corrected " << tpc_z_corrected + << " z_mismatch_corrected " << z_mismatch_corrected << " drift velocity " << vdrift << std::endl; + } + } + else + { + if (Verbosity() > 1) + { + std::cout << " FAILURE: crossing " << crossing << " tpcid " << tpcid << " si id " << si_id + << " tpc z " << tpc_z << " si z " << si_z << " z_mismatch " << z_mismatch << "tpc_z_corrected " << tpc_z_corrected + << " z_mismatch_corrected " << z_mismatch_corrected << std::endl; + } + + bad_map.insert(std::make_pair(tpcid, si_id)); + } + } + + // remove bad entries from tpc_matches + for (auto [tpcid, si_id] : bad_map) + { + // Have to iterate over tpc_matches and examine each pair to find the one matching bad_map + // this logic works because we call the equal range on vertex_map for every id_pair + // so we only delete one entry per equal range call + auto ret = tpc_matches.equal_range(tpcid); + for (auto it = ret.first; it != ret.second; ++it) + { + if (it->first == tpcid && it->second == si_id) + { + if (Verbosity() > 1) + { + std::cout << " erasing tpc_matches entry for tpcid " << tpcid << " si_id " << si_id << std::endl; + } + tpc_matches.erase(it); + break; // the iterator is no longer valid + } + } + } + + return; +} + +std::vector PHSiliconTpcTrackMatchingDummy::getTrackletClusterList(TrackSeed* tracklet) +{ + std::vector cluskey_vec; + for (auto clusIter = tracklet->begin_cluster_keys(); + clusIter != tracklet->end_cluster_keys(); + ++clusIter) + { + auto key = *clusIter; + auto cluster = _cluster_map->findCluster(key); + if (!cluster) + { + if(Verbosity() > 1) + { + std::cout << PHWHERE << "Failed to get cluster with key " << key << std::endl; + } + continue; + } + + /// Make a safety check for clusters that couldn't be attached to a surface + auto surf = _tGeometry->maps().getSurface(key, cluster); + if (!surf) + { + continue; + } + + // drop some bad layers in the TPC completely + unsigned int layer = TrkrDefs::getLayer(key); + if (layer == 7 || layer == 22 || layer == 23 || layer == 38 || layer == 39) + { + continue; + } + + cluskey_vec.push_back(key); + } // end loop over clusters for this track + return cluskey_vec; +} diff --git a/offline/packages/trackreco/PHSiliconTpcTrackMatchingDummy.h b/offline/packages/trackreco/PHSiliconTpcTrackMatchingDummy.h new file mode 100644 index 0000000000..3109f482c5 --- /dev/null +++ b/offline/packages/trackreco/PHSiliconTpcTrackMatchingDummy.h @@ -0,0 +1,216 @@ +// Tell emacs that this is a C++ source +// -*- C++ -*-. +#ifndef PHSILICONTPCTRACKMATCHINGDUMMY_H +#define PHSILICONTPCTRACKMATCHINGDUMMY_H + +#include +#include +#include +#include + +#include +#include + +class PHCompositeNode; +class TrackSeedContainer; +class TrackSeed; +class TrkrClusterContainer; +class TF1; +class TrkrClusterCrossingAssoc; +class TFile; +class TNtuple; + +class PHSiliconTpcTrackMatchingDummy : public SubsysReco, public PHParameterInterface +{ + public: + PHSiliconTpcTrackMatchingDummy(const std::string &name = "PHSiliconTpcTrackMatchingDummy"); + + ~PHSiliconTpcTrackMatchingDummy() override; + + void SetDefaultParameters() override; + + void set_phi_search_window(const double win) { _phi_search_win = win; } + void set_eta_search_window(const double win) { _eta_search_win = win; } + void set_x_search_window(const double win) { _x_search_win = win; } + void set_y_search_window(const double win) { _y_search_win = win; } + void set_z_search_window(const double win) { _z_search_win = win; } + void set_crossing_deltaz_max(const double dz) {_crossing_deltaz_max = dz;} + void set_crossing_deltaz_min(const double dz) {_crossing_deltaz_min = dz;} + void set_deltaeta_min(const double deta) {_deltaeta_min = deta;} + + float get_phi_search_window() const { return _phi_search_win; } + float get_eta_search_window() const { return _eta_search_win; } + float get_x_search_window() const { return _x_search_win; } + float get_y_search_window() const { return _y_search_win; } + float get_z_search_window() const { return _z_search_win; } + + // 2024/01/22 update + struct WindowMatcher { + // --- new method, comparing to a+b*exp(c/pT) + // Each Arr3D object contains a,b,c in order + // Up to four curved are needed: + // - for for positive and negative tracks + // - only one curve if |dX|; + std::string print_fn(const Arr3D&); + Arr3D posLo { 100, 0, 0 }; // (above a2,b2,c2), 100 for |dX| + Arr3D posHi { 100, 0, 0 }; // (above a3,b3,c3), 100 for use _use_legacy_windowing + Arr3D negLo { 100, 0, 0 }; // (above a0,b0,c0), 100 for |dX| + Arr3D negHi { 100, 0, 0 }; // (above a1,b1,c1), 100 for treat all tracks pos Q + + // efficiency flags set during PHSiliconTpcTrackMatchingDummy::InitRun() + bool fabs_max_posQ = true; + bool fabs_max_negQ = true; + bool negLo_b0 = true; + bool negHi_b0 = true; + bool posLo_b0 = true; + bool posHi_b0 = true; + double min_pt_posQ = 0.25; // only grow function windows down to 150 MeV + double min_pt_negQ = 0.25; // only grow function windows down to 150 MeV + + WindowMatcher( + const Arr3D& _posLo={100,0,0}, + const Arr3D& _posHi={100,0,0}, + const Arr3D& _negLo={100,0,0}, + const Arr3D& _negHi={100,0,0}, + const double _min_pt_posQ=0.25, + const double _min_pt_negQ=0.25) + : posLo{_posLo}, posHi{_posHi}, negLo{_negLo}, negHi{_negHi}, + min_pt_posQ{_min_pt_posQ}, min_pt_negQ{_min_pt_negQ} {}; + + inline double fn_exp(const Arr3D& arr, const bool& b_is_0, double pT) { + return (b_is_0 ? arr[0] : arr[0]+arr[1]*exp(arr[2]/pT)); + } + + void init_bools(const std::string& which_window="", const bool print=false); + + bool in_window(bool posQ, const double tpc_pt, const double tpc_X, const double si_X); + + // initialize to fn_lo < deltaX < fn_hi for +Q, and fn_lo < deltaX < fn_hi for -Q + + void reset_fns() { + posLo={100,0,0}; + posHi={100,0,0}; + negLo={100,0,0}; + negHi={100,0,0}; + }; + + // same max for |deltaX| for pos and neg Q + void set_QoverpT_maxabs (const Arr3D& _posHi, const double _min_pt=0.25) + { reset_fns(); posHi=_posHi; min_pt_posQ = _min_pt; }; + + // same range for deltaX for pos and neg Q + void set_QoverpT_range (const Arr3D& _posLo, const Arr3D& _posHi, const double _min_pt=0.25) + { reset_fns(); posLo=_posLo; posHi=_posHi; min_pt_posQ=_min_pt; }; + + // max for |deltaX| for pos Q + void set_posQoverpT_maxabs (const Arr3D& _posHi, const double _min_pt=0.25) + { posLo={100.,0.,0.}; posHi=_posHi; min_pt_posQ = _min_pt; }; + + // max for |deltaX| for neg Q + void set_negQoverpT_maxabs (const Arr3D& _negHi, const double _min_pt=0.25) + { posLo={100.,0.,0.}; negHi=_negHi; min_pt_negQ = _min_pt; }; + + // range for deltaX for pos Q + void set_posQoverpT_range (const Arr3D& _posLo, const Arr3D& _posHi, const double _min_pt=0.25) + { posLo=_posLo; posHi=_posHi; min_pt_posQ = _min_pt; }; + + // range for deltaX for neg Q + void set_negQoverpT_range (const Arr3D& _negLo, const Arr3D& _negHi, const double _min_pt=0.25) + { negLo=_negLo; negHi=_negHi; min_pt_negQ = _min_pt; }; + + }; + + // initialize the window matchers with default values + WindowMatcher window_dx { {100.,0.,0.}, {5.3, 0., 0.} }; + WindowMatcher window_dy { {100.,0.,0.}, {5.2, 0., 0.} }; + WindowMatcher window_dz { {100.,0.,0.}, {0., 2.6, 0.38}, {100,0.,0.,}, {0., 1.45, 0.49} }; + WindowMatcher window_dphi { {-0.25, 0., 0.}, {0.05, 0., 0.} }; + WindowMatcher window_deta { {100.,0.,0.}, {0.050, 0.0064, 1.1}, {100,0.,0.,}, {0.045, 0.0031, 1.0} }; + + bool _print_windows = false; + void print_windows(bool print=true) { _print_windows = print; } + + void zeroField(const bool flag) { _zero_field = flag; } + + // void set_use_old_matching(const bool flag) { _use_old_matching = flag; } + + void set_test_windows_printout(const bool test) { _test_windows = test; } + void set_file_name(const std::string &name) { _file_name = name; } + void set_pp_mode(const bool flag) { _pp_mode = flag; } + void set_use_intt_crossing(const bool flag) { _use_intt_crossing = flag; } + void set_cluster_map_name(const std::string &name) + { + _cluster_map_name = name; + } + int InitRun(PHCompositeNode *topNode) override; + + int process_event(PHCompositeNode *) override; + + int End(PHCompositeNode *) override; + + void fieldMap(std::string &fieldmap) { m_fieldMap = fieldmap; } + + void set_silicon_track_map_name(const std::string &map_name) { _silicon_track_map_name = map_name; } + void set_track_map_name(const std::string &map_name) { _track_map_name = map_name; } + void SetIteration(int iter) { _n_iteration = iter; } + + private: + int GetNodes(PHCompositeNode *topNode); + + std::vector getInttCrossings(TrackSeed *si_track); + void checkZMatches(std::multimap &tpc_matches, + std::multimap &bad_map); + short int getCrossingIntt(TrackSeed *_tracklet_si); + double getBunchCrossing(unsigned int trid, double z_mismatch); + + TFile *_file = nullptr; + TNtuple *_tree = nullptr; + + std::string _file_name = "track_match.root"; + + // default values, can be replaced from the macro + double _phi_search_win = 0.01; + double _eta_search_win = 0.004; + double _x_search_win = 0.3; + double _y_search_win = 0.3; + double _z_search_win = 0.4; + + // bool _use_old_matching = false; // normally false + + bool _zero_field = false; // fit straight lines if true + + TrackSeedContainer *_svtx_seed_map{nullptr}; + TrackSeedContainer *_track_map{nullptr}; + TrackSeedContainer *_track_map_silicon{nullptr}; + TrackSeed *_tracklet_tpc{nullptr}; + TrackSeed *_tracklet_si{nullptr}; + TrkrClusterContainer *_cluster_map{nullptr}; + ActsGeometry *_tGeometry{nullptr}; + TrkrClusterCrossingAssoc *_cluster_crossing_map{nullptr}; + int m_event = 0; + std::map _z_mismatch_map; + + TpcClusterZCrossingCorrection _clusterCrossingCorrection; + float _crossing_deltaz_max = 10.0; + float _crossing_deltaz_min = 1.5; + float _deltaeta_min = 0.03; + + // double _collision_rate = 50e3; // input rate for phi correction + // double _reference_collision_rate = 50e3; // reference rate for phi correction + // double _si_vertex_dzmax = 0.25; // mm + double fieldstrength{std::numeric_limits::quiet_NaN()}; + + bool _test_windows = false; + bool _pp_mode = false; + bool _use_intt_crossing = true; // should always be true except for testing + + int _n_iteration = 0; + std::string _track_map_name = "TpcTrackSeedContainer"; + std::string _silicon_track_map_name = "SiliconTrackSeedContainer"; + std::string _cluster_map_name = "TRKR_CLUSTER"; + std::string m_fieldMap = "1.4"; + std::vector getTrackletClusterList(TrackSeed* tracklet); +}; + +#endif // PHSILICONTPCTRACKMATCHINGDUMMY_H