QA/SimulationModules/QAG4SimulationDistortions.cc

0001 #include "QAG4SimulationDistortions.h"
0002
0003 #include <micromegas/MicromegasDefs.h>
0004
0005 #include <trackbase/ActsGeometry.h>
0006 #include <trackbase/TrkrCluster.h>
0007 #include <trackbase/TrkrClusterContainer.h>
0008 #include <trackbase/TrkrDefs.h>
0009 #include <trackbase_historic/SvtxTrack.h>
0010 #include <trackbase_historic/SvtxTrackMap.h>
0011 #include <trackbase_historic/SvtxTrackState.h>
0012 #include <trackbase_historic/TrackAnalysisUtils.h>
0013 #include <trackbase_historic/TrackSeed.h>
0014
0015 #include <g4detectors/PHG4TpcGeom.h>
0016 #include <g4detectors/PHG4TpcGeomContainer.h>
0017
0018 #include <qautils/QAHistManagerDef.h>
0019
0020 #include <fun4all/Fun4AllHistoManager.h>
0021 #include <fun4all/Fun4AllReturnCodes.h>
0022 #include <fun4all/SubsysReco.h>
0023
0024 #include <phool/PHCompositeNode.h>
0025 #include <phool/getClass.h>
0026 #include <phool/phool.h>  // for PHWHERE
0027
0028 #include <Acts/Definitions/Algebra.hpp>
0029
0030 #include <TH1.h>
0031 #include <TH2.h>
0032 #include <TTree.h>
0033
0034 #include <algorithm>
0035 #include <cassert>
0036 #include <cmath>
0037 #include <iostream>
0038 #include <iterator>
0039 #include <limits>
0040 #include <string>
0041 #include <utility>  // for pair
0042 #include <vector>
0043
0044 namespace
0045 {
0046
0047   // square
0048   template <class T>
0049   constexpr T square(const T& x)
0050   {
0051     return x * x;
0052   }
0053
0054   // radius
0055   template <class T>
0056   T get_r(const T& x, const T& y)
0057   {
0058     return std::sqrt(square(x) + square(y));
0059   }
0060
0061   template <class T>
0062   constexpr T deltaPhi(const T& phi)
0063   {
0064     if (phi > M_PI)
0065     {
0066       return phi - 2. * M_PI;
0067     }
0068     if (phi <= -M_PI)
0069     {
0070       return phi + 2. * M_PI;
0071     }
0072
0073     return phi;
0074   }
0075
0076   /// return number of clusters of a given type that belong to a tracks
0077   template <int type>
0078   int count_clusters(const std::vector<TrkrDefs::cluskey>& keys)
0079   {
0080     return std::count_if(keys.begin(), keys.end(),
0081                          [](const TrkrDefs::cluskey& key)
0082                          { return TrkrDefs::getTrkrId(key) == type; });
0083   }
0084 }  // namespace
0085
0086 //____________________________________________________________________________..
0087 QAG4SimulationDistortions::QAG4SimulationDistortions(const std::string& name)
0088   : SubsysReco(name)
0089 {
0090 }
0091
0092 //____________________________________________________________________________..
0093 int QAG4SimulationDistortions::Init(PHCompositeNode* /*unused*/)
0094 {
0095   // reset counters
0096   m_total_tracks = 0;
0097   m_accepted_tracks = 0;
0098
0099   m_total_states = 0;
0100   m_accepted_states = 0;
0101
0102   // histogram manager
0103   auto* hm = QAHistManagerDef::getHistoManager();
0104   assert(hm);
0105
0106   h_beta = new TH2F(std::string(get_histo_prefix() + "betadz").c_str(), ";tan#beta; #Deltaz [cm]", 100, -0.5, 0.5, 100, -0.5, 0.5);
0107   hm->registerHisto(h_beta);
0108
0109   h_alpha = new TH2F(std::string(get_histo_prefix() + "alphardphi").c_str(), ";tan#alpha; r#Delta#phi [cm]", 100, -0.5, 0.5, 100, -0.5, 0.5);
0110   hm->registerHisto(h_alpha);
0111
0112   h_rphiResid = new TH2F(std::string(get_histo_prefix() + "rphiResid").c_str(), ";r [cm]; #Deltar#phi [cm]", 60, 20, 80, 500, -2, 2);
0113   hm->registerHisto(h_rphiResid);
0114
0115   h_zResid = new TH2F(std::string(get_histo_prefix() + "zResid").c_str(), ";z [cm]; #Deltaz [cm]", 200, -100, 100, 1000, -2, 2);
0116   hm->registerHisto(h_zResid);
0117
0118   h_etaResid = new TH2F(std::string(get_histo_prefix() + "etaResid").c_str(), ";#eta;#Delta#eta", 20, -1, 1, 500, -0.2, 0.2);
0119   hm->registerHisto(h_etaResid);
0120
0121   h_etaResidLayer = new TH2F(std::string(get_histo_prefix() + "etaResidLayer").c_str(), ";r [cm]; #Delta#eta", 60, 20, 80, 500, -0.2, 0.2);
0122   hm->registerHisto(h_etaResidLayer);
0123
0124   h_zResidLayer = new TH2F(std::string(get_histo_prefix() + "zResidLayer").c_str(), ";r [cm]; #Deltaz [cm]", 60, 20, 80, 1000, -2, 2);
0125   hm->registerHisto(h_zResidLayer);
0126
0127   h_deltarphi_layer = new TH2F(std::string(get_histo_prefix() + "deltarphi_layer").c_str(), ";layer; r.#Delta#phi_{cluster-track} (cm)", 57, 0, 57, 500, -2, 2);
0128   hm->registerHisto(h_deltarphi_layer);
0129
0130   h_deltaz_layer = new TH2F(std::string(get_histo_prefix() + "deltaz_layer").c_str(), ";layer; #Deltaz_{cluster-track} (cm)", 57, 0, 57, 100, -2, 2);
0131   hm->registerHisto(h_deltaz_layer);
0132
0133   h_statez_pulls = new TH2F(std::string(get_histo_prefix() + "statez_pulls").c_str(), ";layer; #Deltaz_{cluster-track}/#sigma_{z}^{state}", 57, 0, 57, 100, -5, 5);
0134   hm->registerHisto(h_statez_pulls);
0135
0136   h_staterphi_pulls = new TH2F(std::string(get_histo_prefix() + "staterphi_pulls").c_str(), ";layer; #Deltar#phi_{cluster-track}/#sigma_{rphi}^{state}", 57, 0, 57, 100, -5, 5);
0137   hm->registerHisto(h_staterphi_pulls);
0138
0139   h_clusz_pulls = new TH2F(std::string(get_histo_prefix() + "clusz_pulls").c_str(), ";layer; #Deltaz_{cluster-track}/#sigma_{z}^{clus}", 57, 0, 57, 100, -5, 5);
0140   hm->registerHisto(h_clusz_pulls);
0141
0142   h_clusrphi_pulls = new TH2F(std::string(get_histo_prefix() + "clusrphi_pulls").c_str(), ";layer; #Deltar#phi_{cluster-track}/#sigma_{r#phi}^{clus}", 57, 0, 57, 100, -5, 5);
0143   hm->registerHisto(h_clusrphi_pulls);
0144
0145   h_nstates_vs_nclus = new TH2F(std::string(get_histo_prefix() + "nstates_vs_nclus").c_str(), ";Number of states on track; Number of clusters on track", 57, 0, 57, 57, 0, 57);
0146   hm->registerHisto(h_nstates_vs_nclus);
0147
0148   h_tpot_deltarphi = new TH1F(std::string(get_histo_prefix() + "tpot_deltarphi").c_str(), ";TPOT  r.#Delta#phi_{cluster-track} (cm); Number of TPOT clusters", 100, -1, 1);
0149   hm->registerHisto(h_tpot_deltarphi);
0150
0151   h_tpot_deltaz = new TH1F(std::string(get_histo_prefix() + "tpot_deltaz").c_str(), ";TPOT  #Deltaz_{cluster-track} (cm); Number of TPOT clusters", 100, -2, 2);
0152   hm->registerHisto(h_tpot_deltaz);
0153
0154   t_tree = new TTree(std::string(get_histo_prefix() + "residTree").c_str(), "tpc residual info");
0155   t_tree->Branch("tanAlpha", &m_tanAlpha, "tanAlpha/F");
0156   t_tree->Branch("tanBeta", &m_tanBeta, "tanBeta/F");
0157   t_tree->Branch("drphi", &m_drphi, "drphi/F");
0158   t_tree->Branch("dz", &m_dz, "dz/F");
0159   t_tree->Branch("clusR", &m_clusR, "clusR/F");
0160   t_tree->Branch("clusPhi", &m_clusPhi, "clusPhi/F");
0161   t_tree->Branch("clusZ", &m_clusZ, "clusZ/F");
0162   t_tree->Branch("clusEta", &m_clusEta, "clusEta/F");
0163   t_tree->Branch("stateR", &m_stateR, "stateR/F");
0164   t_tree->Branch("statePhi", &m_statePhi, "statePhi/F");
0165   t_tree->Branch("stateZ", &m_stateZ, "stateZ/F");
0166   t_tree->Branch("stateEta", &m_stateEta, "stateEta/F");
0167   t_tree->Branch("stateRPhiErr", &m_stateRPhiErr, "stateRPhiErr/F");
0168   t_tree->Branch("stateZErr", &m_stateZErr, "stateZErr/F");
0169   t_tree->Branch("clusRPhiErr", &m_clusRPhiErr, "clusRPhiErr/F");
0170   t_tree->Branch("clusZErr", &m_clusZErr, "clusZErr/F");
0171   t_tree->Branch("cluskey", &m_cluskey, "cluskey/l");
0172   t_tree->Branch("event", &m_event, "event/I");
0173
0174   t_tree->Branch("layer", &m_layer, "layer/I");
0175   t_tree->Branch("statePt", &m_statePt, "statePt/F");
0176   t_tree->Branch("statePz", &m_statePz, "statePz/F");
0177
0178   t_tree->Branch("trackPt", &m_trackPt, "trackPt/F");
0179   t_tree->Branch("charge", &m_charge, "charge/I");
0180   t_tree->Branch("crossing", &m_crossing, "crossing/I");
0181   t_tree->Branch("trackdEdx", &m_trackdEdx, "trackdEdx/F");
0182
0183   hm->registerHisto(t_tree);
0184
0185   return Fun4AllReturnCodes::EVENT_OK;
0186 }
0187
0188 //____________________________________________________________________________..
0189 int QAG4SimulationDistortions::InitRun(PHCompositeNode* topNode)
0190 {
0191   // track map
0192   m_trackMap = findNode::getClass<SvtxTrackMap>(topNode, m_trackmapname);
0193
0194   // cluster map
0195   m_clusterContainer = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0196
0197   // find tpc geometry
0198   m_tpcGeom = findNode::getClass<PHG4TpcGeomContainer>(topNode, "TPCGEOMCONTAINER");
0199
0200   // load geometry
0201   m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0202
0203   // load distortion corrections
0204   m_globalPositionWrapper.loadNodes(topNode);
0205   if (m_disable_module_edge_corr)
0206   {
0207     m_globalPositionWrapper.set_enable_module_edge_corr(false);
0208   }
0209   if (m_disable_static_corr)
0210   {
0211     m_globalPositionWrapper.set_enable_static_corr(false);
0212   }
0213   if (m_disable_average_corr)
0214   {
0215     m_globalPositionWrapper.set_enable_average_corr(false);
0216   }
0217   if (m_disable_fluctuation_corr)
0218   {
0219     m_globalPositionWrapper.set_enable_fluctuation_corr(false);
0220   }
0221
0222   if (!m_trackMap || !m_clusterContainer || !m_tGeometry)
0223   {
0224     std::cout << PHWHERE << "Necessary distortion container not on node tree. Bailing."
0225               << std::endl;
0226
0227     return Fun4AllReturnCodes::ABORTRUN;
0228   }
0229
0230   return Fun4AllReturnCodes::EVENT_OK;
0231 }
0232
0233 //____________________________________________________________________________..
0234 int QAG4SimulationDistortions::process_event(PHCompositeNode* /*unused*/)
0235 {
0236   Fun4AllHistoManager* hm = QAHistManagerDef::getHistoManager();
0237   assert(hm);
0238
0239   std::cout << "QAG4SimulationDistortions::process_event - tracks: " << m_trackMap->size() << std::endl;
0240   for (const auto& [key, track] : *m_trackMap)
0241   {
0242     // total track counter
0243     ++m_total_tracks;
0244
0245     // get track crossing and check
0246     const auto crossing = track->get_crossing();
0247     if (crossing == std::numeric_limits<short>::max())
0248     {
0249       std::cout << "QAG4SimulationDistortions::process_event - invalid crossing. Track skipped." << std::endl;
0250       continue;
0251     }
0252
0253     // check track quality
0254     if (!checkTrack(track))
0255     {
0256       if (Verbosity() > 1)
0257       {
0258         std::cout << "failed track selection" << std::endl;
0259       }
0260       continue;
0261     }
0262     if (Verbosity() > 1)
0263     {
0264       std::cout << "pass track selection" << std::endl;
0265     }
0266
0267     // get seeeds
0268     auto* tpcSeed = track->get_tpc_seed();
0269     auto* siliconSeed = track->get_silicon_seed();
0270
0271     /// Should have never been added to the map...
0272     if (!tpcSeed || !siliconSeed)
0273     {
0274       continue;
0275     }
0276
0277     if (Verbosity() > 0)
0278     {
0279       std::cout << " tpc seed " << tpcSeed->get_tpc_seed_index()
0280                 << " with si seed " << siliconSeed->get_silicon_seed_index()
0281                 << " crossing " << siliconSeed->get_crossing()
0282                 << std::endl;
0283     }
0284
0285     // accepted track counter
0286     ++m_accepted_tracks;
0287
0288     for (auto iter = track->begin_states(); iter != track->end_states(); ++iter)
0289     {
0290       ++m_total_states;
0291
0292       auto& state = iter->second;
0293       const auto ckey = state->get_cluskey();
0294       const auto trkrId = TrkrDefs::getTrkrId(ckey);
0295
0296       if (trkrId != TrkrDefs::tpcId)
0297       {
0298         continue;
0299       }
0300
0301       ++m_accepted_states;
0302
0303       auto* cluster = m_clusterContainer->findCluster(ckey);
0304
0305       const auto clusGlobPosition = m_globalPositionWrapper.getGlobalPositionDistortionCorrected(ckey, cluster, crossing);
0306
0307       const float clusR = get_r(clusGlobPosition(0), clusGlobPosition(1));
0308       const float clusPhi = std::atan2(clusGlobPosition(1), clusGlobPosition(0));
0309       const float clusZ = clusGlobPosition(2);
0310
0311       // cluster errors
0312       const float clusRPhiErr = cluster->getRPhiError();
0313       const float clusZErr = cluster->getZError();
0314
0315       const Acts::Vector3 stateGlobPosition = Acts::Vector3(state->get_x(),
0316                                                             state->get_y(),
0317                                                             state->get_z());
0318       const Acts::Vector3 stateGlobMom = Acts::Vector3(state->get_px(),
0319                                                        state->get_py(),
0320                                                        state->get_pz());
0321
0322       const float stateRPhiErr = state->get_rphi_error();
0323       const float stateZErr = state->get_z_error();
0324
0325       const float stateR = get_r(stateGlobPosition(0), stateGlobPosition(1));
0326
0327       const auto dr = clusR - stateR;
0328       const auto trackDrDt = (stateGlobPosition(0) * stateGlobMom(0) + stateGlobPosition(1) * stateGlobMom(1)) / stateR;
0329       const auto trackDxDr = stateGlobMom(0) / trackDrDt;
0330       const auto trackDyDr = stateGlobMom(1) / trackDrDt;
0331       const auto trackDzDr = stateGlobMom(2) / trackDrDt;
0332
0333       const auto trackX = stateGlobPosition(0) + dr * trackDxDr;
0334       const auto trackY = stateGlobPosition(1) + dr * trackDyDr;
0335       const auto trackZ = stateGlobPosition(2) + dr * trackDzDr;
0336       const float statePhi = std::atan2(trackY, trackX);
0337       const float stateZ = trackZ;
0338
0339       // Calculate residuals
0340       const float drphi = clusR * deltaPhi(clusPhi - statePhi);
0341       const float dz = clusZ - stateZ;
0342
0343       const auto trackPPhi = -stateGlobMom(0) * std::sin(statePhi) + stateGlobMom(1) * std::cos(statePhi);
0344       const auto trackPR = stateGlobMom(0) * std::cos(statePhi) + stateGlobMom(1) * std::sin(statePhi);
0345       const auto trackPZ = stateGlobMom(2);
0346
0347       const auto trackAlpha = -trackPPhi / trackPR;
0348       const auto trackBeta = -trackPZ / trackPR;
0349       const auto trackEta = std::atanh(stateGlobMom(2) / stateGlobMom.norm());
0350       const auto clusEta = std::atanh(clusZ / clusGlobPosition.norm());
0351
0352       h_alpha->Fill(trackAlpha, drphi);
0353       h_beta->Fill(trackBeta, dz);
0354       h_rphiResid->Fill(clusR, drphi);
0355       h_zResid->Fill(stateZ, dz);
0356       h_etaResid->Fill(trackEta, clusEta - trackEta);
0357       h_zResidLayer->Fill(clusR, dz);
0358       h_etaResidLayer->Fill(clusR, clusEta - trackEta);
0359
0360       const auto layer = TrkrDefs::getLayer(ckey);
0361       h_deltarphi_layer->Fill(layer, drphi);
0362       h_deltaz_layer->Fill(layer, dz);
0363
0364       h_statez_pulls->Fill(layer, dz / stateZErr);
0365       h_staterphi_pulls->Fill(layer, drphi / stateRPhiErr);
0366       h_clusz_pulls->Fill(layer, dz / clusZErr);
0367       h_clusrphi_pulls->Fill(layer, drphi / clusRPhiErr);
0368
0369       m_tanAlpha = trackAlpha;
0370       m_tanBeta = trackBeta;
0371       m_drphi = drphi;
0372       m_dz = dz;
0373       m_clusR = clusR;
0374       m_clusPhi = clusPhi;
0375       m_clusZ = clusZ;
0376       m_statePhi = statePhi;
0377       m_stateZ = stateZ;
0378       m_stateR = stateR;
0379       m_stateRPhiErr = stateRPhiErr;
0380       m_stateZErr = stateZErr;
0381       m_clusRPhiErr = clusRPhiErr;
0382       m_clusZErr = clusZErr;
0383       m_cluskey = ckey;
0384
0385       m_clusEta = clusEta;
0386       m_stateEta = trackEta;
0387       m_layer = layer;
0388       m_statePt = sqrt(stateGlobMom(0) * stateGlobMom(0) + stateGlobMom(1) * stateGlobMom(1));
0389       m_statePz = stateGlobMom(2);
0390       m_trackPt = track->get_pt();
0391       m_charge = track->get_charge();
0392       m_crossing = track->get_crossing();
0393
0394       float layerThicknesses[4] = {0.0, 0.0, 0.0, 0.0};
0395       // These are randomly chosen layer thicknesses for the TPC, to get the
0396       // correct region thicknesses in an easy to pass way to the helper fxn
0397       layerThicknesses[0] = m_tpcGeom->GetLayerCellGeom(7)->get_thickness();
0398       layerThicknesses[1] = m_tpcGeom->GetLayerCellGeom(8)->get_thickness();
0399       layerThicknesses[2] = m_tpcGeom->GetLayerCellGeom(27)->get_thickness();
0400       layerThicknesses[3] = m_tpcGeom->GetLayerCellGeom(50)->get_thickness();
0401
0402       m_trackdEdx = TrackAnalysisUtils::calc_dedx(tpcSeed, m_clusterContainer, m_tGeometry, layerThicknesses);
0403
0404       t_tree->Fill();
0405     }
0406     int nstates = track->size_states();
0407     int nclus = (track->get_silicon_seed()->size_cluster_keys()) + (track->get_tpc_seed()->size_cluster_keys());
0408     h_nstates_vs_nclus->Fill(nstates, nclus);
0409   }
0410
0411   m_event++;
0412
0413   return Fun4AllReturnCodes::EVENT_OK;
0414 }
0415
0416 //___________________________________________________________________________________
0417 int QAG4SimulationDistortions::End(PHCompositeNode* /*topNode*/)
0418 {
0419   // print counters
0420   std::cout
0421       << "QAG4SimulationDistortions::End -"
0422       << " track statistics total: " << m_total_tracks
0423       << " accepted: " << m_accepted_tracks
0424       << " fraction: " << 100. * m_accepted_tracks / m_total_tracks << "%"
0425       << std::endl;
0426
0427   std::cout
0428       << "QAG4SimulationDistortions::End -"
0429       << " state statistics total: " << m_total_states
0430       << " accepted: " << m_accepted_states << " fraction: "
0431       << 100. * m_accepted_states / m_total_states << "%"
0432       << std::endl;
0433
0434   return Fun4AllReturnCodes::EVENT_OK;
0435 }
0436
0437 //_____________________________________________________________________________
0438 bool QAG4SimulationDistortions::checkTrack(SvtxTrack* track)
0439 {
0440   if (track->get_pt() < m_minPT)
0441   {
0442     return false;
0443   }
0444
0445   // ignore tracks with too few mvtx, intt, tpc and micromegas hits
0446   const auto cluster_keys(get_cluster_keys(track));
0447   if (count_clusters<TrkrDefs::mvtxId>(cluster_keys) < 3)
0448   {
0449     return false;
0450   }
0451   if (count_clusters<TrkrDefs::inttId>(cluster_keys) < 2)
0452   {
0453     return false;
0454   }
0455   if (count_clusters<TrkrDefs::tpcId>(cluster_keys) < 35)
0456   {
0457     return false;
0458   }
0459
0460   const auto state_keys(get_state_keys(track));
0461   if (count_clusters<TrkrDefs::mvtxId>(state_keys) < 3)
0462   {
0463     return false;
0464   }
0465   if (count_clusters<TrkrDefs::inttId>(state_keys) < 2)
0466   {
0467     return false;
0468   }
0469
0470   if (m_useMicromegas && checkTPOTResidual(track) == false)
0471   {
0472     return false;
0473   }
0474
0475   return true;
0476 }
0477
0478 bool QAG4SimulationDistortions::checkTPOTResidual(SvtxTrack* track)
0479 {
0480   Fun4AllHistoManager* hm = QAHistManagerDef::getHistoManager();
0481   assert(hm);
0482
0483   bool flag = true;
0484
0485   int nTPOTcluster = 0;
0486   int nTPOTstate = 0;
0487   int TPOTtileID = -1;
0488   for (const auto& cluskey : get_cluster_keys(track))
0489   {
0490     // make sure cluster is from TPOT
0491     const auto detId = TrkrDefs::getTrkrId(cluskey);
0492     if (detId != TrkrDefs::micromegasId)
0493     {
0494       continue;
0495     }
0496     TPOTtileID = MicromegasDefs::getTileId(cluskey);
0497     nTPOTcluster++;
0498
0499     auto* cluster = m_clusterContainer->findCluster(cluskey);
0500
0501     SvtxTrackState* state = nullptr;
0502
0503     // the track states from the Acts fit are fitted to fully corrected clusters, and are on the surface
0504     for (auto state_iter = track->begin_states();
0505          state_iter != track->end_states();
0506          ++state_iter)
0507     {
0508       SvtxTrackState* tstate = state_iter->second;
0509       auto stateckey = tstate->get_cluskey();
0510       if (stateckey == cluskey)
0511       {
0512         state = tstate;
0513         break;
0514       }
0515     }
0516
0517     const auto layer = TrkrDefs::getLayer(cluskey);
0518
0519     if (!state)
0520     {
0521       if (Verbosity() > 1)
0522       {
0523         std::cout << "   no state for cluster " << cluskey << "  in layer " << layer << std::endl;
0524       }
0525       continue;
0526     }
0527     nTPOTstate++;
0528
0529     const auto crossing = track->get_crossing();
0530     assert(crossing != std::numeric_limits<short>::max());
0531
0532     // calculate residuals with respect to cluster
0533     // Get all the relevant information for residual calculation
0534     const auto globClusPos = m_globalPositionWrapper.getGlobalPositionDistortionCorrected(cluskey, cluster, crossing);
0535     const double clusR = get_r(globClusPos(0), globClusPos(1));
0536     const double clusPhi = std::atan2(globClusPos(1), globClusPos(0));
0537     const double clusZ = globClusPos(2);
0538
0539     const double globStateX = state->get_x();
0540     const double globStateY = state->get_y();
0541     const double globStateZ = state->get_z();
0542     const double globStatePx = state->get_px();
0543     const double globStatePy = state->get_py();
0544     const double globStatePz = state->get_pz();
0545
0546     const double trackR = std::sqrt(square(globStateX) + square(globStateY));
0547
0548     const double dr = clusR - trackR;
0549     const double trackDrDt = (globStateX * globStatePx + globStateY * globStatePy) / trackR;
0550     const double trackDxDr = globStatePx / trackDrDt;
0551     const double trackDyDr = globStatePy / trackDrDt;
0552     const double trackDzDr = globStatePz / trackDrDt;
0553
0554     const double trackX = globStateX + dr * trackDxDr;
0555     const double trackY = globStateY + dr * trackDyDr;
0556     const double trackZ = globStateZ + dr * trackDzDr;
0557     const double trackPhi = std::atan2(trackY, trackX);
0558
0559     // Calculate residuals
0560     // need to be calculated in local coordinates in the future
0561     const double drphi = clusR * deltaPhi(clusPhi - trackPhi);
0562     if (std::isnan(drphi))
0563     {
0564       continue;
0565     }
0566
0567     const double dz = clusZ - trackZ;
0568     if (std::isnan(dz))
0569     {
0570       continue;
0571     }
0572
0573     h_tpot_deltarphi->Fill(drphi);
0574     h_tpot_deltaz->Fill(dz);
0575
0576     if (Verbosity() > 3)
0577     {
0578       std::cout << "QAG4SimulationDistortions::checkTPOTResidual -"
0579                 << " drphi: " << drphi
0580                 << " dz: " << dz
0581                 << std::endl;
0582     }
0583
0584     // check rphi residual for layer 55
0585     if (layer == 55 && std::fabs(drphi) > 0.1)
0586     {
0587       flag = false;
0588       break;
0589     }
0590
0591     // check z residual for layer 56
0592     if (layer == 56 && std::fabs(dz) > 1)
0593     {
0594       flag = false;
0595       break;
0596     }
0597   }
0598
0599   if (flag)
0600   {
0601     // SCOZ has a half dead tile
0602     // only require one TPOT cluster/state from SCOP
0603     if (TPOTtileID == 0)
0604     {
0605       if (nTPOTcluster < 1 || nTPOTstate < 1)
0606       {
0607         flag = false;
0608       }
0609     }
0610     else if (TPOTtileID > 0)
0611     {
0612       if (nTPOTcluster < 2 || nTPOTstate < 2)
0613       {
0614         flag = false;
0615       }
0616     }
0617     else if (TPOTtileID < 0)
0618     {
0619       flag = false;
0620     }
0621   }
0622
0623   return flag;
0624 }
0625
0626 std::vector<TrkrDefs::cluskey> QAG4SimulationDistortions::get_cluster_keys(SvtxTrack* track)
0627 {
0628   std::vector<TrkrDefs::cluskey> out;
0629   for (const auto& seed : {track->get_silicon_seed(), track->get_tpc_seed()})
0630   {
0631     if (seed)
0632     {
0633       std::copy(seed->begin_cluster_keys(), seed->end_cluster_keys(), std::back_inserter(out));
0634     }
0635   }
0636   return out;
0637 }
0638
0639 std::vector<TrkrDefs::cluskey> QAG4SimulationDistortions::get_state_keys(SvtxTrack* track)
0640 {
0641   std::vector<TrkrDefs::cluskey> out;
0642   for (auto state_iter = track->begin_states();
0643        state_iter != track->end_states();
0644        ++state_iter)
0645   {
0646     SvtxTrackState* tstate = state_iter->second;
0647     auto stateckey = tstate->get_cluskey();
0648     out.push_back(stateckey);
0649   }
0650   return out;
0651 }