packages/trackbase_historic/TrackAnalysisUtils.cc

0001 #include "TrackAnalysisUtils.h"
0002
0003 #include <trackbase/ActsGeometry.h>
0004 #include <trackbase/TpcDefs.h>
0005 #include <trackbase/TrkrCluster.h>
0006 #include <trackbase/TrkrClusterContainer.h>
0007 #include <globalvertex/GlobalVertex.h>
0008
0009 #include <tpc/TpcClusterZCrossingCorrection.h>
0010 #include "SvtxTrack.h"
0011 #include "TrackSeed.h"
0012
0013 #include <cmath>
0014
0015 namespace TrackAnalysisUtils
0016 {
0017
0018   float calc_dedx(TrackSeed* tpcseed,
0019                   TrkrClusterContainer* clustermap,
0020                   ActsGeometry* tgeometry,
0021                   float thickness_per_region[4])
0022   {
0023     std::vector<TrkrDefs::cluskey> clusterKeys;
0024     clusterKeys.insert(clusterKeys.end(), tpcseed->begin_cluster_keys(),
0025                        tpcseed->end_cluster_keys());
0026
0027     std::vector<float> dedxlist;
0028     for (unsigned long cluster_key : clusterKeys)
0029     {
0030       auto detid = TrkrDefs::getTrkrId(cluster_key);
0031       if (detid != TrkrDefs::TrkrId::tpcId)
0032       {
0033         continue;  // the micromegas clusters are added to the TPC seeds
0034       }
0035       unsigned int layer_local = TrkrDefs::getLayer(cluster_key);
0036       TrkrCluster* cluster = clustermap->findCluster(cluster_key);
0037       float adc = cluster->getAdc();
0038       float thick = 0;
0039       if (layer_local < 23)
0040       {
0041         if (layer_local % 2 == 0)
0042         {
0043           thick = thickness_per_region[1];
0044         }
0045         else
0046         {
0047           thick = thickness_per_region[0];
0048         }
0049       }
0050       else if (layer_local < 39)
0051       {
0052         thick = thickness_per_region[2];
0053       }
0054       else
0055       {
0056         thick = thickness_per_region[3];
0057       }
0058       auto cglob = tgeometry->getGlobalPosition(cluster_key, cluster);
0059       float r = std::sqrt(cglob(0) * cglob(0) + cglob(1) * cglob(1));
0060       float alpha = (r * r) / (2 * r * std::abs(1.0 / tpcseed->get_qOverR()));
0061       float beta = std::atan(tpcseed->get_slope());
0062       float alphacorr = std::cos(alpha);
0063       if (alphacorr < 0 || alphacorr > 4)
0064       {
0065         alphacorr = 4;
0066       }
0067       float betacorr = std::cos(beta);
0068       if (betacorr < 0 || betacorr > 4)
0069       {
0070         betacorr = 4;
0071       }
0072       adc /= thick;
0073       adc *= alphacorr;
0074       adc *= betacorr;
0075       dedxlist.push_back(adc);
0076       sort(dedxlist.begin(), dedxlist.end());
0077     }
0078     int trunc_min = 0;
0079     if (dedxlist.empty())
0080     {
0081       return std::numeric_limits<float>::quiet_NaN();
0082     }
0083     int trunc_max = (int) dedxlist.size() * 0.7;
0084     float sumdedx = 0;
0085     int ndedx = 0;
0086     for (int j = trunc_min; j <= trunc_max; j++)
0087     {
0088       sumdedx += dedxlist.at(j);
0089       ndedx++;
0090     }
0091     sumdedx /= ndedx;
0092     return sumdedx;
0093   }
0094
0095   float calc_dedx_calib(SvtxTrack* track,
0096                         TrkrClusterContainer* cluster_map,
0097                         ActsGeometry* tgeometry,
0098                         float thickness_per_region[4])
0099   {
0100     auto clusterKeys = get_cluster_keys(track->get_tpc_seed());
0101
0102     std::vector<float> dedxlist;
0103     for (unsigned long cluster_key : clusterKeys)
0104     {
0105       unsigned int layer_local = TrkrDefs::getLayer(cluster_key);
0106       if (TrkrDefs::getTrkrId(cluster_key) != TrkrDefs::TrkrId::tpcId)
0107       {
0108         continue;
0109       }
0110       TrkrCluster* cluster = cluster_map->findCluster(cluster_key);
0111       Acts::Vector3 cglob;
0112       cglob = tgeometry->getGlobalPosition(cluster_key, cluster);
0113
0114       int csegment = 0;
0115       float thickness = 0;
0116       if (layer_local >= 7 + 32)
0117       {
0118         csegment = 2;
0119         thickness = thickness_per_region[3];
0120       }
0121       else if (layer_local >= 7 + 16)
0122       {
0123         csegment = 1;
0124         thickness = thickness_per_region[2];
0125       }
0126       if (csegment == 0)
0127       {
0128         if (layer_local % 2 == 0)
0129         {
0130           thickness = thickness_per_region[1];
0131         }
0132         else
0133         {
0134           thickness = thickness_per_region[0];
0135         }
0136       }
0137       float adc = cluster->getAdc();
0138
0139       float r = std::sqrt(cglob(0) * cglob(0) + cglob(1) * cglob(1));
0140       auto tpcseed = track->get_tpc_seed();
0141       float alpha = (r * r) / (2 * r * std::abs(1.0 / tpcseed->get_qOverR()));
0142       float beta = std::atan(tpcseed->get_slope());
0143       float alphacorr = std::cos(alpha);
0144       if (alphacorr < 0 || alphacorr > 4)
0145       {
0146         alphacorr = 4;
0147       }
0148       float betacorr = std::cos(beta);
0149       if (betacorr < 0 || betacorr > 4)
0150       {
0151         betacorr = 4;
0152       }
0153       if(track->get_crossing() < SHRT_MAX)
0154       {
0155         double z_crossing_corrected =
0156           TpcClusterZCrossingCorrection::correctZ(cglob.z(),
0157           TpcDefs::getSide(cluster_key), track->get_crossing());
0158
0159     double maxz = tgeometry->get_max_driftlength() + tgeometry->get_CM_halfwidth();
0160     adc /= (1 - ((maxz - abs(z_crossing_corrected)) * 0.50 / maxz));
0161       }
0162       adc /= thickness;
0163       adc *= alphacorr;
0164       adc *= betacorr;
0165       dedxlist.push_back(adc);
0166       sort(dedxlist.begin(), dedxlist.end());
0167     }
0168     int trunc_min = 0;
0169     int trunc_max = (int) dedxlist.size() * 0.7;
0170     float sumdedx = 0;
0171     int ndedx = 0;
0172     for (int j = trunc_min; j <= trunc_max; j++)
0173     {
0174       sumdedx += dedxlist.at(j);
0175       ndedx++;
0176     }
0177     sumdedx /= ndedx;
0178     return sumdedx;
0179   }
0180
0181   TrackAnalysisUtils::DCAPair get_dca(SvtxTrack *track,
0182                                       GlobalVertex* vertex)
0183   {
0184     Acts::Vector3 vpos(vertex->get_x(),
0185                        vertex->get_y(),
0186                        vertex->get_z());
0187     Acts::Vector3 mom(track->get_px(),
0188                       track->get_py(),
0189                       track->get_pz());
0190     auto dca = get_dca(track, vpos);
0191     auto rot = rotationMatrixToLocal(mom);
0192     Acts::RotationMatrix3 rot_T = rot.transpose();
0193
0194     Acts::ActsSquareMatrix<3> posCov;
0195     Acts::ActsSquareMatrix<3> vertexCov;
0196     for (int i = 0; i < 3; ++i)
0197     {
0198       for (int j = 0; j < 3; ++j)
0199       {
0200         posCov(i, j) = track->get_error(i, j);
0201         vertexCov(i, j) = vertex->get_error(i, j);
0202       }
0203     }
0204
0205
0206     Acts::ActsSquareMatrix<3> rotCov = rot * (posCov+vertexCov) * rot_T;
0207     dca.first.second = sqrt(rotCov(0, 0));
0208     dca.second.second = sqrt(rotCov(2, 2));
0209
0210     return dca;
0211   }
0212
0213   TrackAnalysisUtils::DCAPair get_dca(SvtxTrack* track,
0214                                       Acts::Vector3& vertex)
0215   {
0216     TrackAnalysisUtils::DCAPair pair;
0217     if (!track)
0218     {
0219       std::cout << "You passed a nullptr, can't calculate DCA" << std::endl;
0220       return pair;
0221     }
0222
0223     Acts::Vector3 pos(track->get_x(),
0224                       track->get_y(),
0225                       track->get_z());
0226     Acts::Vector3 mom(track->get_px(),
0227                       track->get_py(),
0228                       track->get_pz());
0229
0230     pos -= vertex;
0231
0232     Acts::ActsSquareMatrix<3> posCov;
0233     for (int i = 0; i < 3; ++i)
0234     {
0235       for (int j = 0; j < 3; ++j)
0236       {
0237         posCov(i, j) = track->get_error(i, j);
0238       }
0239     }
0240
0241     auto rot = rotationMatrixToLocal(mom);
0242     Acts::RotationMatrix3 rot_T = rot.transpose();
0243
0244     Acts::Vector3 pos_R = rot * pos;
0245     Acts::ActsSquareMatrix<3> rotCov = rot * posCov * rot_T;
0246     //! First pair is DCA_xy +/- DCA_xy_err
0247     pair.first.first = pos_R(0);
0248     pair.first.second = sqrt(rotCov(0, 0));
0249
0250     //! Second pair is DCA_z +/- DCA_z_err
0251     pair.second.first = pos_R(2);
0252     pair.second.second = sqrt(rotCov(2, 2));
0253
0254     return pair;
0255   }
0256   Acts::RotationMatrix3 rotationMatrixToLocal(const Acts::Vector3& mom)
0257   {
0258     Acts::Vector3 r = mom.cross(Acts::Vector3(0., 0., 1.));
0259     float phi = atan2(r(1), r(0));
0260     phi *= -1;
0261     Acts::RotationMatrix3 rot;
0262
0263     rot(0, 0) = std::cos(phi);
0264     rot(0, 1) = -std::sin(phi);
0265     rot(0, 2) = 0;
0266     rot(1, 0) = std::sin(phi);
0267     rot(1, 1) = std::cos(phi);
0268     rot(1, 2) = 0;
0269     rot(2, 0) = 0;
0270     rot(2, 1) = 0;
0271     rot(2, 2) = 1;
0272     return rot;
0273   }
0274   std::vector<TrkrDefs::cluskey> get_cluster_keys(TrackSeed* seed)
0275   {
0276     std::vector<TrkrDefs::cluskey> out;
0277
0278       if (seed)
0279       {
0280         std::copy(seed->begin_cluster_keys(), seed->end_cluster_keys(), std::back_inserter(out));
0281       }
0282
0283     return out;
0284   }
0285
0286   std::vector<TrkrDefs::cluskey> get_cluster_keys(SvtxTrack* track)
0287   {
0288     std::vector<TrkrDefs::cluskey> out;
0289     for (const auto& seed : {track->get_silicon_seed(), track->get_tpc_seed()})
0290     {
0291       if (seed)
0292       {
0293         std::copy(seed->begin_cluster_keys(), seed->end_cluster_keys(), std::back_inserter(out));
0294       }
0295     }
0296     return out;
0297   }
0298
0299 }  // namespace TrackAnalysisUtils