packages/TrackingDiagnostics/TrackResiduals.cc

0001
0002 #include "TrackResiduals.h"
0003
0004 #include <trackbase/ClusterErrorPara.h>
0005 #include <trackbase/InttDefs.h>
0006 #include <trackbase/MvtxDefs.h>
0007 #include <trackbase/TpcDefs.h>
0008 #include <trackbase/TrackFitUtils.h>
0009 #include <trackbase/TrkrCluster.h>
0010 #include <trackbase/TrkrClusterContainer.h>
0011 #include <trackbase/TrkrHit.h>
0012 #include <trackbase/TrkrHitSet.h>
0013 #include <trackbase/TrkrHitSetContainer.h>
0014
0015 #include <g4detectors/PHG4CylinderGeomContainer.h>
0016 #include <g4detectors/PHG4TpcCylinderGeom.h>
0017 #include <g4detectors/PHG4TpcCylinderGeomContainer.h>
0018
0019 #include <globalvertex/MbdVertex.h>
0020 #include <globalvertex/MbdVertexMap.h>
0021
0022 #include <intt/CylinderGeomIntt.h>
0023 #include <intt/CylinderGeomInttHelper.h>
0024
0025 #include <mbd/MbdPmtContainer.h>
0026 #include <mbd/MbdPmtHit.h>
0027
0028 #include <micromegas/CylinderGeomMicromegas.h>
0029 #include <micromegas/MicromegasDefs.h>
0030
0031 #include <mvtx/CylinderGeom_Mvtx.h>
0032 #include <mvtx/CylinderGeom_MvtxHelper.h>
0033
0034 #include <trackbase_historic/ActsTransformations.h>
0035 #include <trackbase_historic/SvtxAlignmentState.h>
0036 #include <trackbase_historic/SvtxAlignmentStateMap.h>
0037 #include <trackbase_historic/SvtxTrack.h>
0038 #include <trackbase_historic/SvtxTrackMap.h>
0039 #include <trackbase_historic/TrackAnalysisUtils.h>
0040 #include <trackbase_historic/TrackSeed.h>
0041 #include <trackbase_historic/TrackSeedContainer.h>
0042 #include <trackbase_historic/TrackSeedHelper.h>
0043
0044 #include <tpc/LaserEventInfo.h>
0045
0046 #include <ffarawobjects/Gl1Packet.h>
0047 #include <ffarawobjects/Gl1RawHit.h>
0048
0049 #include <globalvertex/GlobalVertex.h>
0050 #include <globalvertex/GlobalVertexMap.h>
0051 #include <globalvertex/SvtxVertex.h>
0052 #include <globalvertex/SvtxVertexMap.h>
0053
0054 #include <fun4all/Fun4AllReturnCodes.h>
0055 #include <fun4all/Fun4AllServer.h>
0056
0057 #include <phool/PHCompositeNode.h>
0058 #include <phool/PHNodeIterator.h>
0059 #include <phool/getClass.h>
0060
0061 #include <TSystem.h>
0062
0063 #include <cmath>
0064 #include <algorithm>
0065 #include <limits>
0066
0067 namespace
0068 {
0069   template <class T>
0070   inline T square(const T& t)
0071   {
0072     return t * t;
0073   }
0074   template <class T>
0075   inline T r(const T& x, const T& y)
0076   {
0077     return std::sqrt(square(x) + square(y));
0078   }
0079
0080   std::vector<TrkrDefs::cluskey> get_cluster_keys(SvtxTrack* track)
0081   {
0082     std::vector<TrkrDefs::cluskey> out;
0083
0084     if(!track)
0085       {
0086     return out;
0087       }
0088
0089     for (const auto& seed : {track->get_silicon_seed(), track->get_tpc_seed()})
0090     {
0091       if (seed)
0092       {
0093         std::copy(seed->begin_cluster_keys(), seed->end_cluster_keys(), std::back_inserter(out));
0094       }
0095     }
0096
0097     return out;
0098   }
0099 }  // namespace
0100
0101 //____________________________________________________________________________..
0102 TrackResiduals::TrackResiduals(const std::string& name)
0103   : SubsysReco(name)
0104 {
0105 }
0106
0107 //____________________________________________________________________________..
0108 int TrackResiduals::InitRun(PHCompositeNode* topNode)
0109 {
0110   m_outfile = new TFile(m_outfileName.c_str(), "RECREATE");
0111   createBranches();
0112
0113   // global position wrapper
0114   m_globalPositionWrapper.loadNodes(topNode);
0115   m_globalPositionWrapper.set_suppressCrossing(m_convertSeeds);
0116   // clusterMover needs the correct radii of the TPC layers
0117   auto *tpccellgeo = findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
0118   m_clusterMover.initialize_geometry(tpccellgeo);
0119   m_clusterMover.set_verbosity(0);
0120
0121   auto *se = Fun4AllServer::instance();
0122   m_runnumber = se->RunNumber();
0123
0124   return Fun4AllReturnCodes::EVENT_OK;
0125 }
0126 void TrackResiduals::clearClusterStateVectors()
0127 {
0128   m_cluskeys.clear();
0129   m_clusphisize.clear();
0130   m_cluszsize.clear();
0131   m_idealsurfcenterx.clear();
0132   m_idealsurfcentery.clear();
0133   m_idealsurfcenterz.clear();
0134   m_idealsurfnormx.clear();
0135   m_clsector.clear();
0136   m_clside.clear();
0137   m_idealsurfnormy.clear();
0138   m_idealsurfnormz.clear();
0139   m_missurfcenterx.clear();
0140   m_missurfcentery.clear();
0141   m_missurfcenterz.clear();
0142   m_missurfnormx.clear();
0143   m_missurfnormy.clear();
0144   m_missurfnormz.clear();
0145   m_clusgxideal.clear();
0146   m_clusgyideal.clear();
0147   m_clusgzideal.clear();
0148   m_missurfalpha.clear();
0149   m_missurfbeta.clear();
0150   m_missurfgamma.clear();
0151   m_idealsurfalpha.clear();
0152   m_idealsurfbeta.clear();
0153   m_idealsurfgamma.clear();
0154
0155   m_statelxglobderivdx.clear();
0156   m_statelxglobderivdy.clear();
0157   m_statelxglobderivdz.clear();
0158   m_statelxglobderivdalpha.clear();
0159   m_statelxglobderivdbeta.clear();
0160   m_statelxglobderivdgamma.clear();
0161
0162   m_statelxlocderivd0.clear();
0163   m_statelxlocderivz0.clear();
0164   m_statelxlocderivphi.clear();
0165   m_statelxlocderivtheta.clear();
0166   m_statelxlocderivqop.clear();
0167
0168   m_statelzglobderivdx.clear();
0169   m_statelzglobderivdy.clear();
0170   m_statelzglobderivdz.clear();
0171   m_statelzglobderivdalpha.clear();
0172   m_statelzglobderivdbeta.clear();
0173   m_statelzglobderivdgamma.clear();
0174
0175   m_statelzlocderivd0.clear();
0176   m_statelzlocderivz0.clear();
0177   m_statelzlocderivphi.clear();
0178   m_statelzlocderivtheta.clear();
0179   m_statelzlocderivqop.clear();
0180
0181   m_clusedge.clear();
0182   m_clusoverlap.clear();
0183   m_cluslx.clear();
0184   m_cluslz.clear();
0185   m_cluselx.clear();
0186   m_cluselz.clear();
0187   m_clusgr.clear();
0188   m_clusgx.clear();
0189   m_clusgy.clear();
0190   m_clusgz.clear();
0191   m_clusgxunmoved.clear();
0192   m_clusgyunmoved.clear();
0193   m_clusgzunmoved.clear();
0194   m_clusAdc.clear();
0195   m_clusMaxAdc.clear();
0196   m_cluslayer.clear();
0197
0198   m_statelx.clear();
0199   m_statelz.clear();
0200   m_stateelx.clear();
0201   m_stateelz.clear();
0202   m_stategx.clear();
0203   m_stategy.clear();
0204   m_stategz.clear();
0205   m_statepx.clear();
0206   m_statepy.clear();
0207   m_statepz.clear();
0208   m_statepl.clear();
0209 }
0210 //____________________________________________________________________________..
0211 int TrackResiduals::process_event(PHCompositeNode* topNode)
0212 {
0213   auto *trackmap = findNode::getClass<SvtxTrackMap>(topNode, m_trackMapName);
0214   auto *clustermap = findNode::getClass<TrkrClusterContainer>(topNode, m_clusterContainerName);
0215   auto *geometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0216   auto *hitmap = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0217   auto *tpcGeom =
0218       findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
0219   auto *mvtxGeom = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MVTX");
0220   auto *inttGeom = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_INTT");
0221   auto *mmGeom = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MICROMEGAS_FULL");
0222
0223   if (!mmGeom)
0224   {
0225     mmGeom = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MICROMEGAS");
0226   }
0227   if (!trackmap or !clustermap or !geometry or (!hitmap && m_doHits))
0228   {
0229     std::cout << "Missing node, can't continue" << std::endl;
0230     return Fun4AllReturnCodes::ABORTEVENT;
0231   }
0232   auto *gl1 = findNode::getClass<Gl1RawHit>(topNode, "GL1RAWHIT");
0233   if (gl1)
0234   {
0235     m_bco = gl1->get_bco();
0236     auto lbshift = m_bco << 24U;
0237     m_bcotr = lbshift >> 24U;
0238   }
0239   else
0240   {
0241     Gl1Packet* gl1PacketInfo = findNode::getClass<Gl1Packet>(topNode, "GL1RAWHIT");
0242     if (!gl1PacketInfo)
0243     {
0244       m_bco = std::numeric_limits<uint64_t>::quiet_NaN();
0245       m_bcotr = std::numeric_limits<uint64_t>::quiet_NaN();
0246     }
0247     m_firedTriggers.clear();
0248
0249     if (gl1PacketInfo)
0250     {
0251       m_gl1BunchCrossing = gl1PacketInfo->getBunchNumber();
0252       uint64_t triggervec = gl1PacketInfo->getScaledVector();
0253       m_bco = gl1PacketInfo->getBCO();
0254       auto lbshift = m_bco << 24U;
0255       m_bcotr = lbshift >> 24U;
0256       for (int i = 0; i < 64; i++)
0257       {
0258         bool trig_decision = ((triggervec & 0x1U) == 0x1U);
0259         if (trig_decision)
0260         {
0261           m_firedTriggers.push_back(i);
0262         }
0263         triggervec = (triggervec >> 1U) & 0xffffffffU;
0264       }
0265     }
0266   }
0267
0268   m_mbdvtxz = std::numeric_limits<float>::quiet_NaN();
0269
0270   MbdVertexMap *mbdvertexmap = findNode::getClass<MbdVertexMap>(topNode, "MbdVertexMap");
0271   if(mbdvertexmap)
0272   {
0273     for (auto & it : *mbdvertexmap)
0274     {
0275       MbdVertex* mbdvertex = it.second;
0276       if (mbdvertex)
0277       {
0278         m_mbdvtxz = mbdvertex->get_z();
0279       }
0280     }
0281   }
0282   MbdPmtContainer* bbcpmts = findNode::getClass<MbdPmtContainer>(topNode, "MbdPmtContainer");
0283   m_totalmbd = 0;
0284   if (bbcpmts)
0285   {
0286     int nPMTs = bbcpmts->get_npmt();
0287     for (int i = 0; i < nPMTs; i++)
0288     {
0289       m_totalmbd += bbcpmts->get_pmt(i)->get_q();
0290     }
0291   }
0292   else
0293   {
0294     if (Verbosity() > 0)
0295     {
0296       std::cout << "TrackResiduals::process_event: Could not find MbdPmtContainer," << std::endl;
0297     }
0298   }
0299
0300   m_ntpcclus = 0;
0301   if (Verbosity() > 1)
0302   {
0303     std::cout << "Track map size is " << trackmap->size() << std::endl;
0304   }
0305
0306   if (m_doHits)
0307   {
0308     fillHitTree(hitmap, geometry, tpcGeom, mvtxGeom, inttGeom, mmGeom);
0309   }
0310
0311   if (m_doClusters)
0312   {
0313     fillClusterTree(clustermap, geometry);
0314   }
0315
0316   if (m_convertSeeds)
0317   {
0318     fillResidualTreeSeeds(topNode);
0319   }
0320   else
0321   {
0322     fillResidualTreeKF(topNode);
0323   }
0324
0325   if (m_doVertex)
0326   {
0327     fillVertexTree(topNode);
0328   }
0329   if (m_doEventTree)
0330   {
0331     fillEventTree(topNode);
0332   }
0333   m_event++;
0334   clearClusterStateVectors();
0335   return Fun4AllReturnCodes::EVENT_OK;
0336 }
0337
0338 float TrackResiduals::calc_dedx(TrackSeed* tpcseed, TrkrClusterContainer* clustermap, PHG4TpcCylinderGeomContainer* tpcGeom)
0339 {
0340   std::vector<TrkrDefs::cluskey> clusterKeys;
0341   clusterKeys.insert(clusterKeys.end(), tpcseed->begin_cluster_keys(),
0342                      tpcseed->end_cluster_keys());
0343
0344   std::vector<float> dedxlist;
0345   for (unsigned long cluster_key : clusterKeys)
0346   {
0347     auto detid = TrkrDefs::getTrkrId(cluster_key);
0348     if (detid != TrkrDefs::TrkrId::tpcId)
0349     {
0350       continue;  // the micromegas clusters are added to the TPC seeds
0351     }
0352     unsigned int layer_local = TrkrDefs::getLayer(cluster_key);
0353     TrkrCluster* cluster = clustermap->findCluster(cluster_key);
0354     float adc = cluster->getAdc();
0355     PHG4TpcCylinderGeom* GeoLayer_local = tpcGeom->GetLayerCellGeom(layer_local);
0356     float thick = GeoLayer_local->get_thickness();
0357     float r = GeoLayer_local->get_radius();
0358     float alpha = (r * r) / (2 * r * TMath::Abs(1.0 / tpcseed->get_qOverR()));
0359     float beta = std::atan(tpcseed->get_slope());
0360     float alphacorr = std::cos(alpha);
0361     if (alphacorr < 0 || alphacorr > 4)
0362     {
0363       alphacorr = 4;
0364     }
0365     float betacorr = std::cos(beta);
0366     if (betacorr < 0 || betacorr > 4)
0367     {
0368       betacorr = 4;
0369     }
0370     adc /= thick;
0371     adc *= alphacorr;
0372     adc *= betacorr;
0373     dedxlist.push_back(adc);
0374     sort(dedxlist.begin(), dedxlist.end());
0375   }
0376   int trunc_min = 0;
0377   if (dedxlist.empty())
0378   {
0379     return std::numeric_limits<float>::quiet_NaN();
0380   }
0381   int trunc_max = (int) dedxlist.size() * 0.7;
0382   float sumdedx = 0;
0383   int ndedx = 0;
0384   for (int j = trunc_min; j <= trunc_max; j++)
0385   {
0386     sumdedx += dedxlist.at(j);
0387     ndedx++;
0388   }
0389   sumdedx /= ndedx;
0390   return sumdedx;
0391 }
0392
0393 void TrackResiduals::fillFailedSeedTree(PHCompositeNode* topNode, std::set<unsigned int>& tpc_seed_ids)
0394 {
0395   auto *tpcseedmap = findNode::getClass<TrackSeedContainer>(topNode, "TpcTrackSeedContainer");
0396   auto *trackmap = findNode::getClass<SvtxTrackMap>(topNode, m_trackMapName);
0397   auto *clustermap = findNode::getClass<TrkrClusterContainer>(topNode, m_clusterContainerName);
0398   auto *geometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0399   auto *silseedmap = findNode::getClass<TrackSeedContainer>(topNode, "SiliconTrackSeedContainer");
0400   auto *svtxseedmap = findNode::getClass<TrackSeedContainer>(topNode, "SvtxTrackSeedContainer");
0401   auto *tpcGeo = findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
0402
0403   if (!tpcseedmap or !trackmap or !clustermap or !silseedmap or !svtxseedmap or !geometry)
0404   {
0405     std::cout << "Missing node, can't continue" << std::endl;
0406     return;
0407   }
0408
0409   for (const auto& seed : *svtxseedmap)
0410   {
0411     if (!seed)
0412     {
0413       continue;
0414     }
0415     m_trackid = svtxseedmap->find(seed);
0416     auto tpcseedindex = seed->get_tpc_seed_index();
0417     if (tpc_seed_ids.find(tpcseedindex) != tpc_seed_ids.end())
0418     {
0419       continue;
0420     }
0421     auto siliconseedindex = seed->get_silicon_seed_index();
0422     auto *tpcseed = tpcseedmap->get(tpcseedindex);
0423     auto *silseed = silseedmap->get(siliconseedindex);
0424
0425     int crossing = SHRT_MAX;
0426     if (silseed)
0427     {
0428       const auto si_pos = TrackSeedHelper::get_xyz(silseed);
0429       m_silseedx = si_pos.x();
0430       m_silseedy = si_pos.y();
0431       m_silseedz = si_pos.z();
0432       crossing = silseed->get_crossing();
0433     }
0434     else
0435     {
0436       const auto tpc_pos = TrackSeedHelper::get_xyz(tpcseed);
0437       m_tpcseedx = tpc_pos.x();
0438       m_tpcseedy = tpc_pos.y();
0439       m_tpcseedz = tpc_pos.z();
0440     }
0441
0442     if (m_zeroField)
0443     {
0444       float pt = fabs(1. / tpcseed->get_qOverR()) * (0.3 / 100) * 0.01;
0445       float phi = tpcseed->get_phi();
0446       m_tpcseedpx = pt * std::cos(phi);
0447       m_tpcseedpy = pt * std::sin(phi);
0448       m_tpcseedpz = pt * std::cosh(tpcseed->get_eta()) * std::cos(tpcseed->get_theta());
0449     }
0450     else
0451     {
0452       m_tpcseedpx = tpcseed->get_px();
0453       m_tpcseedpy = tpcseed->get_py();
0454       m_tpcseedpz = tpcseed->get_pz();
0455     }
0456     m_tpcseedcharge = tpcseed->get_qOverR() > 0 ? 1 : -1;
0457     m_dedx = calc_dedx(tpcseed, clustermap, tpcGeo);
0458     m_nmaps = 0;
0459     m_nmapsstate = 0;
0460     m_nintt = 0;
0461     m_ninttstate = 0;
0462     m_ntpc = 0;
0463     m_ntpcstate = 0;
0464     m_nmms = 0;
0465     m_nmmsstate = 0;
0466     clearClusterStateVectors();
0467     for (auto *tseed : {silseed, tpcseed})
0468     {
0469       if (!tseed)
0470       {
0471         continue;
0472       }
0473       for (auto it = tseed->begin_cluster_keys(); it != tseed->end_cluster_keys(); ++it)
0474       {
0475         auto ckey = *it;
0476         auto *cluster = clustermap->findCluster(ckey);
0477         const Acts::Vector3 global = m_globalPositionWrapper.getGlobalPositionDistortionCorrected(ckey, cluster, crossing);
0478         const auto local = geometry->getLocalCoords(ckey, cluster);
0479         m_cluslx.push_back(local.x());
0480         m_cluslz.push_back(local.y());
0481         m_clusgx.push_back(global.x());
0482         m_clusgy.push_back(global.y());
0483         m_clusgz.push_back(global.z());
0484         float cr = r(global.x(), global.y());
0485         if (global.y() < 0)
0486         {
0487           cr = -cr;
0488         }
0489         m_clusgr.push_back(cr);
0490         auto detid = TrkrDefs::getTrkrId(ckey);
0491         if (detid == TrkrDefs::TrkrId::mvtxId)
0492         {
0493           m_nmaps++;
0494         }
0495         else if (detid == TrkrDefs::TrkrId::inttId)
0496         {
0497           m_nintt++;
0498         }
0499         else if (detid == TrkrDefs::TrkrId::tpcId)
0500         {
0501           m_ntpc++;
0502         }
0503         else if (detid == TrkrDefs::TrkrId::micromegasId)
0504         {
0505           m_nmms++;
0506         }
0507       }
0508     }
0509     m_failedfits->Fill();
0510   }
0511 }
0512 void TrackResiduals::fillVertexTree(PHCompositeNode* topNode)
0513 {
0514   auto *svtxvertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
0515   auto *trackmap = findNode::getClass<SvtxTrackMap>(topNode, m_trackMapName);
0516   auto *clustermap = findNode::getClass<TrkrClusterContainer>(topNode, m_clusterContainerName);
0517   if (svtxvertexmap)
0518   {
0519     m_nvertices = svtxvertexmap->size();
0520     clearClusterStateVectors();
0521
0522     for (const auto& [key, vertex] : *svtxvertexmap)
0523     {
0524       m_vertexid = key;
0525       m_vertex_crossing = vertex->get_beam_crossing();
0526       m_vx = vertex->get_x();
0527       m_vy = vertex->get_y();
0528       m_vz = vertex->get_z();
0529       m_ntracks = vertex->size_tracks();
0530
0531       for (auto it = vertex->begin_tracks(); it != vertex->end_tracks(); ++it)
0532       {
0533         auto id = *it;
0534         auto *track = trackmap->find(id)->second;
0535         if (!track)
0536         {
0537           continue;
0538         }
0539         for (const auto& ckey : get_cluster_keys(track))
0540         {
0541           TrkrCluster* cluster = clustermap->findCluster(ckey);
0542
0543           Acts::Vector3 clusglob = m_globalPositionWrapper.getGlobalPositionDistortionCorrected(key, cluster, track->get_crossing());
0544
0545           m_clusgx.push_back(clusglob.x());
0546           m_clusgy.push_back(clusglob.y());
0547           m_clusgz.push_back(clusglob.z());
0548           float clusr = r(clusglob.x(), clusglob.y());
0549           if (clusglob.y() < 0)
0550           {
0551             clusr = -clusr;
0552           }
0553           m_clusgr.push_back(clusr);
0554         }
0555       }
0556
0557       m_vertextree->Fill();
0558     }
0559   }
0560 }
0561
0562 float TrackResiduals::convertTimeToZ(ActsGeometry* geometry, TrkrDefs::cluskey cluster_key, TrkrCluster* cluster)
0563 {
0564   // must convert local Y from cluster average time of arival to local cluster z position
0565   double drift_velocity = geometry->get_drift_velocity();
0566   double zdriftlength = cluster->getLocalY() * drift_velocity;
0567   double surfCenterZ = 52.89;                // 52.89 is where G4 thinks the surface center is
0568   double zloc = surfCenterZ - zdriftlength;  // converts z drift length to local z position in the TPC in north
0569   unsigned int side = TpcDefs::getSide(cluster_key);
0570   if (side == 0)
0571   {
0572     zloc = -zloc;
0573   }
0574   float z = zloc;  // in cm
0575
0576   return z;
0577 }
0578 void TrackResiduals::circleFitClusters(
0579     std::vector<TrkrDefs::cluskey>& keys,
0580     TrkrClusterContainer* clusters,
0581     const short int& crossing)
0582 {
0583   std::vector<Acts::Vector3> clusPos;
0584   std::vector<Acts::Vector3> global_vec;
0585   for (auto& key : keys)
0586   {
0587     auto *cluster = clusters->findCluster(key);
0588     const Acts::Vector3 pos = m_globalPositionWrapper.getGlobalPositionDistortionCorrected(key, cluster, crossing);
0589     clusPos.push_back(pos);
0590   }
0591   TrackFitUtils::position_vector_t yzpoints;
0592   TrackFitUtils::position_vector_t xypoints;
0593
0594   for (auto& pos : clusPos)
0595   {
0596     float clusr = r(pos.x(), pos.y());
0597     // exclude silicon and tpot clusters for now
0598     if (std::fabs(clusr) > 80 || (m_linefitTPCOnly && std::fabs(clusr) < 20.))
0599     {
0600       continue;
0601     }
0602     xypoints.emplace_back(pos.x(), pos.y());
0603     yzpoints.emplace_back(pos.z(), pos.y());
0604     global_vec.push_back(pos);
0605   }
0606
0607   auto xyparams = TrackFitUtils::line_fit(xypoints);
0608   auto yzLineParams = TrackFitUtils::line_fit(yzpoints);
0609   auto fitpars = TrackFitUtils::fitClusters(global_vec, keys, false);
0610   // auto fitpars = TrackFitUtils::fitClusters(global_vec, keys, !m_linefitTPCOnly);
0611   m_xyint = std::get<1>(xyparams);
0612   m_xyslope = std::get<0>(xyparams);
0613   m_yzint = std::get<1>(yzLineParams);
0614   m_yzslope = std::get<0>(yzLineParams);
0615   if (!fitpars.empty())
0616   {
0617     m_R = fitpars[0];
0618     m_X0 = fitpars[1];
0619     m_Y0 = fitpars[2];
0620     m_rzslope = fitpars[3];
0621     m_rzint = fitpars[4];
0622   }
0623   else
0624   {
0625     m_R = std::numeric_limits<float>::quiet_NaN();
0626     m_X0 = std::numeric_limits<float>::quiet_NaN();
0627     m_Y0 = std::numeric_limits<float>::quiet_NaN();
0628     m_rzslope = std::numeric_limits<float>::quiet_NaN();
0629     m_rzint = std::numeric_limits<float>::quiet_NaN();
0630   }
0631 }
0632
0633 void TrackResiduals::lineFitClusters(std::vector<TrkrDefs::cluskey>& keys,
0634                                      TrkrClusterContainer* clusters,
0635                                      const short int& crossing)
0636 {
0637   std::vector<Acts::Vector3> clusPos;
0638   for (auto& key : keys)
0639   {
0640     auto *cluster = clusters->findCluster(key);
0641     const Acts::Vector3 pos = m_globalPositionWrapper.getGlobalPositionDistortionCorrected(key, cluster, crossing);
0642     clusPos.push_back(pos);
0643   }
0644   TrackFitUtils::position_vector_t xypoints;
0645   TrackFitUtils::position_vector_t rzpoints;
0646   TrackFitUtils::position_vector_t yzpoints;
0647   for (auto& pos : clusPos)
0648   {
0649     float clusr = r(pos.x(), pos.y());
0650     if (pos.y() < 0)
0651     {
0652       clusr *= -1;
0653     }
0654
0655     // exclude 1d tpot clusters for now
0656
0657     if (std::fabs(clusr) > 80 || (m_linefitTPCOnly && std::fabs(clusr) < 20.))
0658     {
0659       continue;
0660     }
0661
0662     rzpoints.emplace_back(pos.z(), clusr);
0663     xypoints.emplace_back(pos.x(), pos.y());
0664     yzpoints.emplace_back(pos.z(), pos.y());
0665   }
0666
0667   auto xyparams = TrackFitUtils::line_fit(xypoints);
0668   auto rzparams = TrackFitUtils::line_fit(rzpoints);
0669   auto yzparams = TrackFitUtils::line_fit(yzpoints);
0670   m_xyint = std::get<1>(xyparams);
0671   m_xyslope = std::get<0>(xyparams);
0672   m_rzint = std::get<1>(rzparams);
0673   m_rzslope = std::get<0>(rzparams);
0674   m_yzint = std::get<1>(yzparams);
0675   m_yzslope = std::get<0>(yzparams);
0676 }
0677
0678 void TrackResiduals::fillClusterTree(TrkrClusterContainer* clusters,
0679                                      ActsGeometry* geometry)
0680 {
0681   if (clusters->size() < m_min_cluster_size)
0682   {
0683     return;
0684   }
0685   for (const auto& det : {TrkrDefs::TrkrId::mvtxId, TrkrDefs::TrkrId::inttId,
0686                     TrkrDefs::TrkrId::tpcId, TrkrDefs::TrkrId::micromegasId})
0687   {
0688     for (const auto& hitsetkey : clusters->getHitSetKeys(det))
0689     {
0690       m_scluslayer = TrkrDefs::getLayer(hitsetkey);
0691       auto range = clusters->getClusters(hitsetkey);
0692       for (auto iter = range.first; iter != range.second; ++iter)
0693       {
0694         auto key = iter->first;
0695         auto *cluster = clusters->findCluster(key);
0696         Acts::Vector3 glob;
0697         // NOT IMPLEMENTED YET
0698         // if (TrkrDefs::getTrkrId(key) == TrkrDefs::tpcId)
0699         // {
0700         //   glob = geometry->getGlobalPosition(key, cluster);  // corrections make no sense if crossing is not known
0701         // }
0702         // else
0703         // {
0704         //   glob = geometry->getGlobalPosition(key, cluster);
0705         // }
0706         glob = geometry->getGlobalPosition(key, cluster);
0707         m_sclusgx = glob.x();
0708         m_sclusgy = glob.y();
0709         m_sclusgz = glob.z();
0710         m_sclusgr = r(m_sclusgx, m_sclusgy);
0711         m_sclusphi = atan2(glob.y(), glob.x());
0712         m_scluseta = acos(glob.z() / std::sqrt(square(glob.x()) + square(glob.y()) + square(glob.z())));
0713         m_adc = cluster->getAdc();
0714         m_clusmaxadc = cluster->getMaxAdc();
0715         m_scluslx = cluster->getLocalX();
0716         m_scluslz = cluster->getLocalY();
0717         auto para_errors = m_clusErrPara.get_clusterv5_modified_error(cluster, m_sclusgr, key);
0718         m_phisize = cluster->getPhiSize();
0719         m_zsize = cluster->getZSize();
0720         m_scluselx = std::sqrt(para_errors.first);
0721         m_scluselz = std::sqrt(para_errors.second);
0722
0723         //! Fill relevant geom info that is specific to subsystem
0724         switch (det)
0725         {
0726         case TrkrDefs::TrkrId::mvtxId:
0727           m_staveid = MvtxDefs::getStaveId(key);
0728           m_chipid = MvtxDefs::getChipId(key);
0729           m_strobeid = MvtxDefs::getStrobeId(key);
0730
0731           m_ladderzid = std::numeric_limits<int>::quiet_NaN();
0732           m_ladderphiid = std::numeric_limits<int>::quiet_NaN();
0733           m_timebucket = std::numeric_limits<int>::quiet_NaN();
0734           m_clussector = std::numeric_limits<int>::quiet_NaN();
0735           m_side = std::numeric_limits<int>::quiet_NaN();
0736           m_segtype = std::numeric_limits<int>::quiet_NaN();
0737           m_tileid = std::numeric_limits<int>::quiet_NaN();
0738           break;
0739         case TrkrDefs::TrkrId::inttId:
0740           m_ladderzid = InttDefs::getLadderZId(key);
0741           m_ladderphiid = InttDefs::getLadderPhiId(key);
0742           m_timebucket = InttDefs::getTimeBucketId(key);
0743
0744           m_staveid = std::numeric_limits<int>::quiet_NaN();
0745           m_chipid = std::numeric_limits<int>::quiet_NaN();
0746           m_strobeid = std::numeric_limits<int>::quiet_NaN();
0747           m_clussector = std::numeric_limits<int>::quiet_NaN();
0748           m_side = std::numeric_limits<int>::quiet_NaN();
0749           m_segtype = std::numeric_limits<int>::quiet_NaN();
0750           m_tileid = std::numeric_limits<int>::quiet_NaN();
0751           break;
0752         case TrkrDefs::TrkrId::tpcId:
0753           m_clussector = TpcDefs::getSectorId(key);
0754           m_side = TpcDefs::getSide(key);
0755
0756           m_staveid = std::numeric_limits<int>::quiet_NaN();
0757           m_chipid = std::numeric_limits<int>::quiet_NaN();
0758           m_strobeid = std::numeric_limits<int>::quiet_NaN();
0759           m_ladderzid = std::numeric_limits<int>::quiet_NaN();
0760           m_ladderphiid = std::numeric_limits<int>::quiet_NaN();
0761           m_timebucket = std::numeric_limits<int>::quiet_NaN();
0762           m_segtype = std::numeric_limits<int>::quiet_NaN();
0763           m_tileid = std::numeric_limits<int>::quiet_NaN();
0764           break;
0765         case TrkrDefs::TrkrId::micromegasId:
0766           m_segtype = (int) MicromegasDefs::getSegmentationType(key);
0767           m_tileid = MicromegasDefs::getTileId(key);
0768
0769           m_staveid = std::numeric_limits<int>::quiet_NaN();
0770           m_chipid = std::numeric_limits<int>::quiet_NaN();
0771           m_strobeid = std::numeric_limits<int>::quiet_NaN();
0772           m_ladderzid = std::numeric_limits<int>::quiet_NaN();
0773           m_ladderphiid = std::numeric_limits<int>::quiet_NaN();
0774           m_timebucket = std::numeric_limits<int>::quiet_NaN();
0775           m_clussector = std::numeric_limits<int>::quiet_NaN();
0776           m_side = std::numeric_limits<int>::quiet_NaN();
0777           break;
0778         default:
0779           break;
0780         }
0781
0782         m_clustree->Fill();
0783       }
0784     }
0785   }
0786 }
0787
0788 //____________________________________________________________________________..
0789 int TrackResiduals::End(PHCompositeNode* /*unused*/)
0790 {
0791   m_outfile->cd();
0792   m_tree->Write();
0793   if (m_doClusters)
0794   {
0795     m_clustree->Write();
0796   }
0797   if (m_doHits)
0798   {
0799     m_hittree->Write();
0800   }
0801   if (m_doVertex)
0802   {
0803     m_vertextree->Write();
0804   }
0805   if (m_doFailedSeeds)
0806   {
0807     m_failedfits->Write();
0808   }
0809   if (m_doEventTree)
0810   {
0811     m_eventtree->Write();
0812   }
0813   m_outfile->Close();
0814
0815   return Fun4AllReturnCodes::EVENT_OK;
0816 }
0817 void TrackResiduals::fillHitTree(TrkrHitSetContainer* hitmap,
0818                                  ActsGeometry* geometry,
0819                                  PHG4TpcCylinderGeomContainer* tpcGeom,
0820                                  PHG4CylinderGeomContainer* mvtxGeom,
0821                                  PHG4CylinderGeomContainer* inttGeom,
0822                                  PHG4CylinderGeomContainer* mmGeom)
0823 {
0824   if (!tpcGeom or !mvtxGeom or !inttGeom or !mmGeom)
0825   {
0826     std::cout << PHWHERE << "missing hit map, can't continue with hit tree"
0827               << std::endl;
0828     return;
0829   }
0830   TrkrHitSetContainer::ConstRange all_hitsets = hitmap->getHitSets();
0831   for (TrkrHitSetContainer::ConstIterator hitsetiter = all_hitsets.first;
0832        hitsetiter != all_hitsets.second;
0833        ++hitsetiter)
0834   {
0835     m_hitsetkey = hitsetiter->first;
0836     TrkrHitSet* hitset = hitsetiter->second;
0837
0838     m_hitlayer = TrkrDefs::getLayer(m_hitsetkey);
0839     auto det = TrkrDefs::getTrkrId(m_hitsetkey);
0840     //! Fill relevant geom info that is specific to subsystem
0841     switch (det)
0842     {
0843     case TrkrDefs::TrkrId::mvtxId:
0844     {
0845       m_staveid = MvtxDefs::getStaveId(m_hitsetkey);
0846       m_chipid = MvtxDefs::getChipId(m_hitsetkey);
0847       m_strobeid = MvtxDefs::getStrobeId(m_hitsetkey);
0848
0849       m_ladderzid = std::numeric_limits<int>::quiet_NaN();
0850       m_ladderphiid = std::numeric_limits<int>::quiet_NaN();
0851       m_timebucket = std::numeric_limits<int>::quiet_NaN();
0852       m_sector = std::numeric_limits<int>::quiet_NaN();
0853       m_side = std::numeric_limits<int>::quiet_NaN();
0854       m_segtype = std::numeric_limits<int>::quiet_NaN();
0855       m_tileid = std::numeric_limits<int>::quiet_NaN();
0856       break;
0857     }
0858     case TrkrDefs::TrkrId::inttId:
0859     {
0860       m_ladderzid = InttDefs::getLadderZId(m_hitsetkey);
0861       m_ladderphiid = InttDefs::getLadderPhiId(m_hitsetkey);
0862       m_timebucket = InttDefs::getTimeBucketId(m_hitsetkey);
0863
0864       m_staveid = std::numeric_limits<int>::quiet_NaN();
0865       m_chipid = std::numeric_limits<int>::quiet_NaN();
0866       m_strobeid = std::numeric_limits<int>::quiet_NaN();
0867       m_sector = std::numeric_limits<int>::quiet_NaN();
0868       m_side = std::numeric_limits<int>::quiet_NaN();
0869       m_segtype = std::numeric_limits<int>::quiet_NaN();
0870       m_tileid = std::numeric_limits<int>::quiet_NaN();
0871       break;
0872     }
0873     case TrkrDefs::TrkrId::tpcId:
0874     {
0875       m_sector = TpcDefs::getSectorId(m_hitsetkey);
0876       m_side = TpcDefs::getSide(m_hitsetkey);
0877
0878       m_staveid = std::numeric_limits<int>::quiet_NaN();
0879       m_chipid = std::numeric_limits<int>::quiet_NaN();
0880       m_strobeid = std::numeric_limits<int>::quiet_NaN();
0881       m_ladderzid = std::numeric_limits<int>::quiet_NaN();
0882       m_ladderphiid = std::numeric_limits<int>::quiet_NaN();
0883       m_timebucket = std::numeric_limits<int>::quiet_NaN();
0884       m_segtype = std::numeric_limits<int>::quiet_NaN();
0885       m_tileid = std::numeric_limits<int>::quiet_NaN();
0886
0887       break;
0888     }
0889     case TrkrDefs::TrkrId::micromegasId:
0890     {
0891       m_segtype = (int) MicromegasDefs::getSegmentationType(m_hitsetkey);
0892       m_tileid = MicromegasDefs::getTileId(m_hitsetkey);
0893
0894       m_staveid = std::numeric_limits<int>::quiet_NaN();
0895       m_chipid = std::numeric_limits<int>::quiet_NaN();
0896       m_strobeid = std::numeric_limits<int>::quiet_NaN();
0897       m_ladderzid = std::numeric_limits<int>::quiet_NaN();
0898       m_ladderphiid = std::numeric_limits<int>::quiet_NaN();
0899       m_timebucket = std::numeric_limits<int>::quiet_NaN();
0900       m_sector = std::numeric_limits<int>::quiet_NaN();
0901       m_side = std::numeric_limits<int>::quiet_NaN();
0902       break;
0903     }
0904     default:
0905       break;
0906     }
0907
0908     // Got all stave/ladder/sector/tile info, now get the actual hit info
0909     auto hitrangei = hitset->getHits();
0910     for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
0911          hitr != hitrangei.second;
0912          ++hitr)
0913     {
0914       auto hitkey = hitr->first;
0915       auto *hit = hitr->second;
0916       m_adc = hit->getAdc();
0917
0918       switch (det)
0919       {
0920       case TrkrDefs::TrkrId::mvtxId:
0921       {
0922         m_row = MvtxDefs::getRow(hitkey);
0923         m_col = MvtxDefs::getCol(hitkey);
0924         auto *layergeom = dynamic_cast<CylinderGeom_Mvtx*>(mvtxGeom->GetLayerGeom(m_hitlayer));
0925         auto local_coords = layergeom->get_local_coords_from_pixel(m_row, m_col);
0926         TVector2 local;
0927         local.SetX(local_coords.X());
0928         local.SetY(local_coords.Z());
0929         auto surf = geometry->maps().getSiliconSurface(m_hitsetkey);
0930         auto glob = CylinderGeom_MvtxHelper::get_world_from_local_coords(surf, geometry, local);
0931         m_hitgx = glob.X();
0932         m_hitgy = glob.Y();
0933         m_hitgz = glob.Z();
0934
0935         m_segtype = std::numeric_limits<int>::quiet_NaN();
0936         m_tileid = std::numeric_limits<int>::quiet_NaN();
0937         m_strip = std::numeric_limits<int>::quiet_NaN();
0938         m_hitpad = std::numeric_limits<int>::quiet_NaN();
0939         m_hittbin = std::numeric_limits<int>::quiet_NaN();
0940
0941         m_zdriftlength = std::numeric_limits<float>::quiet_NaN();
0942         break;
0943       }
0944       case TrkrDefs::TrkrId::inttId:
0945       {
0946         m_row = InttDefs::getRow(hitkey);
0947         m_col = InttDefs::getCol(hitkey);
0948         auto *geom = dynamic_cast<CylinderGeomIntt*>(inttGeom->GetLayerGeom(m_hitlayer));
0949         double local_hit_loc[3] = {0, 0, 0};
0950         geom->find_strip_center_localcoords(m_ladderzid, m_row, m_col, local_hit_loc);
0951         auto surf = geometry->maps().getSiliconSurface(m_hitsetkey);
0952         TVector2 local;
0953         local.SetX(local_hit_loc[1]);
0954         local.SetY(local_hit_loc[2]);
0955         auto glob = CylinderGeomInttHelper::get_world_from_local_coords(surf, geometry, local);
0956
0957         m_hitgx = glob.X();
0958         m_hitgy = glob.Y();
0959         m_hitgz = glob.Z();
0960         m_segtype = std::numeric_limits<int>::quiet_NaN();
0961         m_tileid = std::numeric_limits<int>::quiet_NaN();
0962         m_strip = std::numeric_limits<int>::quiet_NaN();
0963         m_hitpad = std::numeric_limits<int>::quiet_NaN();
0964         m_hittbin = std::numeric_limits<int>::quiet_NaN();
0965         m_zdriftlength = std::numeric_limits<float>::quiet_NaN();
0966         break;
0967       }
0968       case TrkrDefs::TrkrId::tpcId:
0969       {
0970         m_row = std::numeric_limits<int>::quiet_NaN();
0971         m_col = std::numeric_limits<int>::quiet_NaN();
0972         m_segtype = std::numeric_limits<int>::quiet_NaN();
0973         m_tileid = std::numeric_limits<int>::quiet_NaN();
0974         m_strip = std::numeric_limits<int>::quiet_NaN();
0975
0976         m_hitpad = TpcDefs::getPad(hitkey);
0977         m_hittbin = TpcDefs::getTBin(hitkey);
0978
0979         auto *geoLayer = tpcGeom->GetLayerCellGeom(m_hitlayer);
0980         auto phi = geoLayer->get_phicenter(m_hitpad, m_side);
0981         auto radius = geoLayer->get_radius();
0982         float AdcClockPeriod = geoLayer->get_zstep();
0983         auto glob = geometry->getGlobalPositionTpc(m_hitsetkey, hitkey, phi, radius, AdcClockPeriod);
0984         m_hitgx = glob.x();
0985         m_hitgy = glob.y();
0986         m_hitgz = glob.z();
0987
0988         break;
0989       }
0990       case TrkrDefs::TrkrId::micromegasId:
0991       {
0992         auto *const layergeom = dynamic_cast<CylinderGeomMicromegas*>(mmGeom->GetLayerGeom(m_hitlayer));
0993         m_strip = MicromegasDefs::getStrip(hitkey);
0994         const auto global_coord = layergeom->get_world_coordinates(m_tileid, geometry, m_strip);
0995         m_hitgx = global_coord.X();
0996         m_hitgy = global_coord.Y();
0997         m_hitgz = global_coord.Z();
0998         m_row = std::numeric_limits<int>::quiet_NaN();
0999         m_col = std::numeric_limits<int>::quiet_NaN();
1000         m_segtype = std::numeric_limits<int>::quiet_NaN();
1001         m_tileid = std::numeric_limits<int>::quiet_NaN();
1002         m_hitpad = std::numeric_limits<int>::quiet_NaN();
1003         m_hittbin = std::numeric_limits<int>::quiet_NaN();
1004
1005         m_zdriftlength = std::numeric_limits<float>::quiet_NaN();
1006       }
1007       default:
1008         break;
1009       }
1010
1011       m_hittree->Fill();
1012     }
1013   }
1014 }
1015
1016 void TrackResiduals::fillClusterBranchesKF(TrkrDefs::cluskey ckey, SvtxTrack* track,
1017                                            const std::vector<std::pair<TrkrDefs::cluskey, Acts::Vector3>>& global,
1018                                            PHCompositeNode* topNode)
1019 {
1020   auto *clustermap = findNode::getClass<TrkrClusterContainer>(topNode, m_clusterContainerName);
1021   auto *geometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
1022
1023   auto global_moved = global;  // if use transient transforms for distortion correction
1024   if (m_use_clustermover)
1025   {
1026     // move the corrected cluster positions back to the original readout surface
1027     global_moved = m_clusterMover.processTrack(global);
1028   }
1029
1030   ActsTransformations transformer;
1031   TrkrCluster* cluster = clustermap->findCluster(ckey);
1032
1033   // loop over global vectors and get this cluster
1034   Acts::Vector3 clusglob(0, 0, 0);
1035   for (const auto& pair : global)
1036   {
1037     auto thiskey = pair.first;
1038     clusglob = pair.second;
1039     if (thiskey == ckey)
1040     {
1041       break;
1042     }
1043   }
1044
1045   Acts::Vector3 clusglob_moved(0, 0, 0);
1046   for (const auto& pair : global_moved)
1047   {
1048     auto thiskey = pair.first;
1049     clusglob_moved = pair.second;
1050     if (thiskey == ckey)
1051     {
1052       break;
1053     }
1054   }
1055
1056   unsigned int layer = TrkrDefs::getLayer(ckey);
1057
1058   if (Verbosity() > 1)
1059   {
1060     std::cout << "Called :fillClusterBranchesKF for ckey " << ckey << " layer " << layer << " trackid " << track->get_id() << " clusglob x " << clusglob(0) << " y " << clusglob(1) << " z " << clusglob(2) << std::endl;
1061   }
1062
1063   switch (TrkrDefs::getTrkrId(ckey))
1064   {
1065   case TrkrDefs::mvtxId:
1066     m_nmaps++;
1067     break;
1068   case TrkrDefs::inttId:
1069     m_nintt++;
1070     break;
1071   case TrkrDefs::tpcId:
1072     m_ntpc++;
1073     m_clsector.push_back(TpcDefs::getSectorId(ckey));
1074     m_clside.push_back(TpcDefs::getSide(ckey));
1075     break;
1076   case TrkrDefs::micromegasId:
1077     m_nmms++;
1078     m_tileid = MicromegasDefs::getTileId(ckey);
1079     break;
1080   default:
1081     std::cout << PHWHERE << " unknown key " << ckey << std::endl;
1082     gSystem->Exit(1);
1083     exit(1);
1084   }
1085
1086   SvtxTrackState* state = nullptr;
1087
1088   // the track states from the Acts fit are fitted to fully corrected clusters, and are on the surface
1089   for (auto state_iter = track->begin_states();
1090        state_iter != track->end_states();
1091        ++state_iter)
1092   {
1093     SvtxTrackState* tstate = state_iter->second;
1094     auto stateckey = tstate->get_cluskey();
1095     if (stateckey == ckey)
1096     {
1097       state = tstate;
1098       break;
1099     }
1100   }
1101
1102   if (!state)
1103   {
1104     if (Verbosity() > 1)
1105     {
1106       std::cout << "   no state for cluster " << ckey << "  in layer " << layer << std::endl;
1107     }
1108   }
1109   else
1110   {
1111     switch (TrkrDefs::getTrkrId(ckey))
1112     {
1113     case TrkrDefs::mvtxId:
1114       m_nmapsstate++;
1115       break;
1116     case TrkrDefs::inttId:
1117       m_ninttstate++;
1118       break;
1119     case TrkrDefs::tpcId:
1120       m_ntpcstate++;
1121       break;
1122     case TrkrDefs::micromegasId:
1123       m_nmmsstate++;
1124       break;
1125     default:
1126       std::cout << PHWHERE << " unknown key " << ckey << std::endl;
1127       gSystem->Exit(1);
1128       exit(1);
1129     }
1130   }
1131
1132   m_cluskeys.push_back(ckey);
1133
1134   //! have cluster and state, fill vectors
1135   m_clusedge.push_back(cluster->getEdge());
1136   m_clusoverlap.push_back(cluster->getOverlap());
1137
1138   // get new local coords from moved cluster
1139   Surface surf = geometry->maps().getSurface(ckey, cluster);
1140   Surface surf_ideal = geometry->maps().getSurface(ckey, cluster);  // Unchanged by distortion corrections
1141   // if this is a TPC cluster, the crossing correction may have moved it across the central membrane, check the surface
1142   auto trkrid = TrkrDefs::getTrkrId(ckey);
1143   if (trkrid == TrkrDefs::tpcId)
1144   {
1145     TrkrDefs::hitsetkey hitsetkey = TrkrDefs::getHitSetKeyFromClusKey(ckey);
1146     TrkrDefs::subsurfkey new_subsurfkey = 0;
1147     surf = geometry->get_tpc_surface_from_coords(hitsetkey, clusglob_moved, new_subsurfkey);
1148   }
1149   if (!surf)
1150   {
1151     if (Verbosity() > 2)
1152     {
1153       std::cout << " Failed to find surface for cluskey " << ckey << std::endl;
1154     }
1155     return;
1156   }
1157
1158   // get local coordinates
1159   Acts::Vector2 loc;
1160   loc = geometry->getLocalCoords(ckey, cluster, m_crossing);
1161   if (m_use_clustermover)
1162   {
1163     // in this case we get local coords from transform of corrected global coords
1164     clusglob_moved *= Acts::UnitConstants::cm;  // we want mm for transformations
1165     Acts::Vector3 normal = surf->normal(geometry->geometry().getGeoContext(),
1166                                         Acts::Vector3(1, 1, 1), Acts::Vector3(1, 1, 1));
1167     auto local = surf->globalToLocal(geometry->geometry().getGeoContext(),
1168                                      clusglob_moved, normal);
1169     if (local.ok())
1170     {
1171       loc = local.value() / Acts::UnitConstants::cm;
1172     }
1173     else
1174     {
1175       // otherwise take the manual calculation for the TPC
1176       // doing it this way just avoids the bounds check that occurs in the surface class method
1177       Acts::Vector3 loct = surf->transform(geometry->geometry().getGeoContext()).inverse() * clusglob_moved;  // global is in mm
1178       loct /= Acts::UnitConstants::cm;
1179
1180       loc(0) = loct(0);
1181       loc(1) = loct(1);
1182     }
1183     clusglob_moved /= Acts::UnitConstants::cm;  // we want cm for the tree
1184   }
1185
1186   m_cluslx.push_back(loc.x());
1187   m_cluslz.push_back(loc.y());
1188
1189   if (Verbosity() > 2)
1190   {
1191     std::cout << "Trackresiduals cluster (cm): localX " << loc.x() << " localY " << loc.y() << std::endl
1192               << " global.x " << clusglob_moved(0) << " global.y " << clusglob_moved(1) << " global.z " << clusglob_moved(2) << std::endl;
1193   }
1194
1195   float clusr = r(clusglob_moved.x(), clusglob_moved.y());
1196   auto para_errors = m_clusErrPara.get_clusterv5_modified_error(cluster,
1197                                                                 clusr, ckey);
1198   m_cluselx.push_back(sqrt(para_errors.first));
1199   m_cluselz.push_back(sqrt(para_errors.second));
1200   m_clusgx.push_back(clusglob_moved.x());
1201   m_clusgy.push_back(clusglob_moved.y());
1202   m_clusgr.push_back(clusglob_moved.y() > 0 ? clusr : -1 * clusr);
1203   m_clusgz.push_back(clusglob_moved.z());
1204   m_clusgxunmoved.push_back(clusglob.x());
1205   m_clusgyunmoved.push_back(clusglob.y());
1206   m_clusgzunmoved.push_back(clusglob.z());
1207   m_clusAdc.push_back(cluster->getAdc());
1208   m_clusMaxAdc.push_back(cluster->getMaxAdc());
1209   m_cluslayer.push_back(TrkrDefs::getLayer(ckey));
1210   m_clusphisize.push_back(cluster->getPhiSize());
1211   m_cluszsize.push_back(cluster->getZSize());
1212
1213   auto misaligncenter = surf->center(geometry->geometry().getGeoContext());
1214   auto misalignnorm = -1 * surf->normal(geometry->geometry().getGeoContext(), Acts::Vector3(1, 1, 1), Acts::Vector3(1, 1, 1));
1215   auto misrot = surf->transform(geometry->geometry().getGeoContext()).rotation();
1216
1217   float mgamma = atan2(-misrot(1, 0), misrot(0, 0));
1218   float mbeta = -asin(misrot(0, 1));
1219   float malpha = atan2(misrot(1, 1), misrot(2, 1));
1220
1221   //! Switch to get ideal transforms
1222   alignmentTransformationContainer::use_alignment = false;
1223   auto idealcenter = surf_ideal->center(geometry->geometry().getGeoContext());
1224   auto idealnorm = -1 * surf_ideal->normal(geometry->geometry().getGeoContext(), Acts::Vector3(1, 1, 1), Acts::Vector3(1, 1, 1));
1225
1226   // replace the corrected moved cluster local position with the readout position from ideal geometry for now
1227   // This allows us to see the distortion corrections by subtracting this uncorrected position
1228   // revisit this when looking at the alignment case
1229   //  Acts::Vector3 ideal_local(loc.x(), loc.y(), 0.0);
1230   auto nominal_loc = geometry->getLocalCoords(ckey, cluster);
1231   Acts::Vector3 ideal_local(nominal_loc.x(), nominal_loc.y(), 0.0);
1232   Acts::Vector3 ideal_glob = surf_ideal->transform(geometry->geometry().getGeoContext()) * (ideal_local * Acts::UnitConstants::cm);
1233   auto idealrot = surf_ideal->transform(geometry->geometry().getGeoContext()).rotation();
1234
1235   //! These calculations are taken from the wikipedia page for Euler angles,
1236   //! under the Tait-Bryan angle explanation. Formulas for the angles
1237   //! calculated from the rotation matrices depending on what order the
1238   //! rotation matrix is constructed are given
1239   //! They need to be modified to conform to the Acts basis of (x,z,y), for
1240   //! which the wiki page expects (x,y,z). This includes swapping the sign
1241   //! of some elements to account for the permutation
1242   //! https://en.wikipedia.org/wiki/Euler_angles#Conversion_to_other_orientation_representations
1243   float igamma = atan2(-idealrot(1, 0), idealrot(0, 0));
1244   float ibeta = -asin(idealrot(0, 1));
1245   float ialpha = atan2(idealrot(1, 1), idealrot(2, 1));
1246
1247   alignmentTransformationContainer::use_alignment = true;
1248
1249   idealcenter /= Acts::UnitConstants::cm;
1250   misaligncenter /= Acts::UnitConstants::cm;
1251   ideal_glob /= Acts::UnitConstants::cm;
1252
1253   m_idealsurfalpha.push_back(ialpha);
1254   m_idealsurfbeta.push_back(ibeta);
1255   m_idealsurfgamma.push_back(igamma);
1256   m_missurfalpha.push_back(malpha);
1257   m_missurfbeta.push_back(mbeta);
1258   m_missurfgamma.push_back(mgamma);
1259
1260   m_idealsurfcenterx.push_back(idealcenter.x());
1261   m_idealsurfcentery.push_back(idealcenter.y());
1262   m_idealsurfcenterz.push_back(idealcenter.z());
1263   m_idealsurfnormx.push_back(idealnorm.x());
1264   m_idealsurfnormy.push_back(idealnorm.y());
1265   m_idealsurfnormz.push_back(idealnorm.z());
1266   m_missurfcenterx.push_back(misaligncenter.x());
1267   m_missurfcentery.push_back(misaligncenter.y());
1268   m_missurfcenterz.push_back(misaligncenter.z());
1269   m_missurfnormx.push_back(misalignnorm.x());
1270   m_missurfnormy.push_back(misalignnorm.y());
1271   m_missurfnormz.push_back(misalignnorm.z());
1272   m_clusgxideal.push_back(ideal_glob.x());
1273   m_clusgyideal.push_back(ideal_glob.y());
1274   m_clusgzideal.push_back(ideal_glob.z());
1275
1276   if (state)
1277   {
1278     Acts::Vector3 stateglob(state->get_x(), state->get_y(), state->get_z());
1279     Acts::Vector2 stateloc(state->get_localX(), state->get_localY());
1280
1281     if (Verbosity() > 2)
1282     {
1283       std::cout << "Trackresiduals state (cm): localX " << stateloc(0) << " localY " << stateloc(1) << std::endl
1284                 << " stateglobx " << stateglob(0) << " stategloby " << stateglob(1) << " stateglobz " << stateglob(2) << std::endl;
1285     }
1286
1287     const auto actscov =
1288         transformer.rotateSvtxTrackCovToActs(state);
1289
1290     m_statelx.push_back(stateloc(0));
1291     m_statelz.push_back(stateloc(1));
1292     m_stateelx.push_back(std::sqrt(actscov(Acts::eBoundLoc0, Acts::eBoundLoc0)) / Acts::UnitConstants::cm);
1293     m_stateelz.push_back(std::sqrt(actscov(Acts::eBoundLoc1, Acts::eBoundLoc1)) / Acts::UnitConstants::cm);
1294     m_stategx.push_back(state->get_x());
1295     m_stategy.push_back(state->get_y());
1296     m_stategz.push_back(state->get_z());
1297     m_statepx.push_back(state->get_px());
1298     m_statepy.push_back(state->get_py());
1299     m_statepz.push_back(state->get_pz());
1300     m_statepl.push_back(state->get_pathlength());
1301   }
1302   else
1303   {
1304     // cluster has no corresponding state, set state variables to NaNs
1305     m_statelx.push_back(std::numeric_limits<float>::quiet_NaN());
1306     m_statelz.push_back(std::numeric_limits<float>::quiet_NaN());
1307     m_stategx.push_back(std::numeric_limits<float>::quiet_NaN());
1308     m_stategy.push_back(std::numeric_limits<float>::quiet_NaN());
1309     m_stategz.push_back(std::numeric_limits<float>::quiet_NaN());
1310     m_statepx.push_back(std::numeric_limits<float>::quiet_NaN());
1311     m_statepy.push_back(std::numeric_limits<float>::quiet_NaN());
1312     m_statepz.push_back(std::numeric_limits<float>::quiet_NaN());
1313     m_statepl.push_back(std::numeric_limits<float>::quiet_NaN());
1314   }
1315
1316   if (Verbosity() > 2)
1317   {
1318     if (ideal_glob(2) > 0)
1319     {
1320       double xideal = ideal_glob.x();
1321       double yideal = ideal_glob.y();
1322       double zideal = ideal_glob.z();
1323
1324       double xunmoved = clusglob.x();
1325       double yunmoved = clusglob.y();
1326       double zunmoved = clusglob.z();
1327
1328       double xmoved = clusglob_moved.x();
1329       double ymoved = clusglob_moved.y();
1330       double zmoved = clusglob_moved.z();
1331
1332       double this_radius_ideal = sqrt((xideal * xideal) + (yideal * yideal));
1333       double this_radius_unmoved = sqrt((xunmoved * xunmoved) + (yunmoved * yunmoved));
1334
1335       double this_phi_unmoved = atan2(yunmoved, xunmoved);
1336       double this_phi_ideal = atan2(yideal, xideal);
1337       double this_phi_moved = atan2(ymoved, xmoved);
1338
1339       std::cout << " global:  unmoved " << xunmoved << "  " << yunmoved << "  " << zunmoved << " this_phi " << this_phi_unmoved << " radius_unmoved " << this_radius_unmoved
1340                 << " r*phi " << this_radius_ideal * this_phi_unmoved << std::endl;
1341       std::cout << "              global: ideal " << xideal << "  " << yideal << "  " << zideal << " this_phi_ideal " << this_phi_ideal << " radius_ideal " << this_radius_ideal
1342                 << " r*phi " << this_radius_ideal * this_phi_ideal << std::endl;
1343       std::cout << "              global: moved " << xmoved << "  " << ymoved << "  " << zmoved << " phi_moved " << this_phi_moved
1344                 << " r*phi " << this_radius_ideal * this_phi_moved << std::endl;
1345       std::cout << "                        d_radius " << this_radius_unmoved - this_radius_ideal << " clusgz " << zideal << " r*dphi " << this_radius_ideal * (this_phi_unmoved - this_phi_ideal) << std::endl;
1346     }
1347   }
1348 }
1349
1350 void TrackResiduals::fillClusterBranchesSeeds(TrkrDefs::cluskey ckey,  // SvtxTrack* track,
1351                                               const std::vector<std::pair<TrkrDefs::cluskey, Acts::Vector3>>& global,
1352                                               PHCompositeNode* topNode)
1353 {
1354   // The input map global contains the corrected cluster positions - NOT moved back to the surfacer.
1355   // When filling the residualtree:
1356   //    clusgx etc are the corrected - but not moved back to the surface - cluster positions
1357   //    clugxideal etc are the completely uncorrected cluster positions - they do not even have crossing corrections
1358   // CircleFitClusters is called in this method. It applies TOF, crossing, and all distortion corrections before fitting
1359   //    stategx etc are at the intersection point of the helical fit with the cluster surface
1360
1361   auto *clustermap = findNode::getClass<TrkrClusterContainer>(topNode, m_clusterContainerName);
1362   auto *geometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
1363
1364   auto global_moved = global;
1365   if (m_use_clustermover)
1366   {
1367     // move the corrected cluster positions back to the original readout surface
1368     global_moved = m_clusterMover.processTrack(global);
1369   }
1370
1371   TrkrCluster* cluster = clustermap->findCluster(ckey);
1372
1373   // loop over global vectors and get this cluster
1374   Acts::Vector3 clusglob(0, 0, 0);
1375   for (const auto& pair : global)
1376   {
1377     auto thiskey = pair.first;
1378     clusglob = pair.second;
1379     if (thiskey == ckey)
1380     {
1381       break;
1382     }
1383   }
1384   Acts::Vector3 clusglob_moved(0, 0, 0);
1385   for (const auto& pair : global_moved)
1386   {
1387     auto thiskey = pair.first;
1388     clusglob_moved = pair.second;
1389     if (thiskey == ckey)
1390     {
1391       break;
1392     }
1393   }
1394
1395   switch (TrkrDefs::getTrkrId(ckey))
1396   {
1397   case TrkrDefs::mvtxId:
1398     m_nmaps++;
1399     m_clsector.push_back(-1);
1400     m_clside.push_back(-1);
1401     break;
1402   case TrkrDefs::inttId:
1403     m_clsector.push_back(-1);
1404     m_clside.push_back(-1);
1405     m_nintt++;
1406     break;
1407   case TrkrDefs::tpcId:
1408     m_ntpc++;
1409     m_clsector.push_back(TpcDefs::getSectorId(ckey));
1410     m_clside.push_back(TpcDefs::getSide(ckey));
1411     break;
1412   case TrkrDefs::micromegasId:
1413     m_nmms++;
1414     m_tileid = MicromegasDefs::getTileId(ckey);
1415     m_clsector.push_back(-1);
1416     m_clside.push_back(-1);
1417     break;
1418   default:
1419     std::cout << PHWHERE << " unknown key " << ckey << std::endl;
1420     gSystem->Exit(1);
1421     exit(1);
1422   }
1423
1424   m_cluskeys.push_back(ckey);
1425
1426   if (Verbosity() > 1)
1427   {
1428     std::cout << "Called fillClusterBranchesSeeds for ckey " << ckey << " clusglob x " << clusglob(0) << " y " << clusglob(1) << " z " << clusglob(2) << std::endl;
1429   }
1430
1431   //! have cluster and state, fill vectors
1432   m_clusedge.push_back(cluster->getEdge());
1433   m_clusoverlap.push_back(cluster->getOverlap());
1434
1435   // This is the nominal position of the cluster in local coords, completely uncorrected - is that what we want?
1436   auto loc = geometry->getLocalCoords(ckey, cluster);
1437
1438   m_cluslx.push_back(loc.x());
1439   m_cluslz.push_back(loc.y());
1440
1441   float clusr = r(clusglob_moved.x(), clusglob_moved.y());
1442   auto para_errors = m_clusErrPara.get_clusterv5_modified_error(cluster,
1443                                                                 clusr, ckey);
1444   m_cluselx.push_back(sqrt(para_errors.first));
1445   m_cluselz.push_back(sqrt(para_errors.second));
1446   m_clusgx.push_back(clusglob_moved.x());
1447   m_clusgy.push_back(clusglob_moved.y());
1448   m_clusgr.push_back(clusglob_moved.y() > 0 ? clusr : -1 * clusr);
1449   m_clusgz.push_back(clusglob_moved.z());
1450   m_clusgxunmoved.push_back(clusglob.x());
1451   m_clusgyunmoved.push_back(clusglob.y());
1452   m_clusgzunmoved.push_back(clusglob.z());
1453   m_clusAdc.push_back(cluster->getAdc());
1454   m_clusMaxAdc.push_back(cluster->getMaxAdc());
1455   m_cluslayer.push_back(TrkrDefs::getLayer(ckey));
1456   m_clusphisize.push_back(cluster->getPhiSize());
1457   m_cluszsize.push_back(cluster->getZSize());
1458
1459   if (Verbosity() > 1)
1460   {
1461     std::cout << "Track state/clus in layer "
1462               << (unsigned int) TrkrDefs::getLayer(ckey) << " with pos "
1463               << clusglob.transpose() << std::endl;
1464   }
1465
1466   auto surf = geometry->maps().getSurface(ckey, cluster);
1467
1468   auto misaligncenter = surf->center(geometry->geometry().getGeoContext());
1469   auto misalignnorm = -1 * surf->normal(geometry->geometry().getGeoContext(), Acts::Vector3(1, 1, 1), Acts::Vector3(1, 1, 1));
1470   auto misrot = surf->transform(geometry->geometry().getGeoContext()).rotation();
1471
1472   float mgamma = atan2(-misrot(1, 0), misrot(0, 0));
1473   float mbeta = -asin(misrot(0, 1));
1474   float malpha = atan2(misrot(1, 1), misrot(2, 1));
1475
1476   //! Switch to get ideal transforms
1477   alignmentTransformationContainer::use_alignment = false;
1478   auto idealcenter = surf->center(geometry->geometry().getGeoContext());
1479   auto idealnorm = -1 * surf->normal(geometry->geometry().getGeoContext(), Acts::Vector3(1, 1, 1), Acts::Vector3(1, 1, 1));
1480   Acts::Vector3 ideal_local(loc.x(), loc.y(), 0.0);
1481   Acts::Vector3 ideal_glob = surf->transform(geometry->geometry().getGeoContext()) * (ideal_local * Acts::UnitConstants::cm);
1482   auto idealrot = surf->transform(geometry->geometry().getGeoContext()).rotation();
1483
1484   //! These calculations are taken from the wikipedia page for Euler angles,
1485   //! under the Tait-Bryan angle explanation. Formulas for the angles
1486   //! calculated from the rotation matrices depending on what order the
1487   //! rotation matrix is constructed are given
1488   //! They need to be modified to conform to the Acts basis of (x,z,y), for
1489   //! which the wiki page expects (x,y,z). This includes swapping the sign
1490   //! of some elements to account for the permutation
1491   //! https://en.wikipedia.org/wiki/Euler_angles#Conversion_to_other_orientation_representations
1492   float igamma = atan2(-idealrot(1, 0), idealrot(0, 0));
1493   float ibeta = -asin(idealrot(0, 1));
1494   float ialpha = atan2(idealrot(1, 1), idealrot(2, 1));
1495
1496   alignmentTransformationContainer::use_alignment = true;
1497
1498   idealcenter /= Acts::UnitConstants::cm;
1499   misaligncenter /= Acts::UnitConstants::cm;
1500   ideal_glob /= Acts::UnitConstants::cm;
1501
1502   m_idealsurfalpha.push_back(ialpha);
1503   m_idealsurfbeta.push_back(ibeta);
1504   m_idealsurfgamma.push_back(igamma);
1505   m_missurfalpha.push_back(malpha);
1506   m_missurfbeta.push_back(mbeta);
1507   m_missurfgamma.push_back(mgamma);
1508
1509   m_idealsurfcenterx.push_back(idealcenter.x());
1510   m_idealsurfcentery.push_back(idealcenter.y());
1511   m_idealsurfcenterz.push_back(idealcenter.z());
1512   m_idealsurfnormx.push_back(idealnorm.x());
1513   m_idealsurfnormy.push_back(idealnorm.y());
1514   m_idealsurfnormz.push_back(idealnorm.z());
1515   m_missurfcenterx.push_back(misaligncenter.x());
1516   m_missurfcentery.push_back(misaligncenter.y());
1517   m_missurfcenterz.push_back(misaligncenter.z());
1518   m_missurfnormx.push_back(misalignnorm.x());
1519   m_missurfnormy.push_back(misalignnorm.y());
1520   m_missurfnormz.push_back(misalignnorm.z());
1521   m_clusgxideal.push_back(ideal_glob.x());
1522   m_clusgyideal.push_back(ideal_glob.y());
1523   m_clusgzideal.push_back(ideal_glob.z());
1524
1525   if (m_zeroField)
1526   {
1527     fillStatesWithLineFit(ckey, cluster, geometry);
1528   }
1529   else
1530   {
1531     fillStatesWithCircleFit(ckey, cluster, clusglob, geometry);
1532   }
1533
1534   //! skip filling the state information if a state is not there
1535   //! or we just ran the seeding. Fill with Nans to maintain the
1536   //! 1-to-1 mapping between cluster/state vectors
1537   m_statepx.push_back(std::numeric_limits<float>::quiet_NaN());
1538   m_statepy.push_back(std::numeric_limits<float>::quiet_NaN());
1539   m_statepz.push_back(std::numeric_limits<float>::quiet_NaN());
1540   m_statepl.push_back(std::numeric_limits<float>::quiet_NaN());
1541   return;
1542 }
1543
1544 void TrackResiduals::fillStatesWithCircleFit(const TrkrDefs::cluskey& key,
1545                                              TrkrCluster* cluster, Acts::Vector3& glob, ActsGeometry* geometry)
1546 {
1547   auto surf = geometry->maps().getSurface(key, cluster);
1548   std::vector<float> fitpars;
1549   fitpars.push_back(m_R);
1550   fitpars.push_back(m_X0);
1551   fitpars.push_back(m_Y0);
1552   fitpars.push_back(m_rzslope);
1553   fitpars.push_back(m_rzint);
1554
1555   auto intersection = TrackFitUtils::get_helix_surface_intersection(surf, fitpars, glob, geometry);
1556   m_stategx.push_back(intersection.x());
1557   m_stategy.push_back(intersection.y());
1558   m_stategz.push_back(intersection.z());
1559
1560   auto result = surf->globalToLocal(geometry->geometry().getGeoContext(), intersection, Acts::Vector3(1, 1, 1));
1561   if (result.ok())
1562   {
1563     auto loc = result.value() / Acts::UnitConstants::cm;
1564     m_statelx.push_back(loc.x());
1565     m_statelz.push_back(loc.y());
1566   }
1567   else
1568   {
1569     auto local = (surf->transform(geometry->geometry().getGeoContext())).inverse() * (intersection * Acts::UnitConstants::cm);
1570     local /= Acts::UnitConstants::cm;
1571     m_statelx.push_back(local.x());
1572     m_statelz.push_back(local.y());
1573   }
1574 }
1575 void TrackResiduals::fillStatesWithLineFit(const TrkrDefs::cluskey& key,
1576                                            TrkrCluster* cluster, ActsGeometry* geometry)
1577 {
1578   auto intersection = TrackFitUtils::surface_3Dline_intersection(key, cluster, geometry, m_xyslope,
1579                                                                  m_xyint, m_yzslope, m_yzint);
1580
1581   auto surf = geometry->maps().getSurface(key, cluster);
1582   Acts::Vector3 surfnorm = surf->normal(geometry->geometry().getGeoContext(), Acts::Vector3(1, 1, 1), Acts::Vector3(1, 1, 1));
1583   if (!std::isnan(intersection.x()))
1584   {
1585     auto locstateres = surf->globalToLocal(geometry->geometry().getGeoContext(),
1586                                            intersection * Acts::UnitConstants::cm,
1587                                            surfnorm);
1588     if (locstateres.ok())
1589     {
1590       Acts::Vector2 loc = locstateres.value() / Acts::UnitConstants::cm;
1591       m_statelx.push_back(loc(0));
1592       m_statelz.push_back(loc(1));
1593     }
1594     else
1595     {
1596       Acts::Vector3 loct = surf->transform(geometry->geometry().getGeoContext()).inverse() * (intersection * Acts::UnitConstants::cm);
1597       loct /= Acts::UnitConstants::cm;
1598       m_statelx.push_back(loct(0));
1599       m_statelz.push_back(loct(1));
1600     }
1601     m_stategx.push_back(intersection.x());
1602     m_stategy.push_back(intersection.y());
1603     m_stategz.push_back(intersection.z());
1604   }
1605   else
1606   {
1607     //! otherwise the line is parallel to the surface, should not happen if
1608     //! we have a cluster on the surface but just fill the state vecs with nan
1609     m_statelx.push_back(std::numeric_limits<float>::quiet_NaN());
1610     m_statelz.push_back(std::numeric_limits<float>::quiet_NaN());
1611     m_stategx.push_back(std::numeric_limits<float>::quiet_NaN());
1612     m_stategy.push_back(std::numeric_limits<float>::quiet_NaN());
1613     m_stategz.push_back(std::numeric_limits<float>::quiet_NaN());
1614   }
1615 }
1616 void TrackResiduals::createBranches()
1617 {
1618   if (m_doEventTree)
1619   {
1620     m_eventtree = new TTree("eventtree", "A tree with all hits");
1621     m_eventtree->Branch("run", &m_runnumber, "m_runnumber/I");
1622     m_eventtree->Branch("segment", &m_segment, "m_segment/I");
1623     m_eventtree->Branch("event", &m_event, "m_event/I");
1624     m_eventtree->Branch("gl1bco", &m_bco, "m_bco/I");
1625     m_eventtree->Branch("nmvtx", &m_nmvtx_all, "m_nmvtx_all/I");
1626     m_eventtree->Branch("nintt", &m_nintt_all, "m_nintt_all/I");
1627     m_eventtree->Branch("nhittpc0", &m_ntpc_hits0, "m_ntpc_hits0/I");
1628     m_eventtree->Branch("nhittpc1", &m_ntpc_hits1, "m_ntpc_hits1/I");
1629     m_eventtree->Branch("nclustpc0", &m_ntpc_clus0, "m_ntpc_clus0/I");
1630     m_eventtree->Branch("nclustpc1", &m_ntpc_clus1, "m_ntpc_clus1/I");
1631     m_eventtree->Branch("nmms", &m_nmms_all, "m_nmms_all/I");
1632     m_eventtree->Branch("nsiseed", &m_nsiseed, "m_nsiseed/I");
1633     m_eventtree->Branch("ntpcseed", &m_ntpcseed, "m_ntpcseed/I");
1634     m_eventtree->Branch("ntracks", &m_ntracks_all, "m_ntracks_all/I");
1635     m_eventtree->Branch("mbdcharge",&m_totalmbd, "m_totalmbd/F");
1636     m_eventtree->Branch("ntpcClusSector", &m_ntpc_clus_sector);
1637   }
1638
1639   m_failedfits = new TTree("failedfits", "tree with seeds from failed Acts fits");
1640   m_failedfits->Branch("run", &m_runnumber, "m_runnumber/I");
1641   m_failedfits->Branch("segment", &m_segment, "m_segment/I");
1642   m_failedfits->Branch("trackid", &m_trackid, "m_trackid/I");
1643   m_failedfits->Branch("event", &m_event, "m_event/I");
1644   m_failedfits->Branch("silseedx", &m_silseedx, "m_silseedx/F");
1645   m_failedfits->Branch("silseedy", &m_silseedy, "m_silseedy/F");
1646   m_failedfits->Branch("silseedz", &m_silseedz, "m_silseedz/F");
1647   m_failedfits->Branch("tpcseedx", &m_tpcseedx, "m_tpcseedx/F");
1648   m_failedfits->Branch("tpcseedy", &m_tpcseedy, "m_tpcseedy/F");
1649   m_failedfits->Branch("tpcseedz", &m_tpcseedz, "m_tpcseedz/F");
1650   m_failedfits->Branch("tpcseedpx", &m_tpcseedpx, "m_tpcseedpx/F");
1651   m_failedfits->Branch("tpcseedpy", &m_tpcseedpy, "m_tpcseedpy/F");
1652   m_failedfits->Branch("tpcseedpz", &m_tpcseedpz, "m_tpcseedpz/F");
1653   m_failedfits->Branch("tpcseedcharge", &m_tpcseedcharge, "m_tpcseedcharge/I");
1654   m_failedfits->Branch("dedx", &m_dedx, "m_dedx/F");
1655   m_failedfits->Branch("nmaps", &m_nmaps, "m_nmaps/I");
1656   m_failedfits->Branch("nintt", &m_nintt, "m_nintt/I");
1657   m_failedfits->Branch("ntpc", &m_ntpc, "m_ntpc/I");
1658   m_failedfits->Branch("nmms", &m_nmms, "m_nmms/I");
1659   m_failedfits->Branch("gx", &m_clusgx);
1660   m_failedfits->Branch("gy", &m_clusgy);
1661   m_failedfits->Branch("gz", &m_clusgz);
1662   m_failedfits->Branch("gr", &m_clusgr);
1663   m_failedfits->Branch("lx", &m_cluslx);
1664   m_failedfits->Branch("lz", &m_cluslz);
1665
1666   m_vertextree = new TTree("vertextree", "tree with vertices");
1667   m_vertextree->Branch("run", &m_runnumber, "m_runnumber/I");
1668   m_vertextree->Branch("segment", &m_segment, "m_segment/I");
1669   m_vertextree->Branch("event", &m_event, "m_event/I");
1670   m_vertextree->Branch("firedTriggers", &m_firedTriggers);
1671   m_vertextree->Branch("gl1BunchCrossing", &m_gl1BunchCrossing, "m_gl1BunchCrossing/l");
1672   m_vertextree->Branch("gl1bco", &m_bco, "m_bco/l");
1673   m_vertextree->Branch("trbco", &m_bcotr, "m_bcotr/l");
1674   m_vertextree->Branch("vertexid", &m_vertexid);
1675   m_vertextree->Branch("vertex_crossing", &m_vertex_crossing, "m_vertex_crossing/I");
1676   m_vertextree->Branch("mbdzvtx", &m_mbdvtxz, "m_mbdvtxz/F");
1677   m_vertextree->Branch("vx", &m_vx, "m_vx/F");
1678   m_vertextree->Branch("vy", &m_vy, "m_vy/F");
1679   m_vertextree->Branch("vz", &m_vz, "m_vz/F");
1680   m_vertextree->Branch("ntracks", &m_ntracks, "m_ntracks/I");
1681   m_vertextree->Branch("nvertices", &m_nvertices, "m_nvertices/I");
1682   m_vertextree->Branch("gx", &m_clusgx);
1683   m_vertextree->Branch("gy", &m_clusgy);
1684   m_vertextree->Branch("gz", &m_clusgz);
1685   m_vertextree->Branch("gr", &m_clusgr);
1686   m_vertextree->Branch("mbdcharge", &m_totalmbd, "m_totalmbd/F");
1687
1688   m_hittree = new TTree("hittree", "A tree with all hits");
1689   m_hittree->Branch("run", &m_runnumber, "m_runnumber/I");
1690   m_hittree->Branch("segment", &m_segment, "m_segment/I");
1691   m_hittree->Branch("event", &m_event, "m_event/I");
1692   m_hittree->Branch("gl1bco", &m_bco, "m_bco/l");
1693   m_hittree->Branch("hitsetkey", &m_hitsetkey, "m_hitsetkey/i");
1694   m_hittree->Branch("gx", &m_hitgx, "m_hitgx/F");
1695   m_hittree->Branch("gy", &m_hitgy, "m_hitgy/F");
1696   m_hittree->Branch("gz", &m_hitgz, "m_hitgz/F");
1697   m_hittree->Branch("layer", &m_hitlayer, "m_hitlayer/I");
1698   m_hittree->Branch("sector", &m_sector, "m_sector/I");
1699   m_hittree->Branch("side", &m_side, "m_side/I");
1700   m_hittree->Branch("stave", &m_staveid, "m_staveid/I");
1701   m_hittree->Branch("chip", &m_chipid, "m_chipid/I");
1702   m_hittree->Branch("strobe", &m_strobeid, "m_strobeid/I");
1703   m_hittree->Branch("ladderz", &m_ladderzid, "m_ladderzid/I");
1704   m_hittree->Branch("ladderphi", &m_ladderphiid, "m_ladderphiid/I");
1705   m_hittree->Branch("timebucket", &m_timebucket, "m_timebucket/I");
1706   m_hittree->Branch("pad", &m_hitpad, "m_hitpad/I");
1707   m_hittree->Branch("tbin", &m_hittbin, "m_hittbin/I");
1708   m_hittree->Branch("col", &m_col, "m_col/I");
1709   m_hittree->Branch("row", &m_row, "m_row/I");
1710   m_hittree->Branch("segtype", &m_segtype, "m_segtype/I");
1711   m_hittree->Branch("tile", &m_tileid, "m_tileid/I");
1712   m_hittree->Branch("strip", &m_strip, "m_strip/I");
1713   m_hittree->Branch("adc", &m_adc, "m_adc/F");
1714   m_hittree->Branch("zdriftlength", &m_zdriftlength, "m_zdriftlength/F");
1715   m_hittree->Branch("mbdcharge",&m_totalmbd, "m_totalmbd/F");
1716
1717   m_clustree = new TTree("clustertree", "A tree with all clusters");
1718   m_clustree->Branch("run", &m_runnumber, "m_runnumber/I");
1719   m_clustree->Branch("segment", &m_segment, "m_segment/I");
1720   m_clustree->Branch("event", &m_event, "m_event/I");
1721   m_clustree->Branch("gl1bco", &m_bco, "m_bco/l");
1722   m_clustree->Branch("lx", &m_scluslx, "m_scluslx/F");
1723   m_clustree->Branch("lz", &m_scluslz, "m_scluslz/F");
1724   m_clustree->Branch("gx", &m_sclusgx, "m_sclusgx/F");
1725   m_clustree->Branch("gy", &m_sclusgy, "m_sclusgy/F");
1726   m_clustree->Branch("gz", &m_sclusgz, "m_sclusgz/F");
1727   m_clustree->Branch("phi", &m_sclusphi, "m_sclusphi/F");
1728   m_clustree->Branch("eta", &m_scluseta, "m_scluseta/F");
1729   m_clustree->Branch("adc", &m_adc, "m_adc/F");
1730   m_clustree->Branch("phisize", &m_phisize, "m_phisize/I");
1731   m_clustree->Branch("zsize", &m_zsize, "m_zsize/I");
1732   m_clustree->Branch("erphi", &m_scluselx, "m_scluselx/F");
1733   m_clustree->Branch("ez", &m_scluselz, "m_scluselz/F");
1734   m_clustree->Branch("maxadc", &m_clusmaxadc, "m_clusmaxadc/F");
1735   m_clustree->Branch("sector", &m_clussector, "m_clussector/I");
1736   m_clustree->Branch("side", &m_side, "m_side/I");
1737   m_clustree->Branch("stave", &m_staveid, "m_staveid/I");
1738   m_clustree->Branch("chip", &m_chipid, "m_chipid/I");
1739   m_clustree->Branch("strobe", &m_strobeid, "m_strobeid/I");
1740   m_clustree->Branch("ladderz", &m_ladderzid, "m_ladderzid/I");
1741   m_clustree->Branch("ladderphi", &m_ladderphiid, "m_ladderphiid/I");
1742   m_clustree->Branch("timebucket", &m_timebucket, "m_timebucket/I");
1743   m_clustree->Branch("segtype", &m_segtype, "m_segtype/I");
1744   m_clustree->Branch("tile", &m_tileid, "m_tileid/I");
1745
1746   m_tree = new TTree("residualtree", "A tree with track, cluster, and state info");
1747   m_tree->Branch("run", &m_runnumber, "m_runnumber/I");
1748   m_tree->Branch("segment", &m_segment, "m_segment/I");
1749   m_tree->Branch("event", &m_event, "m_event/I");
1750   m_tree->Branch("mbdcharge",&m_totalmbd, "m_totalmbd/F");
1751   m_tree->Branch("mbdzvtx", &m_mbdvtxz, "m_mbdvtxz/F");
1752   m_tree->Branch("firedTriggers", &m_firedTriggers);
1753   m_tree->Branch("gl1BunchCrossing", &m_gl1BunchCrossing, "m_gl1BunchCrossing/l");
1754   m_tree->Branch("trackid", &m_trackid, "m_trackid/I");
1755   m_tree->Branch("tpcid", &m_tpcid, "m_tpcid/I");
1756   m_tree->Branch("silid", &m_silid, "m_silid/I");
1757   m_tree->Branch("gl1bco", &m_bco, "m_bco/l");
1758   m_tree->Branch("crossing", &m_crossing, "m_crossing/I");
1759   m_tree->Branch("crossing_estimate", &m_crossing_estimate, "m_crossing_estimate/I");
1760   m_tree->Branch("silseedx", &m_silseedx, "m_silseedx/F");
1761   m_tree->Branch("silseedy", &m_silseedy, "m_silseedy/F");
1762   m_tree->Branch("silseedz", &m_silseedz, "m_silseedz/F");
1763   m_tree->Branch("silseedpx", &m_silseedpx, "m_silseedpx/F");
1764   m_tree->Branch("silseedpy", &m_silseedpy, "m_silseedpy/F");
1765   m_tree->Branch("silseedpz", &m_silseedpz, "m_silseedpz/F");
1766   m_tree->Branch("silseedphi", &m_silseedphi, "m_silseedphi/F");
1767   m_tree->Branch("silseedeta", &m_silseedeta, "m_silseedeta/F");
1768   m_tree->Branch("silseedcharge", &m_silseedcharge, "m_silseedcharge/I");
1769   m_tree->Branch("tpcseedx", &m_tpcseedx, "m_tpcseedx/F");
1770   m_tree->Branch("tpcseedy", &m_tpcseedy, "m_tpcseedy/F");
1771   m_tree->Branch("tpcseedz", &m_tpcseedz, "m_tpcseedz/F");
1772   m_tree->Branch("tpcseedpx", &m_tpcseedpx, "m_tpcseedpx/F");
1773   m_tree->Branch("tpcseedpy", &m_tpcseedpy, "m_tpcseedpy/F");
1774   m_tree->Branch("tpcseedpz", &m_tpcseedpz, "m_tpcseedpz/F");
1775   m_tree->Branch("tpcseedphi", &m_tpcseedphi, "m_tpcseedphi/F");
1776   m_tree->Branch("tpcseedeta", &m_tpcseedeta, "m_tpcseedeta/F");
1777   m_tree->Branch("tpcseedcharge", &m_tpcseedcharge, "m_tpcseedcharge/I");
1778   m_tree->Branch("dedx", &m_dedx, "m_dedx/F");
1779   m_tree->Branch("tracklength", &m_tracklength, "m_tracklength/F");
1780   m_tree->Branch("px", &m_px, "m_px/F");
1781   m_tree->Branch("py", &m_py, "m_py/F");
1782   m_tree->Branch("pz", &m_pz, "m_pz/F");
1783   m_tree->Branch("pt", &m_pt, "m_pt/F");
1784   m_tree->Branch("eta", &m_eta, "m_eta/F");
1785   m_tree->Branch("phi", &m_phi, "m_phi/F");
1786   m_tree->Branch("deltapt", &m_deltapt, "m_deltapt/F");
1787   m_tree->Branch("charge", &m_charge, "m_charge/I");
1788   m_tree->Branch("quality", &m_quality, "m_quality/F");
1789   m_tree->Branch("ndf", &m_ndf, "m_ndf/F");
1790   m_tree->Branch("nhits", &m_nhits, "m_nhits/I");
1791   m_tree->Branch("nmaps", &m_nmaps, "m_nmaps/I");
1792   m_tree->Branch("nmapsstate", &m_nmapsstate, "m_nmapsstate/I");
1793   m_tree->Branch("nintt", &m_nintt, "m_nintt/I");
1794   m_tree->Branch("ninttstate", &m_ninttstate, "m_ninttstate/I");
1795   m_tree->Branch("ntpc", &m_ntpc, "m_ntpc/I");
1796   m_tree->Branch("ntpcstate", &m_ntpcstate, "m_ntpcstate/I");
1797   m_tree->Branch("nmms", &m_nmms, "m_nmms/I");
1798   m_tree->Branch("nmmsstate", &m_nmmsstate, "m_nmmsstate/I");
1799   m_tree->Branch("tile", &m_tileid, "m_tileid/I");
1800   m_tree->Branch("vertexid", &m_vertexid);
1801   m_tree->Branch("vertex_crossing", &m_vertex_crossing, "m_vertex_crossing/I");
1802   m_tree->Branch("vx", &m_vx, "m_vx/F");
1803   m_tree->Branch("vy", &m_vy, "m_vy/F");
1804   m_tree->Branch("vz", &m_vz, "m_vz/F");
1805   m_tree->Branch("vertex_ntracks", &m_vertex_ntracks, "m_vertex_ntracks/I");
1806   m_tree->Branch("pcax", &m_pcax, "m_pcax/F");
1807   m_tree->Branch("pcay", &m_pcay, "m_pcay/F");
1808   m_tree->Branch("pcaz", &m_pcaz, "m_pcaz/F");
1809   m_tree->Branch("rzslope", &m_rzslope, "m_rzslope/F");
1810   m_tree->Branch("xyslope", &m_xyslope, "m_xyslope/F");
1811   m_tree->Branch("yzslope", &m_yzslope, "m_yzslope/F");
1812   m_tree->Branch("rzint", &m_rzint, "m_rzint/F");
1813   m_tree->Branch("xyint", &m_xyint, "m_xyint/F");
1814   m_tree->Branch("yzint", &m_yzint, "m_yzint/F");
1815   m_tree->Branch("R", &m_R, "m_R/F");
1816   m_tree->Branch("X0", &m_X0, "m_X0/F");
1817   m_tree->Branch("Y0", &m_Y0, "m_Y0/F");
1818   m_tree->Branch("dcaxy", &m_dcaxy, "m_dcaxy/F");
1819   m_tree->Branch("dcaz", &m_dcaz, "m_dcaz/F");
1820
1821   m_tree->Branch("clussector", &m_clsector);
1822   m_tree->Branch("cluslayer", &m_cluslayer);
1823   m_tree->Branch("clusside", &m_clside);
1824   m_tree->Branch("cluskeys", &m_cluskeys);
1825   m_tree->Branch("clusedge", &m_clusedge);
1826   m_tree->Branch("clusoverlap", &m_clusoverlap);
1827   m_tree->Branch("cluslx", &m_cluslx);
1828   m_tree->Branch("cluslz", &m_cluslz);
1829   m_tree->Branch("cluselx", &m_cluselx);
1830   m_tree->Branch("cluselz", &m_cluselz);
1831   m_tree->Branch("clusgx", &m_clusgx);
1832   m_tree->Branch("clusgy", &m_clusgy);
1833   m_tree->Branch("clusgz", &m_clusgz);
1834   m_tree->Branch("clusgr", &m_clusgr);
1835   if (m_doAlignment)
1836   {
1837     m_tree->Branch("clusgxunmoved", &m_clusgxunmoved);
1838     m_tree->Branch("clusgyunmoved", &m_clusgyunmoved);
1839     m_tree->Branch("clusgzunmoved", &m_clusgzunmoved);
1840   }
1841   m_tree->Branch("clusAdc", &m_clusAdc);
1842   m_tree->Branch("clusMaxAdc", &m_clusMaxAdc);
1843   m_tree->Branch("clusphisize", &m_clusphisize);
1844   m_tree->Branch("cluszsize", &m_cluszsize);
1845
1846   if (m_doAlignment)
1847   {
1848     m_tree->Branch("idealsurfcenterx", &m_idealsurfcenterx);
1849     m_tree->Branch("idealsurfcentery", &m_idealsurfcentery);
1850     m_tree->Branch("idealsurfcenterz", &m_idealsurfcenterz);
1851     m_tree->Branch("idealsurfnormx", &m_idealsurfnormx);
1852     m_tree->Branch("idealsurfnormy", &m_idealsurfnormy);
1853     m_tree->Branch("idealsurfnormz", &m_idealsurfnormz);
1854     m_tree->Branch("missurfcenterx", &m_missurfcenterx);
1855     m_tree->Branch("missurfcentery", &m_missurfcentery);
1856     m_tree->Branch("missurfcenterz", &m_missurfcenterz);
1857     m_tree->Branch("missurfnormx", &m_missurfnormx);
1858     m_tree->Branch("missurfnormy", &m_missurfnormy);
1859     m_tree->Branch("missurfnormz", &m_missurfnormz);
1860     m_tree->Branch("clusgxideal", &m_clusgxideal);
1861     m_tree->Branch("clusgyideal", &m_clusgyideal);
1862     m_tree->Branch("clusgzideal", &m_clusgzideal);
1863     m_tree->Branch("missurfalpha", &m_missurfalpha);
1864     m_tree->Branch("missurfbeta", &m_missurfbeta);
1865     m_tree->Branch("missurfgamma", &m_missurfgamma);
1866     m_tree->Branch("idealsurfalpha", &m_idealsurfalpha);
1867     m_tree->Branch("idealsurfbeta", &m_idealsurfbeta);
1868     m_tree->Branch("idealsurfgamma", &m_idealsurfgamma);
1869   }
1870
1871   m_tree->Branch("statelx", &m_statelx);
1872   m_tree->Branch("statelz", &m_statelz);
1873   m_tree->Branch("stateelx", &m_stateelx);
1874   m_tree->Branch("stateelz", &m_stateelz);
1875   m_tree->Branch("stategx", &m_stategx);
1876   m_tree->Branch("stategy", &m_stategy);
1877   m_tree->Branch("stategz", &m_stategz);
1878   m_tree->Branch("statepx", &m_statepx);
1879   m_tree->Branch("statepy", &m_statepy);
1880   m_tree->Branch("statepz", &m_statepz);
1881   m_tree->Branch("statepl", &m_statepl);
1882
1883   if (m_doAlignment)
1884   {
1885     m_tree->Branch("statelxglobderivdx", &m_statelxglobderivdx);
1886     m_tree->Branch("statelxglobderivdy", &m_statelxglobderivdy);
1887     m_tree->Branch("statelxglobderivdz", &m_statelxglobderivdz);
1888     m_tree->Branch("statelxglobderivdalpha", &m_statelxglobderivdalpha);
1889     m_tree->Branch("statelxglobderivdbeta", &m_statelxglobderivdbeta);
1890     m_tree->Branch("statelxglobderivdgamma", &m_statelxglobderivdgamma);
1891
1892     m_tree->Branch("statelxlocderivd0", &m_statelxlocderivd0);
1893     m_tree->Branch("statelxlocderivz0", &m_statelxlocderivz0);
1894     m_tree->Branch("statelxlocderivphi", &m_statelxlocderivphi);
1895     m_tree->Branch("statelxlocderivtheta", &m_statelxlocderivtheta);
1896     m_tree->Branch("statelxlocderivqop", &m_statelxlocderivqop);
1897
1898     m_tree->Branch("statelzglobderivdx", &m_statelzglobderivdx);
1899     m_tree->Branch("statelzglobderivdy", &m_statelzglobderivdy);
1900     m_tree->Branch("statelzglobderivdz", &m_statelzglobderivdz);
1901     m_tree->Branch("statelzglobderivdalpha", &m_statelzglobderivdalpha);
1902     m_tree->Branch("statelzglobderivdbeta", &m_statelzglobderivdbeta);
1903     m_tree->Branch("statelzglobderivdgamma", &m_statelzglobderivdgamma);
1904
1905     m_tree->Branch("statelzlocderivd0", &m_statelzlocderivd0);
1906     m_tree->Branch("statelzlocderivz0", &m_statelzlocderivz0);
1907     m_tree->Branch("statelzlocderivphi", &m_statelzlocderivphi);
1908     m_tree->Branch("statelzlocderivtheta", &m_statelzlocderivtheta);
1909     m_tree->Branch("statelzlocderivqop", &m_statelzlocderivqop);
1910   }
1911 }
1912
1913 void TrackResiduals::fillResidualTreeKF(PHCompositeNode* topNode)
1914 {
1915   auto *silseedmap = findNode::getClass<TrackSeedContainer>(topNode, "SiliconTrackSeedContainer");
1916   auto *tpcseedmap = findNode::getClass<TrackSeedContainer>(topNode, "TpcTrackSeedContainer");
1917   auto *tpcGeom =
1918       findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
1919   auto *trackmap = findNode::getClass<SvtxTrackMap>(topNode, m_trackMapName);
1920   auto *clustermap = findNode::getClass<TrkrClusterContainer>(topNode, m_clusterContainerName);
1921   auto *vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
1922   auto *alignmentmap = findNode::getClass<SvtxAlignmentStateMap>(topNode, m_alignmentMapName);
1923
1924   std::set<unsigned int> tpc_seed_ids;
1925   for (const auto& [key, track] : *trackmap)
1926   {
1927     if (!track)
1928     {
1929       continue;
1930     }
1931     m_trackid = track->get_id();
1932
1933     m_crossing = track->get_crossing();
1934     m_crossing_estimate = SHRT_MAX;
1935     m_px = track->get_px();
1936     m_py = track->get_py();
1937     m_pz = track->get_pz();
1938
1939     m_pt = std::sqrt(square(m_px) + square(m_py));
1940     m_eta = std::atanh(m_pz / std::sqrt(square(m_pt) + square(m_pz)));
1941     m_phi = std::atan2(m_py, m_px);
1942     float CVxx = track->get_error(3, 3);
1943     float CVxy = track->get_error(3, 4);
1944     float CVyy = track->get_error(4, 4);
1945     m_deltapt = std::sqrt((CVxx * square(m_px) + 2 * CVxy * m_px * m_py + CVyy * square(m_py)) / (square(m_px) + square(m_py)));
1946
1947     m_charge = track->get_charge();
1948     m_quality = track->get_quality();
1949     m_chisq = track->get_chisq();
1950     m_ndf = track->get_ndf();
1951
1952     m_nmaps = 0;
1953     m_nmapsstate = 0;
1954     m_nintt = 0;
1955     m_ninttstate = 0;
1956     m_ntpc = 0;
1957     m_ntpcstate = 0;
1958     m_nmms = 0;
1959     m_nmmsstate = 0;
1960     m_silid = std::numeric_limits<unsigned int>::quiet_NaN();
1961     m_tpcid = std::numeric_limits<unsigned int>::quiet_NaN();
1962     m_silseedx = std::numeric_limits<float>::quiet_NaN();
1963     m_silseedy = std::numeric_limits<float>::quiet_NaN();
1964     m_silseedz = std::numeric_limits<float>::quiet_NaN();
1965     m_silseedpx = std::numeric_limits<float>::quiet_NaN();
1966     m_silseedpy = std::numeric_limits<float>::quiet_NaN();
1967     m_silseedpz = std::numeric_limits<float>::quiet_NaN();
1968     m_silseedphi = std::numeric_limits<float>::quiet_NaN();
1969     m_silseedeta = std::numeric_limits<float>::quiet_NaN();
1970     m_silseedcharge = std::numeric_limits<int>::quiet_NaN();
1971     m_tpcseedx = std::numeric_limits<float>::quiet_NaN();
1972     m_tpcseedy = std::numeric_limits<float>::quiet_NaN();
1973     m_tpcseedz = std::numeric_limits<float>::quiet_NaN();
1974     m_tpcseedpx = std::numeric_limits<float>::quiet_NaN();
1975     m_tpcseedpy = std::numeric_limits<float>::quiet_NaN();
1976     m_tpcseedpz = std::numeric_limits<float>::quiet_NaN();
1977     m_tpcseedphi = std::numeric_limits<float>::quiet_NaN();
1978     m_tpcseedeta = std::numeric_limits<float>::quiet_NaN();
1979     m_tpcseedcharge = std::numeric_limits<int>::quiet_NaN();
1980
1981     m_pcax = track->get_x();
1982     m_pcay = track->get_y();
1983     m_pcaz = track->get_z();
1984     m_dcaxy = std::numeric_limits<float>::quiet_NaN();
1985     m_dcaz = std::numeric_limits<float>::quiet_NaN();
1986     m_vertexid = track->get_vertex_id();
1987     if (vertexmap)
1988     {
1989       auto vertexit = vertexmap->find(m_vertexid);
1990       if (vertexit != vertexmap->end())
1991       {
1992         auto *vertex = vertexit->second;
1993         m_vx = vertex->get_x();
1994         m_vy = vertex->get_y();
1995         m_vz = vertex->get_z();
1996         m_vertex_ntracks = vertex->size_tracks();
1997         Acts::Vector3 v(m_vx, m_vy, m_vz);
1998         auto dcapair = TrackAnalysisUtils::get_dca(track, v);
1999         m_dcaxy = dcapair.first.first;
2000         m_dcaz = dcapair.second.first;
2001       }
2002     }
2003
2004     auto *tpcseed = track->get_tpc_seed();
2005     if (tpcseed)
2006     {
2007       m_tpcid = tpcseedmap->find(tpcseed);
2008       tpc_seed_ids.insert(tpcseedmap->find(tpcseed));
2009     }
2010     auto *silseed = track->get_silicon_seed();
2011     if (silseed)
2012     {
2013       m_silid = silseedmap->find(silseed);
2014
2015       const auto si_pos = TrackSeedHelper::get_xyz(silseed);
2016       m_silseedx = si_pos.x();
2017       m_silseedy = si_pos.y();
2018       m_silseedz = si_pos.z();
2019       m_silseedpx = silseed->get_px();
2020       m_silseedpy = silseed->get_py();
2021       m_silseedpz = silseed->get_pz();
2022       m_silseedphi = silseed->get_phi();
2023       m_silseedeta = silseed->get_eta();
2024       m_silseedcharge = silseed->get_qOverR() > 0 ? 1 : -1;
2025     }
2026     if (tpcseed)
2027     {
2028       const auto tpc_pos = TrackSeedHelper::get_xyz(tpcseed);
2029       m_tpcseedx = tpc_pos.x();
2030       m_tpcseedy = tpc_pos.y();
2031       m_tpcseedz = tpc_pos.z();
2032       m_tpcseedpx = tpcseed->get_px();
2033       m_tpcseedpy = tpcseed->get_py();
2034       m_tpcseedpz = tpcseed->get_pz();
2035       m_tpcseedphi = tpcseed->get_phi();
2036       m_tpcseedeta = tpcseed->get_eta();
2037       m_tpcseedcharge = tpcseed->get_qOverR() > 0 ? 1 : -1;
2038     }
2039     if (tpcseed)
2040     {
2041       m_dedx = calc_dedx(tpcseed, clustermap, tpcGeom);
2042     }
2043
2044     if (tpcseed)
2045     {
2046       m_tpcseedpx = tpcseed->get_px();
2047       m_tpcseedpy = tpcseed->get_py();
2048       m_tpcseedpz = tpcseed->get_pz();
2049     }
2050     else if (silseed)
2051     {
2052       m_silseedpx = silseed->get_px();
2053       m_silseedpy = silseed->get_py();
2054       m_silseedpz = silseed->get_pz();
2055     }
2056
2057     clearClusterStateVectors();
2058     if (Verbosity() > 1)
2059     {
2060       std::cout << "Track " << key << " has cluster/states"
2061                 << std::endl;
2062     }
2063
2064     // get the fully corrected cluster global positions
2065     std::vector<std::pair<TrkrDefs::cluskey, Acts::Vector3>> global_raw;
2066     for (const auto& ckey : get_cluster_keys(track))
2067     {
2068       auto *cluster = clustermap->findCluster(ckey);
2069
2070       // Fully correct the cluster positions for the crossing and all distortions
2071       Acts::Vector3 global = m_globalPositionWrapper.getGlobalPositionDistortionCorrected(ckey, cluster, m_crossing);
2072
2073       // add the global positions to a vector to give to the cluster mover
2074       global_raw.emplace_back(ckey, global);
2075     }
2076
2077     // move the cluster positions back to the original readout surface in the fillClusterBranchesKF method
2078
2079     if (!m_doAlignment)
2080     {
2081       for (const auto& ckey : get_cluster_keys(track))
2082       {
2083         fillClusterBranchesKF(ckey, track, global_raw, topNode);
2084       }
2085     }
2086
2087     m_nhits = m_nmaps + m_nintt + m_ntpc + m_nmms;
2088
2089     if (m_doAlignment)
2090     {
2091       /// repopulate with info that is going into alignment
2092       clearClusterStateVectors();
2093
2094       if (alignmentmap and alignmentmap->find(key) != alignmentmap->end())
2095       {
2096         auto& statevec = alignmentmap->find(key)->second;
2097
2098         for (auto& state : statevec)
2099         {
2100           auto ckey = state->get_cluster_key();
2101
2102           fillClusterBranchesKF(ckey, track, global_raw, topNode);
2103
2104           const auto& globderivs = state->get_global_derivative_matrix();
2105           const auto& locderivs = state->get_local_derivative_matrix();
2106
2107           m_statelxglobderivdalpha.push_back(globderivs(0, 0));
2108           m_statelxglobderivdbeta.push_back(globderivs(0, 1));
2109           m_statelxglobderivdgamma.push_back(globderivs(0, 2));
2110           m_statelxglobderivdx.push_back(globderivs(0, 3));
2111           m_statelxglobderivdy.push_back(globderivs(0, 4));
2112           m_statelxglobderivdz.push_back(globderivs(0, 5));
2113
2114           m_statelzglobderivdalpha.push_back(globderivs(1, 0));
2115           m_statelzglobderivdbeta.push_back(globderivs(1, 1));
2116           m_statelzglobderivdgamma.push_back(globderivs(1, 2));
2117           m_statelzglobderivdx.push_back(globderivs(1, 3));
2118           m_statelzglobderivdy.push_back(globderivs(1, 4));
2119           m_statelzglobderivdz.push_back(globderivs(1, 5));
2120
2121           m_statelxlocderivd0.push_back(locderivs(0, 0));
2122           m_statelxlocderivz0.push_back(locderivs(0, 1));
2123           m_statelxlocderivphi.push_back(locderivs(0, 2));
2124           m_statelxlocderivtheta.push_back(locderivs(0, 3));
2125           m_statelxlocderivqop.push_back(locderivs(0, 4));
2126
2127           m_statelzlocderivd0.push_back(locderivs(1, 0));
2128           m_statelzlocderivz0.push_back(locderivs(1, 1));
2129           m_statelzlocderivphi.push_back(locderivs(1, 2));
2130           m_statelzlocderivtheta.push_back(locderivs(1, 3));
2131           m_statelzlocderivqop.push_back(locderivs(1, 4));
2132         }
2133       }
2134     }
2135
2136     if (m_nmms > 0 || !m_doMicromegasOnly)
2137     {
2138       m_tree->Fill();
2139     }
2140
2141   }  // end loop over tracks
2142
2143   if (m_doFailedSeeds)
2144   {
2145     fillFailedSeedTree(topNode, tpc_seed_ids);
2146   }
2147 }
2148 void TrackResiduals::fillEventTree(PHCompositeNode* topNode)
2149 {
2150   auto *silseedmap = findNode::getClass<TrackSeedContainer>(topNode, "SiliconTrackSeedContainer");
2151   auto *tpcseedmap = findNode::getClass<TrackSeedContainer>(topNode, "TpcTrackSeedContainer");
2152   auto *trackmap = findNode::getClass<SvtxTrackMap>(topNode, m_trackMapName);
2153   auto *clustermap = findNode::getClass<TrkrClusterContainer>(topNode, m_clusterContainerName);
2154   auto *hitmap = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
2155
2156   m_ntpc_hits0 = 0;
2157   m_ntpc_hits1 = 0;
2158   m_ntpc_clus0 = 0;
2159   m_ntpc_clus1 = 0;
2160   m_nmvtx_all = 0;
2161   m_nintt_all = 0;
2162   m_nmms_all = 0;
2163   m_nsiseed = 0;
2164   m_ntpcseed = 0;
2165   m_ntpc_clus_sector.resize(24, 0);
2166   m_nsiseed = silseedmap->size();
2167   m_ntpcseed = tpcseedmap->size();
2168   m_ntracks_all = trackmap->size();
2169
2170   // Hits
2171   if (m_doHits)
2172   {
2173     TrkrHitSetContainer::ConstRange tpc_hitsets = hitmap->getHitSets(TrkrDefs::TrkrId::tpcId);
2174     for (TrkrHitSetContainer::ConstIterator hitsetiter = tpc_hitsets.first;
2175          hitsetiter != tpc_hitsets.second;
2176          ++hitsetiter)
2177     {
2178       TrkrDefs::hitsetkey hitsetkey = hitsetiter->first;
2179       TrkrHitSet* hitset = hitsetiter->second;
2180
2181       int thisside = TpcDefs::getSide(hitsetkey);
2182       if (thisside == 0)
2183       {
2184         m_ntpc_hits0 += hitset->size();
2185       }
2186       if (thisside == 1)
2187       {
2188         m_ntpc_hits1 += hitset->size();
2189       }
2190     }
2191   }
2192   for (const auto& det : {TrkrDefs::TrkrId::mvtxId, TrkrDefs::TrkrId::inttId,
2193                     TrkrDefs::TrkrId::tpcId, TrkrDefs::TrkrId::micromegasId})
2194   {
2195     for (const auto& hitsetkey : clustermap->getHitSetKeys(det))
2196     {
2197       auto range = clustermap->getClusters(hitsetkey);
2198       int nclus = std::distance(range.first, range.second);
2199       int tpcside = TrkrDefs::getZElement(hitsetkey);
2200       int sector = TpcDefs::getSectorId(hitsetkey);
2201
2202       switch (det)
2203       {
2204       case TrkrDefs::TrkrId::mvtxId:
2205         m_nmvtx_all += nclus;
2206         break;
2207       case TrkrDefs::TrkrId::inttId:
2208         m_nintt_all += nclus;
2209         break;
2210       case TrkrDefs::TrkrId::tpcId:
2211         if(tpcside == 1)
2212         {
2213           sector += 12;
2214         }
2215         m_ntpc_clus_sector[sector] += nclus;
2216         if (tpcside == 0)
2217         {
2218           m_ntpc_clus0 += nclus;
2219         }
2220         if (tpcside == 1)
2221         {
2222           m_ntpc_clus1 += nclus;
2223         }
2224         break;
2225       case TrkrDefs::TrkrId::micromegasId:
2226         m_nmms_all += nclus;
2227         break;
2228       default:
2229         break;
2230       }
2231     }
2232   }
2233   if (m_doEventTree)
2234   {
2235     if (Verbosity() > 1)
2236     {
2237       std::cout << " m_event:" << m_event << std::endl;
2238       std::cout << " m_ntpc_clus0:" << m_ntpc_clus0 << std::endl;
2239       std::cout << " m_ntpc_clus1: " << m_ntpc_clus1 << std::endl;
2240       std::cout << " m_nmvtx_all:" << m_nmvtx_all << std::endl;
2241       std::cout << " m_nintt_all: " << m_nintt_all << std::endl;
2242       std::cout << " m_nmms_all: " << m_nmms_all << std::endl;
2243     }
2244     m_eventtree->Fill();
2245   }
2246 }
2247
2248 void TrackResiduals::fillResidualTreeSeeds(PHCompositeNode* topNode)
2249 {
2250   auto *silseedmap = findNode::getClass<TrackSeedContainer>(topNode, "SiliconTrackSeedContainer");
2251   auto *tpcseedmap = findNode::getClass<TrackSeedContainer>(topNode, "TpcTrackSeedContainer");
2252   auto *tpcGeom =
2253       findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
2254   auto *trackmap = findNode::getClass<SvtxTrackMap>(topNode, m_trackMapName);
2255   auto *clustermap = findNode::getClass<TrkrClusterContainer>(topNode, m_clusterContainerName);
2256   auto *vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
2257   auto *alignmentmap = findNode::getClass<SvtxAlignmentStateMap>(topNode, m_alignmentMapName);
2258
2259   std::set<unsigned int> tpc_seed_ids;
2260
2261   for (const auto& [key, track] : *trackmap)
2262   {
2263     if (!track)
2264     {
2265       continue;
2266     }
2267     m_trackid = track->get_id();
2268
2269     m_crossing = track->get_crossing();
2270     m_crossing_estimate = SHRT_MAX;
2271     m_px = track->get_px();
2272     m_py = track->get_py();
2273     m_pz = track->get_pz();
2274
2275     m_pt = std::sqrt(square(m_px) + square(m_py));
2276     m_eta = std::atanh(m_pz / std::sqrt(square(m_pt) + square(m_pz)));
2277     m_phi = std::atan2(m_py, m_px);
2278     float CVxx = track->get_error(3, 3);
2279     float CVxy = track->get_error(3, 4);
2280     float CVyy = track->get_error(4, 4);
2281     m_deltapt = std::sqrt((CVxx * square(m_px) + 2 * CVxy * m_px * m_py + CVyy * square(m_py)) / (square(m_px) + square(m_py)));
2282
2283     m_charge = track->get_charge();
2284     m_quality = track->get_quality();
2285     m_chisq = track->get_chisq();
2286     m_ndf = track->get_ndf();
2287
2288     if (Verbosity() > 1)
2289     {
2290       std::cout << "fillResidualTreeSeeds:  track " << m_trackid << " m_crossing " << m_crossing << " m_crossing_estimate " << m_crossing_estimate << " m_pt " << m_pt << std::endl;
2291     }
2292
2293     m_nmaps = 0;
2294     m_nintt = 0;
2295     m_ntpc = 0;
2296     m_nmms = 0;
2297     m_silid = std::numeric_limits<unsigned int>::quiet_NaN();
2298     m_tpcid = std::numeric_limits<unsigned int>::quiet_NaN();
2299     m_silseedx = std::numeric_limits<float>::quiet_NaN();
2300     m_silseedy = std::numeric_limits<float>::quiet_NaN();
2301     m_silseedz = std::numeric_limits<float>::quiet_NaN();
2302     m_silseedpx = std::numeric_limits<float>::quiet_NaN();
2303     m_silseedpy = std::numeric_limits<float>::quiet_NaN();
2304     m_silseedpz = std::numeric_limits<float>::quiet_NaN();
2305     m_silseedphi = std::numeric_limits<float>::quiet_NaN();
2306     m_silseedeta = std::numeric_limits<float>::quiet_NaN();
2307     m_silseedcharge = std::numeric_limits<int>::quiet_NaN();
2308     m_tpcseedx = std::numeric_limits<float>::quiet_NaN();
2309     m_tpcseedy = std::numeric_limits<float>::quiet_NaN();
2310     m_tpcseedz = std::numeric_limits<float>::quiet_NaN();
2311     m_tpcseedpx = std::numeric_limits<float>::quiet_NaN();
2312     m_tpcseedpy = std::numeric_limits<float>::quiet_NaN();
2313     m_tpcseedpz = std::numeric_limits<float>::quiet_NaN();
2314     m_tpcseedphi = std::numeric_limits<float>::quiet_NaN();
2315     m_tpcseedeta = std::numeric_limits<float>::quiet_NaN();
2316     m_tpcseedcharge = std::numeric_limits<int>::quiet_NaN();
2317
2318     m_vertexid = track->get_vertex_id();
2319
2320     m_pcax = track->get_x();
2321     m_pcay = track->get_y();
2322     m_pcaz = track->get_z();
2323     m_dcaxy = std::numeric_limits<float>::quiet_NaN();
2324     m_dcaz = std::numeric_limits<float>::quiet_NaN();
2325     if (vertexmap)
2326     {
2327       auto vertexit = vertexmap->find(m_vertexid);
2328       if (vertexit != vertexmap->end())
2329       {
2330         auto *vertex = vertexit->second;
2331         m_vx = vertex->get_x();
2332         m_vy = vertex->get_y();
2333         m_vz = vertex->get_z();
2334         Acts::Vector3 vert(m_vx, m_vy, m_vz);
2335         auto dcapair = TrackAnalysisUtils::get_dca(track, vert);
2336         m_dcaxy = dcapair.first.first;
2337         m_dcaz = dcapair.second.first;
2338       }
2339     }
2340
2341     auto *tpcseed = track->get_tpc_seed();
2342     if (tpcseed)
2343     {
2344       m_tpcid = tpcseedmap->find(tpcseed);
2345       tpc_seed_ids.insert(tpcseedmap->find(tpcseed));
2346     }
2347     auto *silseed = track->get_silicon_seed();
2348     if (silseed)
2349     {
2350       m_silid = silseedmap->find(silseed);
2351       const auto si_pos = TrackSeedHelper::get_xyz(silseed);
2352       m_silseedx = si_pos.x();
2353       m_silseedy = si_pos.y();
2354       m_silseedz = si_pos.z();
2355       m_silseedpx = silseed->get_px();
2356       m_silseedpy = silseed->get_py();
2357       m_silseedpz = silseed->get_pz();
2358       m_silseedphi = silseed->get_phi();
2359       m_silseedeta = silseed->get_eta();
2360       m_silseedcharge = silseed->get_qOverR() > 0 ? 1 : -1;
2361     }
2362     if (tpcseed)
2363     {
2364       const auto tpc_pos = TrackSeedHelper::get_xyz(tpcseed);
2365       m_tpcseedx = tpc_pos.x();
2366       m_tpcseedy = tpc_pos.y();
2367       m_tpcseedz = tpc_pos.z();
2368       m_tpcseedpx = tpcseed->get_px();
2369       m_tpcseedpy = tpcseed->get_py();
2370       m_tpcseedpz = tpcseed->get_pz();
2371       m_tpcseedphi = tpcseed->get_phi();
2372       m_tpcseedeta = tpcseed->get_eta();
2373       m_tpcseedcharge = tpcseed->get_qOverR() > 0 ? 1 : -1;
2374     }
2375     if (tpcseed)
2376     {
2377       m_dedx = calc_dedx(tpcseed, clustermap, tpcGeom);
2378     }
2379     if (m_zeroField)
2380     {
2381       float qor = std::numeric_limits<float>::quiet_NaN();
2382       float phi = std::numeric_limits<float>::quiet_NaN();
2383       float theta = std::numeric_limits<float>::quiet_NaN();
2384       float eta = std::numeric_limits<float>::quiet_NaN();
2385       if (tpcseed)
2386       {
2387         qor = tpcseed->get_qOverR();
2388         phi = tpcseed->get_phi();
2389         eta = tpcseed->get_eta();
2390         theta = tpcseed->get_theta();
2391       }
2392       else if (silseed)
2393       {
2394         qor = silseed->get_qOverR();
2395         phi = silseed->get_phi();
2396         eta = silseed->get_eta();
2397         theta = silseed->get_theta();
2398       }
2399       float pt = fabs(1. / qor) * (0.3 / 100) * 0.01;
2400       m_tpcseedpx = pt * std::cos(phi);
2401       m_tpcseedpy = pt * std::sin(phi);
2402       m_tpcseedpz = pt * std::cosh(eta) * std::cos(theta);
2403     }
2404     else
2405     {
2406       if (tpcseed)
2407       {
2408         m_tpcseedpx = tpcseed->get_px();
2409         m_tpcseedpy = tpcseed->get_py();
2410         m_tpcseedpz = tpcseed->get_pz();
2411       }
2412       else if (silseed)
2413       {
2414         m_silseedpx = silseed->get_px();
2415         m_silseedpy = silseed->get_py();
2416         m_silseedpz = silseed->get_pz();
2417       }
2418     }
2419     clearClusterStateVectors();
2420     if (Verbosity() > 1)
2421     {
2422       std::cout << "Track " << key << " has cluster/states"
2423                 << std::endl;
2424     }
2425
2426     // get the fully corrected cluster global positions
2427     std::vector<std::pair<TrkrDefs::cluskey, Acts::Vector3>> global_raw;
2428     float minR = std::numeric_limits<float>::max();
2429     float maxR = 0;
2430     for (const auto& ckey : get_cluster_keys(track))
2431     {
2432       auto *cluster = clustermap->findCluster(ckey);
2433
2434       // Fully correct the cluster positions for the crossing and all distortions
2435       Acts::Vector3 global = m_globalPositionWrapper.getGlobalPositionDistortionCorrected(ckey, cluster, m_crossing);
2436
2437       // add the global positions to a vector to give to the cluster mover
2438       global_raw.emplace_back(ckey, global);
2439       minR = std::min<double>(r(global.x(), global.y()), minR);
2440       maxR = std::max<double>(r(global.x(), global.y()), maxR);
2441     }
2442     m_tracklength = maxR - minR;
2443     // ---- we move the global positions back to the surface in fillClusterBranchesSeeds
2444
2445     if (!m_doAlignment)
2446     {
2447       std::vector<TrkrDefs::cluskey> keys;
2448       for (const auto& ckey : get_cluster_keys(track))
2449       {
2450         keys.push_back(ckey);
2451       }
2452       if (m_zeroField)
2453       {
2454         lineFitClusters(keys, clustermap, m_crossing);
2455       }
2456       else
2457       {
2458         // this corrects the cluster positions and fits them, to fill the helical fit parameters
2459         //  that are used to calculate the "state" positions
2460         circleFitClusters(keys, clustermap, m_crossing);
2461       }
2462
2463       for (const auto& ckey : get_cluster_keys(track))
2464       {
2465         fillClusterBranchesSeeds(ckey, global_raw, topNode);
2466       }
2467     }
2468
2469     m_nhits = m_nmaps + m_nintt + m_ntpc + m_nmms;
2470
2471     if (m_doAlignment)
2472     {
2473       /// repopulate with info that is going into alignment
2474       clearClusterStateVectors();
2475
2476       if (alignmentmap and alignmentmap->find(key) != alignmentmap->end())
2477       {
2478         auto& statevec = alignmentmap->find(key)->second;
2479
2480         for (auto& state : statevec)
2481         {
2482           auto ckey = state->get_cluster_key();
2483
2484           fillClusterBranchesSeeds(ckey, global_raw, topNode);
2485
2486           const auto& globderivs = state->get_global_derivative_matrix();
2487           const auto& locderivs = state->get_local_derivative_matrix();
2488
2489           m_statelxglobderivdalpha.push_back(globderivs(0, 0));
2490           m_statelxglobderivdbeta.push_back(globderivs(0, 1));
2491           m_statelxglobderivdgamma.push_back(globderivs(0, 2));
2492           m_statelxglobderivdx.push_back(globderivs(0, 3));
2493           m_statelxglobderivdy.push_back(globderivs(0, 4));
2494           m_statelxglobderivdz.push_back(globderivs(0, 5));
2495
2496           m_statelzglobderivdalpha.push_back(globderivs(1, 0));
2497           m_statelzglobderivdbeta.push_back(globderivs(1, 1));
2498           m_statelzglobderivdgamma.push_back(globderivs(1, 2));
2499           m_statelzglobderivdx.push_back(globderivs(1, 3));
2500           m_statelzglobderivdy.push_back(globderivs(1, 4));
2501           m_statelzglobderivdz.push_back(globderivs(1, 5));
2502
2503           m_statelxlocderivd0.push_back(locderivs(0, 0));
2504           m_statelxlocderivz0.push_back(locderivs(0, 1));
2505           m_statelxlocderivphi.push_back(locderivs(0, 2));
2506           m_statelxlocderivtheta.push_back(locderivs(0, 3));
2507           m_statelxlocderivqop.push_back(locderivs(0, 4));
2508
2509           m_statelzlocderivd0.push_back(locderivs(1, 0));
2510           m_statelzlocderivz0.push_back(locderivs(1, 1));
2511           m_statelzlocderivphi.push_back(locderivs(1, 2));
2512           m_statelzlocderivtheta.push_back(locderivs(1, 3));
2513           m_statelzlocderivqop.push_back(locderivs(1, 4));
2514         }
2515       }
2516     }
2517     if (!m_doMatchedOnly)
2518     {
2519       m_tree->Fill();
2520     }
2521     else
2522     {
2523       if (m_nmaps >= 3 && m_nintt >= 1 && m_ntpc > 32 && abs(m_crossing) < 5 && m_pt > 0.3)
2524       {
2525         m_tree->Fill();
2526       }
2527     }
2528   }
2529 }