packages/TrackingDiagnostics/helixResiduals.cc

0001 #include "helixResiduals.h"
0002
0003 #include <fun4all/Fun4AllReturnCodes.h>
0004 #include <phool/PHCompositeNode.h>
0005 #include <phool/getClass.h>
0006 #include <trackbase/TrkrDefs.h>
0007
0008 #include <TFile.h>
0009 #include <TNtuple.h>
0010
0011 helixResiduals::helixResiduals(const std::string &name)
0012   : SubsysReco(name)
0013 {
0014   _fitter = new HelicalFitter;
0015 }
0016
0017 int helixResiduals::InitRun(PHCompositeNode *topNode)
0018 {
0019   _fitter->InitRun(topNode);
0020   const char *cfilepath = filepath.c_str();
0021   fout = new TFile(cfilepath, "recreate");
0022   ntp_residuals = new TNtuple("ntp_residuals", "Seed Residuals", "seed_id:layer:dphi:dz:x:y:z:pt:px:py:pz:crossing:isSilicon:isTpc");
0023
0024   getNodes(topNode);
0025
0026   return 0;
0027 }
0028
0029 int helixResiduals::End(PHCompositeNode * /**topNode*/)
0030 {
0031   fout->cd();
0032   ntp_residuals->Write();
0033   fout->Close();
0034
0035   delete _fitter;
0036
0037   return 0;
0038 }
0039
0040 int helixResiduals::getNodes(PHCompositeNode *topNode)
0041 {
0042   _tracks = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0043   if (!_tracks)
0044   {
0045     std::cerr << PHWHERE << "No SvtxTrackMap on node tree, exiting." << std::endl;
0046     return Fun4AllReturnCodes::ABORTEVENT;
0047   }
0048   _tpc_seeds = findNode::getClass<TrackSeedContainer>(topNode, "TpcTrackSeedContainer");
0049   if (!_tpc_seeds)
0050   {
0051     std::cerr << PHWHERE << " ERROR: Can't find TpcTrackSeedContainer " << std::endl;
0052     return Fun4AllReturnCodes::ABORTEVENT;
0053   }
0054   _si_seeds = findNode::getClass<TrackSeedContainer>(topNode, "SiliconTrackSeedContainer");
0055   if (!_si_seeds)
0056   {
0057     std::cerr << PHWHERE << " ERROR: Can't find SiliconTrackSeedContainer" << std::endl;
0058     return Fun4AllReturnCodes::ABORTEVENT;
0059   }
0060   _clusters = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0061   if (!_clusters)
0062   {
0063     std::cerr << PHWHERE << " ERROR: Can't find TRKR_CLUSTER" << std::endl;
0064     return Fun4AllReturnCodes::ABORTEVENT;
0065   }
0066   tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0067   if (!tGeometry)
0068   {
0069     std::cerr << PHWHERE << "No acts tracking geometry, can't proceed" << std::endl;
0070     return Fun4AllReturnCodes::ABORTEVENT;
0071   }
0072   return Fun4AllReturnCodes::EVENT_OK;
0073 }
0074
0075 int helixResiduals::process_event(PHCompositeNode * /*topNode*/)
0076 {
0077   for (auto tpcseed_iter = _tpc_seeds->begin(); tpcseed_iter != _tpc_seeds->end(); ++tpcseed_iter)
0078   {
0079     int id = _tpc_seeds->index(tpcseed_iter);
0080     TrackSeed *tpcseed = _tpc_seeds->get(id);
0081     if (!tpcseed)
0082     {
0083       continue;
0084     }
0085     std::cout << "processing tpc seed " << id << std::endl;
0086     fill_residuals(tpcseed, id, true);
0087   }
0088   for (auto siseed_iter = _si_seeds->begin(); siseed_iter != _si_seeds->end(); ++siseed_iter)
0089   {
0090     int id = _si_seeds->index(siseed_iter);
0091     TrackSeed *siseed = _si_seeds->get(id);
0092     if (!siseed)
0093     {
0094       continue;
0095     }
0096     std::cout << "processing si seed " << id << std::endl;
0097     fill_residuals(siseed, id, false);
0098   }
0099   std::cout << "done" << std::endl;
0100   return Fun4AllReturnCodes::EVENT_OK;
0101 }
0102
0103 void helixResiduals::fill_residuals(TrackSeed *seed, int seed_id, bool isTpc)
0104 {
0105   if (seed->size_cluster_keys() == 0)
0106   {
0107     return;
0108   }
0109
0110   std::vector<Acts::Vector3> clusterPositions;
0111   std::vector<TrkrDefs::cluskey> clusterKeys;
0112   _fitter->getTrackletClusters(seed, clusterPositions, clusterKeys);
0113   std::vector<float> fitparams = _fitter->fitClusters(clusterPositions, clusterKeys);
0114
0115   const float pt = seed->get_pt();
0116   const float px = seed->get_px();
0117   const float py = seed->get_py();
0118   const float pz = seed->get_pz();
0119   const unsigned int crossing = seed->get_crossing();
0120
0121   for (size_t i = 0; i < clusterPositions.size(); i++)
0122   {
0123     unsigned int layer = TrkrDefs::getLayer(clusterKeys[i]);
0124     Acts::Vector3 position = clusterPositions[i];
0125     Acts::Vector3 pca = _fitter->get_helix_pca(fitparams, position);
0126     float cluster_phi = atan2(position(1), position(0));
0127     float pca_phi = atan2(pca(1), pca(0));
0128     float dphi = cluster_phi - pca_phi;
0129     if (dphi > M_PI)
0130     {
0131       dphi = 2 * M_PI - dphi;
0132     }
0133     if (dphi < -M_PI)
0134     {
0135       dphi = 2 * M_PI + dphi;
0136     }
0137     float dz = position(2) - pca(2);
0138
0139     ntp_residuals->Fill(seed_id, layer, dphi, dz, position(0), position(1), position(2), pt, px, py, pz, crossing, !isTpc, isTpc);
0140   }
0141 }