From fedf524c66a021ab37d90dbf7dfa9c69ef80878c Mon Sep 17 00:00:00 2001 From: Matevz Tadel <mtadel@ucsd.edu> Date: Thu, 21 Dec 2017 15:26:52 -0800 Subject: [PATCH 01/11] Canonize mplex propagation and update, make b-field and material flags mandatory. --- Config.h | 39 +++ Matrix.h | 1 + mkFit/KalmanUtilsMPlex.cc | 564 ++++++++++++++++++------------------- mkFit/KalmanUtilsMPlex.h | 89 ++++-- mkFit/MkBase.h | 60 +++- mkFit/MkBuilder.cc | 18 +- mkFit/MkFinder.cc | 38 +-- mkFit/MkFinder.h | 9 +- mkFit/MkFitter.cc | 32 +-- mkFit/PropagationMPlex.cc | 101 +++---- mkFit/PropagationMPlex.h | 37 ++- mkFit/PropagationMPlex.icc | 14 +- mkFit/SteeringParams.h | 8 +- 13 files changed, 558 insertions(+), 452 deletions(-) diff --git a/Config.h b/Config.h index fc94c6ea82d4c..0443b67e9a054 100644 --- a/Config.h +++ b/Config.h @@ -15,6 +15,43 @@ class TrackerInfo; #define CUDA_CALLABLE #endif + +//------------------------------------------------------------------------------ + +enum PropagationFlagsEnum +{ + PF_none = 0, + + PF_use_param_b_field = 0x1, + PF_apply_material = 0x2 +}; + +struct PropagationFlags +{ + union + { + struct + { + bool use_param_b_field : 1; + bool apply_material : 1; + // Could add: bool use_trig_approx -- now Config::useTrigApprox = true + // Could add: int n_iter : 8 -- now Config::Niter = 5 + }; + + unsigned int _raw_; + + }; + + PropagationFlags() : _raw_(0) {} + + PropagationFlags(int pfe) : + use_param_b_field ( pfe & PF_use_param_b_field), + apply_material ( pfe & PF_apply_material) + {} +}; + +//------------------------------------------------------------------------------ + // Enum for input seed options enum seedOpts {simSeeds, cmsswSeeds, findSeeds}; typedef std::map<std::string,seedOpts> seedOptsMap; @@ -27,6 +64,8 @@ typedef std::map<std::string,cleanOpts> cleanOptsMap; enum matchOpts {trkParamBased, hitBased, labelBased}; typedef std::map<std::string, matchOpts> matchOptsMap; +//------------------------------------------------------------------------------ + namespace Config { extern TrackerInfo TrkInfo; diff --git a/Matrix.h b/Matrix.h index 903ceeea0f3cd..6c10bbb9367fa 100644 --- a/Matrix.h +++ b/Matrix.h @@ -75,6 +75,7 @@ inline void sincos4(const float x, float& sin, float& cos) cos = 1.f - 0.5f*x2 + 0.04166667f*x2*x2; sin = x - 0.16666667f*x*x2; } + //============================================================================== // This ifdef needs to be changed to something like "use matriplex" and/or diff --git a/mkFit/KalmanUtilsMPlex.cc b/mkFit/KalmanUtilsMPlex.cc index 32f43679337d4..ced1c8b51e32f 100644 --- a/mkFit/KalmanUtilsMPlex.cc +++ b/mkFit/KalmanUtilsMPlex.cc @@ -437,180 +437,119 @@ void KHC(const MPlexL2& A, const MPlexLS& B, MPlexLS& C) //============================================================================== -// updateParametersMPlex +// Kalman operations - common dummy variables //============================================================================== -void updateParametersMPlex(const MPlexLS &psErr, const MPlexLV& psPar, const MPlexQI &inChg, - const MPlexHS &msErr, const MPlexHV& msPar, - MPlexLS &outErr, MPlexLV& outPar, - const int N_proc) +namespace { - // const idx_t N = psErr.N; - // Assert N-s of all parameters are the same. + // Dummy variables for parameter consistency to kalmanOperation. + // Through KalmanFilterOperation enum parameter it is guaranteed that + // those will never get accessed in the code (read from or written into). - // Temporaries -- this is expensive -- should have them allocated outside and reused. - // Can be passed in in a struct, see above. + MPlexLS dummy_err; + MPlexLV dummy_par; + MPlexQF dummy_chi2; +} - // updateParametersContext ctx; - //assert((long long)(&updateCtx.propErr.fArray[0]) % 64 == 0); - // debug = true; +//============================================================================== +// Kalman operations - Barrel +//============================================================================== +void kalmanUpdate(const MPlexLS &psErr, const MPlexLV& psPar, + const MPlexHS &msErr, const MPlexHV& msPar, + MPlexLS &outErr, MPlexLV& outPar, + const int N_proc) +{ + kalmanOperation(KFO_Update_Params, psErr, psPar, msErr, msPar, + outErr, outPar, dummy_chi2, N_proc); +} + +void kalmanPropagateAndUpdate(const MPlexLS &psErr, const MPlexLV& psPar, const MPlexQI &inChg, + const MPlexHS &msErr, const MPlexHV& msPar, + MPlexLS &outErr, MPlexLV& outPar, + const int N_proc, const PropagationFlags propFlags) +{ MPlexLS propErr; MPlexLV propPar; - // do a full propagation step to correct for residual distance from the hit radius - need the charge for this - if (Config::useCMSGeom) { - propagateHelixToRMPlex(psErr, psPar, inChg, msPar, propErr, propPar, N_proc); - } else { - propErr = psErr; - propPar = psPar; - } -#ifdef DEBUG - { - dmutex_guard; - printf("propPar:\n"); - for (int i = 0; i < 6; ++i) { - printf("%8f ", propPar.ConstAt(0,0,i)); printf("\n"); - } printf("\n"); - printf("msPar:\n"); - for (int i = 0; i < 3; ++i) { - printf("%8f ", msPar.ConstAt(0,0,i)); printf("\n"); - } printf("\n"); - printf("propErr:\n"); - for (int i = 0; i < 6; ++i) { for (int j = 0; j < 6; ++j) - printf("%8f ", propErr.At(0,i,j)); printf("\n"); - } printf("\n"); - printf("msErr:\n"); - for (int i = 0; i < 6; ++i) { for (int j = 0; j < 6; ++j) - printf("%8f ", msErr.ConstAt(0,i,j)); printf("\n"); - } printf("\n"); - } -#endif - - // Rotate global point on tangent plane to cylinder - // Tangent point is half way between hit and propagate position - - // Rotation matrix - // rotT00 0 rotT01 - // rotT01 0 -rotT00 - // 0 1 0 - // Minimize temporaries: only two float are needed! + // XXXX The following if could go away if we use proper update_foos in + // steering and set them from geom. For now it is still here (I'm going + // crazy already as things are). - MPlexQF rotT00; - MPlexQF rotT01; + if (Config::useCMSGeom) { + MPlexQF msRad; #pragma simd - for (int n = 0; n < NN; ++n) { - float r = hipo(msPar.ConstAt(n, 0, 0), msPar.ConstAt(n, 1, 0)); - rotT00.At(n, 0, 0) = -(msPar.ConstAt(n, 1, 0)+propPar.ConstAt(n, 1, 0))/(2*r); - rotT01.At(n, 0, 0) = (msPar.ConstAt(n, 0, 0)+propPar.ConstAt(n, 0, 0))/(2*r); - } - - MPlexHV res_glo; //position residual in global coordinates - SubtractFirst3(msPar, propPar, res_glo); - - MPlexHS resErr_glo;//covariance sum in global position coordinates - AddIntoUpperLeft3x3(propErr, msErr, resErr_glo); - - MPlex2V res_loc; //position residual in local coordinates - RotateResidulsOnTangentPlane(rotT00,rotT01,res_glo,res_loc); - MPlex2S resErr_loc;//covariance sum in local position coordinates - MPlexHH tempHH; - ProjectResErr (rotT00, rotT01, resErr_glo, tempHH); - ProjectResErrTransp(rotT00, rotT01, tempHH, resErr_loc); + for (int n = 0; n < NN; ++n) + { + msRad.At(n, 0, 0) = std::hypot(msPar.ConstAt(n, 0, 0), msPar.ConstAt(n, 1, 0)); + } -#ifdef DEBUG - { - dmutex_guard; - printf("resErr:\n"); - for (int i = 0; i < 2; ++i) { for (int j = 0; j < 2; ++j) - printf("%8f ", resErr_loc.At(0,i,j)); printf("\n"); - } printf("\n"); + propagateHelixToRMPlex(psErr, psPar, inChg, msRad, propErr, propPar, N_proc, propFlags); + } else { + propErr = psErr; + propPar = psPar; } -#endif - - //invert the 2x2 matrix - Matriplex::InvertCramerSym(resErr_loc); - - MPlexLH K; // kalman gain, fixme should be L2 - KalmanHTG(rotT00, rotT01, resErr_loc, tempHH); // intermediate term to get kalman gain (H^T*G) - KalmanGain(propErr, tempHH, K); - - MultResidualsAdd(K, propPar, res_loc, outPar); - MPlexLL tempLL; - squashPhiMPlex(outPar,N_proc); // ensure phi is between |pi| + kalmanOperation(KFO_Update_Params, propErr, propPar, msErr, msPar, + outErr, outPar, dummy_chi2, N_proc); +} - KHMult(K, rotT00, rotT01, tempLL); - KHC(tempLL, propErr, outErr); - outErr.Subtract(propErr, outErr); +//------------------------------------------------------------------------------ -#ifdef DEBUG - { - dmutex_guard; - printf("res_glo:\n"); - for (int i = 0; i < 3; ++i) { - printf("%8f ", res_glo.At(0,i,0)); - } printf("\n"); - printf("res_loc:\n"); - for (int i = 0; i < 2; ++i) { - printf("%8f ", res_loc.At(0,i,0)); - } printf("\n"); - printf("resErr_loc (Inv):\n"); - for (int i = 0; i < 2; ++i) { for (int j = 0; j < 2; ++j) - printf("%8f ", resErr_loc.At(0,i,j)); printf("\n"); - } printf("\n"); - printf("K:\n"); - for (int i = 0; i < 6; ++i) { for (int j = 0; j < 3; ++j) - printf("%8f ", K.At(0,i,j)); printf("\n"); - } printf("\n"); - printf("outPar:\n"); - for (int i = 0; i < 6; ++i) { - printf("%8f ", outPar.At(0,i,0)); - } printf("\n"); - printf("outErr:\n"); - for (int i = 0; i < 6; ++i) { for (int j = 0; j < 6; ++j) - printf("%8f ", outErr.At(0,i,j)); printf("\n"); - } printf("\n"); - } -#endif +void kalmanComputeChi2(const MPlexLS &psErr, const MPlexLV& psPar, const MPlexQI &inChg, + const MPlexHS &msErr, const MPlexHV& msPar, + MPlexQF& outChi2, + const int N_proc) +{ + kalmanOperation(KFO_Calculate_Chi2, psErr, psPar, msErr, msPar, + dummy_err, dummy_par, outChi2, N_proc); } -void computeChi2MPlex(const MPlexLS &psErr, const MPlexLV& psPar, const MPlexQI &inChg, - const MPlexHS &msErr, const MPlexHV& msPar, - MPlexQF& outChi2, - const int N_proc) +void kalmanPropagateAndComputeChi2(const MPlexLS &psErr, const MPlexLV& psPar, const MPlexQI &inChg, + const MPlexHS &msErr, const MPlexHV& msPar, + MPlexQF& outChi2, + const int N_proc, const PropagationFlags propFlags) { - - // const idx_t N = psErr.N; - // Assert N-s of all parameters are the same. - - // Temporaries -- this is expensive -- should have them allocated outside and reused. - // Can be passed in in a struct, see above. - - // updateParametersContext ctx; - //assert((long long)(&updateCtx.propErr.fArray[0]) % 64 == 0); - MPlexLS propErr; MPlexLV propPar; - // do a full propagation step to correct for residual distance from the hit radius - need the charge for this + if (Config::useCMSGeom) { - propagateHelixToRMPlex(psErr, psPar, inChg, msPar, propErr, propPar, N_proc); + MPlexQF msRad; +#pragma simd + for (int n = 0; n < NN; ++n) + { + msRad.At(n, 0, 0) = std::hypot(msPar.ConstAt(n, 0, 0), msPar.ConstAt(n, 1, 0)); + } + + propagateHelixToRMPlex(psErr, psPar, inChg, msRad, propErr, propPar, N_proc, propFlags); } else { propErr = psErr; propPar = psPar; } + kalmanOperation(KFO_Calculate_Chi2, propErr, propPar, msErr, msPar, + dummy_err, dummy_par, outChi2, N_proc); +} + +//------------------------------------------------------------------------------ + +void kalmanOperation(const int kfOp, + const MPlexLS &psErr, const MPlexLV& psPar, + const MPlexHS &msErr, const MPlexHV& msPar, + MPlexLS &outErr, MPlexLV& outPar, MPlexQF& outChi2, + const int N_proc) +{ #ifdef DEBUG { dmutex_guard; - printf("propPar:\n"); + printf("psPar:\n"); for (int i = 0; i < 6; ++i) { - printf("%8f ", propPar.ConstAt(0,0,i)); printf("\n"); + printf("%8f ", psPar.ConstAt(0,0,i)); printf("\n"); } printf("\n"); - printf("propErr:\n"); + printf("psErr:\n"); for (int i = 0; i < 6; ++i) { for (int j = 0; j < 6; ++j) - printf("%8f ", propErr.At(0,i,j)); printf("\n"); + printf("%8f ", psErr.At(0,i,j)); printf("\n"); } printf("\n"); printf("msPar:\n"); for (int i = 0; i < 3; ++i) { @@ -635,16 +574,16 @@ void computeChi2MPlex(const MPlexLS &psErr, const MPlexLV& psPar, const MPlexQI MPlexQF rotT00; MPlexQF rotT01; for (int n = 0; n < NN; ++n) { - const float r = hipo(msPar.ConstAt(n, 0, 0), msPar.ConstAt(n, 1, 0)); - rotT00.At(n, 0, 0) = -(msPar.ConstAt(n, 1, 0)+propPar.ConstAt(n, 1, 0))/(2*r); - rotT01.At(n, 0, 0) = (msPar.ConstAt(n, 0, 0)+propPar.ConstAt(n, 0, 0))/(2*r); + const float r = std::hypot(msPar.ConstAt(n, 0, 0), msPar.ConstAt(n, 1, 0)); + rotT00.At(n, 0, 0) = -(msPar.ConstAt(n, 1, 0) + psPar.ConstAt(n, 1, 0)) / (2*r); + rotT01.At(n, 0, 0) = (msPar.ConstAt(n, 0, 0) + psPar.ConstAt(n, 0, 0)) / (2*r); } MPlexHV res_glo; //position residual in global coordinates - SubtractFirst3(msPar, propPar, res_glo); - + SubtractFirst3(msPar, psPar, res_glo); + MPlexHS resErr_glo;//covariance sum in global position coordinates - AddIntoUpperLeft3x3(propErr, msErr, resErr_glo); + AddIntoUpperLeft3x3(psErr, msErr, resErr_glo); MPlex2V res_loc; //position residual in local coordinates RotateResidulsOnTangentPlane(rotT00,rotT01,res_glo,res_loc); @@ -666,179 +605,174 @@ void computeChi2MPlex(const MPlexLS &psErr, const MPlexLV& psPar, const MPlexQI //invert the 2x2 matrix Matriplex::InvertCramerSym(resErr_loc); - //compute chi2 - Chi2Similarity(res_loc, resErr_loc, outChi2); + if (kfOp & KFO_Calculate_Chi2) + { + Chi2Similarity(res_loc, resErr_loc, outChi2); #ifdef DEBUG - { - dmutex_guard; - printf("resErr_loc (Inv):\n"); - for (int i = 0; i < 2; ++i) { - for (int j = 0; j < 2; ++j) - printf("%8f ", resErr_loc.At(0,i,j)); printf("\n"); - } printf("\n"); - printf("chi2: %8f\n", outChi2.At(0,0,0)); - } + { + dmutex_guard; + printf("resErr_loc (Inv):\n"); + for (int i = 0; i < 2; ++i) { + for (int j = 0; j < 2; ++j) + printf("%8f ", resErr_loc.At(0,i,j)); printf("\n"); + } printf("\n"); + printf("chi2: %8f\n", outChi2.At(0,0,0)); + } #endif + } -} + if (kfOp & KFO_Update_Params) + { + MPlexLH K; // kalman gain, fixme should be L2 + KalmanHTG(rotT00, rotT01, resErr_loc, tempHH); // intermediate term to get kalman gain (H^T*G) + KalmanGain(psErr, tempHH, K); + MultResidualsAdd(K, psPar, res_loc, outPar); + MPlexLL tempLL; + squashPhiMPlex(outPar,N_proc); // ensure phi is between |pi| + KHMult(K, rotT00, rotT01, tempLL); + KHC(tempLL, psErr, outErr); + outErr.Subtract(psErr, outErr); +#ifdef DEBUG + { + dmutex_guard; + printf("res_glo:\n"); + for (int i = 0; i < 3; ++i) { + printf("%8f ", res_glo.At(0,i,0)); + } printf("\n"); + printf("res_loc:\n"); + for (int i = 0; i < 2; ++i) { + printf("%8f ", res_loc.At(0,i,0)); + } printf("\n"); + printf("resErr_loc (Inv):\n"); + for (int i = 0; i < 2; ++i) { for (int j = 0; j < 2; ++j) + printf("%8f ", resErr_loc.At(0,i,j)); printf("\n"); + } printf("\n"); + printf("K:\n"); + for (int i = 0; i < 6; ++i) { for (int j = 0; j < 3; ++j) + printf("%8f ", K.At(0,i,j)); printf("\n"); + } printf("\n"); + printf("outPar:\n"); + for (int i = 0; i < 6; ++i) { + printf("%8f ", outPar.At(0,i,0)); + } printf("\n"); + printf("outErr:\n"); + for (int i = 0; i < 6; ++i) { for (int j = 0; j < 6; ++j) + printf("%8f ", outErr.At(0,i,j)); printf("\n"); + } printf("\n"); + } +#endif + } +} -void updateParametersEndcapMPlex(const MPlexLS &psErr, const MPlexLV& psPar, const MPlexQI &inChg, - const MPlexHS &msErr, const MPlexHV& msPar, - MPlexLS &outErr, MPlexLV& outPar, - const int N_proc) -{ - // const idx_t N = psErr.N; - // Assert N-s of all parameters are the same. - // Temporaries -- this is expensive -- should have them allocated outside and reused. - // Can be passed in in a struct, see above. +//============================================================================== +// Kalman operations - Endcap +//============================================================================== - // updateParametersContext ctx; - //assert((long long)(&updateCtx.propErr.fArray[0]) % 64 == 0); +void kalmanUpdateEndcap(const MPlexLS &psErr, const MPlexLV& psPar, + const MPlexHS &msErr, const MPlexHV& msPar, + MPlexLS &outErr, MPlexLV& outPar, + const int N_proc) +{ + kalmanOperationEndcap(KFO_Update_Params, psErr, psPar, msErr, msPar, + outErr, outPar, dummy_chi2, N_proc); +} +void kalmanPropagateAndUpdateEndcap(const MPlexLS &psErr, const MPlexLV& psPar, const MPlexQI &inChg, + const MPlexHS &msErr, const MPlexHV& msPar, + MPlexLS &outErr, MPlexLV& outPar, + const int N_proc, const PropagationFlags propFlags) +{ MPlexLS propErr; MPlexLV propPar; - // do a full propagation step to correct for residual distance from the hit radius - need the charge for this + + // XXXX The following if could go away if we use proper update_foos in + // steering and set them from geom. For now it is still here (I'm going + // crazy already as things are). + if (Config::useCMSGeom) { - propagateHelixToZMPlex(psErr, psPar, inChg, msPar, propErr, propPar, N_proc); + MPlexQF msZ; +#pragma simd + for (int n = 0; n < NN; ++n) + { + msZ.At(n, 0, 0) = msPar.ConstAt(n, 2, 0); + } + + propagateHelixToZMPlex(psErr, psPar, inChg, msZ, propErr, propPar, N_proc, propFlags); } else { propErr = psErr; propPar = psPar; } -#ifdef DEBUG - { - dmutex_guard; - printf("updateParametersEndcapMPlex\n"); - printf("propPar:\n"); - for (int i = 0; i < 6; ++i) { - printf("%8f ", propPar.ConstAt(0,0,i)); printf("\n"); - } printf("\n"); - printf("msPar:\n"); - for (int i = 0; i < 3; ++i) { - printf("%8f ", msPar.ConstAt(0,0,i)); printf("\n"); - } printf("\n"); - printf("propErr:\n"); - for (int i = 0; i < 6; ++i) { for (int j = 0; j < 6; ++j) - printf("%8f ", propErr.At(0,i,j)); printf("\n"); - } printf("\n"); - printf("msErr:\n"); - for (int i = 0; i < 3; ++i) { for (int j = 0; j < 3; ++j) - printf("%8f ", msErr.ConstAt(0,i,j)); printf("\n"); - } printf("\n"); - } -#endif - - MPlex2V res; - SubtractFirst2(msPar, propPar, res); - - MPlex2S resErr; - AddIntoUpperLeft2x2(propErr, msErr, resErr); - -#ifdef DEBUG - { - dmutex_guard; - printf("resErr:\n"); - for (int i = 0; i < 2; ++i) { for (int j = 0; j < 2; ++j) - printf("%8f ", resErr.At(0,i,j)); printf("\n"); - } printf("\n"); - } -#endif - - //invert the 2x2 matrix - Matriplex::InvertCramerSym(resErr); - - MPlexL2 K; - KalmanGain(propErr, resErr, K); - - MultResidualsAdd(K, propPar, res, outPar); - - squashPhiMPlex(outPar,N_proc); // ensure phi is between |pi| - - KHC(K, propErr, outErr); - -#ifdef DEBUG - { - printf("outErr before subtract:\n"); - for (int i = 0; i < 6; ++i) { for (int j = 0; j < 6; ++j) - printf("%8f ", outErr.At(0,i,j)); printf("\n"); - } printf("\n"); - } -#endif + kalmanOperationEndcap(KFO_Update_Params, propErr, propPar, msErr, msPar, + outErr, outPar, dummy_chi2, N_proc); +} - outErr.Subtract(propErr, outErr); +//------------------------------------------------------------------------------ -#ifdef DEBUG - { - dmutex_guard; - printf("res:\n"); - for (int i = 0; i < 2; ++i) { - printf("%8f ", res.At(0,i,0)); - } printf("\n"); - printf("resErr (Inv):\n"); - for (int i = 0; i < 2; ++i) { for (int j = 0; j < 2; ++j) - printf("%8f ", resErr.At(0,i,j)); printf("\n"); - } printf("\n"); - printf("K:\n"); - for (int i = 0; i < 6; ++i) { for (int j = 0; j < 2; ++j) - printf("%8f ", K.At(0,i,j)); printf("\n"); - } printf("\n"); - printf("outPar:\n"); - for (int i = 0; i < 6; ++i) { - printf("%8f ", outPar.At(0,i,0)); - } printf("\n"); - printf("outErr:\n"); - for (int i = 0; i < 6; ++i) { for (int j = 0; j < 6; ++j) - printf("%8f ", outErr.At(0,i,j)); printf("\n"); - } printf("\n"); - } -#endif +void kalmanComputeChi2Endcap(const MPlexLS &psErr, const MPlexLV& psPar, const MPlexQI &inChg, + const MPlexHS &msErr, const MPlexHV& msPar, + MPlexQF& outChi2, + const int N_proc) +{ + kalmanOperationEndcap(KFO_Calculate_Chi2, psErr, psPar, msErr, msPar, + dummy_err, dummy_par, outChi2, N_proc); } -void computeChi2EndcapMPlex(const MPlexLS &psErr, const MPlexLV& psPar, const MPlexQI &inChg, - const MPlexHS &msErr, const MPlexHV& msPar, - MPlexQF& outChi2, - const int N_proc) +void kalmanPropagateAndComputeChi2Endcap(const MPlexLS &psErr, const MPlexLV& psPar, const MPlexQI &inChg, + const MPlexHS &msErr, const MPlexHV& msPar, + MPlexQF& outChi2, + const int N_proc, const PropagationFlags propFlags) { - // const idx_t N = psErr.N; - // Assert N-s of all parameters are the same. - - // Temporaries -- this is expensive -- should have them allocated outside and reused. - // Can be passed in in a struct, see above. - - // updateParametersContext ctx; - //assert((long long)(&updateCtx.propErr.fArray[0]) % 64 == 0); - MPlexLS propErr; MPlexLV propPar; - // do a full propagation step to correct for residual distance from the hit radius - need the charge for this + if (Config::useCMSGeom) { - propagateHelixToZMPlex(psErr, psPar, inChg, msPar, propErr, propPar, N_proc); + MPlexQF msZ; +#pragma simd + for (int n = 0; n < NN; ++n) + { + msZ.At(n, 0, 0) = msPar.ConstAt(n, 2, 0); + } + + propagateHelixToZMPlex(psErr, psPar, inChg, msZ, propErr, propPar, N_proc, propFlags); } else { propErr = psErr; propPar = psPar; } + kalmanOperationEndcap(KFO_Calculate_Chi2, propErr, propPar, msErr, msPar, + dummy_err, dummy_par, outChi2, N_proc); +} + +//------------------------------------------------------------------------------ + +void kalmanOperationEndcap(const int kfOp, + const MPlexLS &psErr, const MPlexLV& psPar, + const MPlexHS &msErr, const MPlexHV& msPar, + MPlexLS &outErr, MPlexLV& outPar, MPlexQF& outChi2, + const int N_proc) +{ #ifdef DEBUG { dmutex_guard; - printf("computeChi2EndcapMPlex\n"); - printf("propPar:\n"); + printf("updateParametersEndcapMPlex\n"); + printf("psPar:\n"); for (int i = 0; i < 6; ++i) { - printf("%8f ", propPar.ConstAt(0,0,i)); printf("\n"); + printf("%8f ", psPar.ConstAt(0,0,i)); printf("\n"); } printf("\n"); printf("msPar:\n"); for (int i = 0; i < 3; ++i) { printf("%8f ", msPar.ConstAt(0,0,i)); printf("\n"); } printf("\n"); - printf("propErr:\n"); + printf("psErr:\n"); for (int i = 0; i < 6; ++i) { for (int j = 0; j < 6; ++j) - printf("%8f ", propErr.At(0,i,j)); printf("\n"); + printf("%8f ", psErr.At(0,i,j)); printf("\n"); } printf("\n"); printf("msErr:\n"); for (int i = 0; i < 3; ++i) { for (int j = 0; j < 3; ++j) @@ -848,18 +782,14 @@ void computeChi2EndcapMPlex(const MPlexLS &psErr, const MPlexLV& psPar, const M #endif MPlex2V res; - SubtractFirst2(msPar, propPar, res); + SubtractFirst2(msPar, psPar, res); MPlex2S resErr; - AddIntoUpperLeft2x2(propErr, msErr, resErr); + AddIntoUpperLeft2x2(psErr, msErr, resErr); #ifdef DEBUG { dmutex_guard; - printf("res:\n"); - for (int i = 0; i < 2; ++i) { - printf("%8f ", res.At(0,0,i)); printf("\n"); - } printf("\n"); printf("resErr:\n"); for (int i = 0; i < 2; ++i) { for (int j = 0; j < 2; ++j) printf("%8f ", resErr.At(0,i,j)); printf("\n"); @@ -870,19 +800,69 @@ void computeChi2EndcapMPlex(const MPlexLS &psErr, const MPlexLV& psPar, const M //invert the 2x2 matrix Matriplex::InvertCramerSym(resErr); - //compute chi2 - Chi2Similarity(res, resErr, outChi2); + if (kfOp & KFO_Calculate_Chi2) + { + Chi2Similarity(res, resErr, outChi2); #ifdef DEBUG - { - dmutex_guard; - printf("resErr_loc (Inv):\n"); - for (int i = 0; i < 2; ++i) { - for (int j = 0; j < 2; ++j) - printf("%8f ", resErr.At(0,i,j)); printf("\n"); - } printf("\n"); - printf("chi2: %8f\n", outChi2.At(0,0,0)); + { + dmutex_guard; + printf("resErr_loc (Inv):\n"); + for (int i = 0; i < 2; ++i) { + for (int j = 0; j < 2; ++j) + printf("%8f ", resErr.At(0,i,j)); printf("\n"); + } printf("\n"); + printf("chi2: %8f\n", outChi2.At(0,0,0)); + } +#endif } + + if (kfOp & KFO_Update_Params) + { + MPlexL2 K; + KalmanGain(psErr, resErr, K); + + MultResidualsAdd(K, psPar, res, outPar); + + squashPhiMPlex(outPar,N_proc); // ensure phi is between |pi| + + KHC(K, psErr, outErr); + +#ifdef DEBUG + { + printf("outErr before subtract:\n"); + for (int i = 0; i < 6; ++i) { for (int j = 0; j < 6; ++j) + printf("%8f ", outErr.At(0,i,j)); printf("\n"); + } printf("\n"); + } #endif + outErr.Subtract(psErr, outErr); + +#ifdef DEBUG + { + dmutex_guard; + printf("res:\n"); + for (int i = 0; i < 2; ++i) { + printf("%8f ", res.At(0,i,0)); + } printf("\n"); + printf("resErr (Inv):\n"); + for (int i = 0; i < 2; ++i) { for (int j = 0; j < 2; ++j) + printf("%8f ", resErr.At(0,i,j)); printf("\n"); + } printf("\n"); + printf("K:\n"); + for (int i = 0; i < 6; ++i) { for (int j = 0; j < 2; ++j) + printf("%8f ", K.At(0,i,j)); printf("\n"); + } printf("\n"); + printf("outPar:\n"); + for (int i = 0; i < 6; ++i) { + printf("%8f ", outPar.At(0,i,0)); + } printf("\n"); + printf("outErr:\n"); + for (int i = 0; i < 6; ++i) { for (int j = 0; j < 6; ++j) + printf("%8f ", outErr.At(0,i,j)); printf("\n"); + } printf("\n"); + } +#endif + } } diff --git a/mkFit/KalmanUtilsMPlex.h b/mkFit/KalmanUtilsMPlex.h index 53234074e2178..ca4d32c791eb4 100644 --- a/mkFit/KalmanUtilsMPlex.h +++ b/mkFit/KalmanUtilsMPlex.h @@ -8,30 +8,75 @@ #include "FitterCU.h" #endif -void updateParametersMPlex(const MPlexLS &psErr, const MPlexLV& psPar, const MPlexQI &inChg, +//------------------------------------------------------------------------------ + +enum KalmanFilterOperation +{ + KFO_Calculate_Chi2 = 1, + KFO_Update_Params = 2 +}; + + +//------------------------------------------------------------------------------ + +void kalmanUpdate(const MPlexLS &psErr, const MPlexLV& psPar, + const MPlexHS &msErr, const MPlexHV& msPar, + MPlexLS &outErr, MPlexLV& outPar, + const int N_proc); + +void kalmanPropagateAndUpdate(const MPlexLS &psErr, const MPlexLV& psPar, const MPlexQI &inChg, + const MPlexHS &msErr, const MPlexHV& msPar, + MPlexLS &outErr, MPlexLV& outPar, + const int N_proc, const PropagationFlags propFlags); + + +void kalmanComputeChi2(const MPlexLS &psErr, const MPlexLV& psPar, const MPlexQI &inChg, + const MPlexHS &msErr, const MPlexHV& msPar, + MPlexQF& outChi2, + const int N_proc); + +void kalmanPropagateAndComputeChi2(const MPlexLS &psErr, const MPlexLV& psPar, const MPlexQI &inChg, + const MPlexHS &msErr, const MPlexHV& msPar, + MPlexQF& outChi2, + const int N_proc, const PropagationFlags propFlags); + + +void kalmanOperation(const int kfOp, + const MPlexLS &psErr, const MPlexLV& psPar, + const MPlexHS &msErr, const MPlexHV& msPar, + MPlexLS &outErr, MPlexLV& outPar, MPlexQF& outChi2, + const int N_proc); + +//------------------------------------------------------------------------------ + + +void kalmanUpdateEndcap(const MPlexLS &psErr, const MPlexLV& psPar, + const MPlexHS &msErr, const MPlexHV& msPar, + MPlexLS &outErr, MPlexLV& outPar, + const int N_proc); + +void kalmanPropagateAndUpdateEndcap(const MPlexLS &psErr, const MPlexLV& psPar, const MPlexQI &inChg, + const MPlexHS &msErr, const MPlexHV& msPar, + MPlexLS &outErr, MPlexLV& outPar, + const int N_proc, const PropagationFlags propFlags); + + +void kalmanComputeChi2Endcap(const MPlexLS &psErr, const MPlexLV& psPar, const MPlexQI &inChg, + const MPlexHS &msErr, const MPlexHV& msPar, + MPlexQF& outChi2, + const int N_proc); + +void kalmanPropagateAndComputeChi2Endcap(const MPlexLS &psErr, const MPlexLV& psPar, const MPlexQI &inChg, + const MPlexHS &msErr, const MPlexHV& msPar, + MPlexQF& outChi2, + const int N_proc, const PropagationFlags propFlags); + + +void kalmanOperationEndcap(const int kfOp, + const MPlexLS &psErr, const MPlexLV& psPar, const MPlexHS &msErr, const MPlexHV& msPar, - MPlexLS &outErr, MPlexLV& outPar, + MPlexLS &outErr, MPlexLV& outPar, MPlexQF& outChi2, const int N_proc); -#ifdef USE_CUDA // FIXME: temporary; move to FitterCU -void computeChi2MPlex_tmp(const MPlexLS &psErr, const MPlexLV& psPar, const MPlexQI &inChg, - const MPlexHS &msErr, const MPlexHV& msPar, - MPlexQF& outChi2, - FitterCU<float>& cuFitter); -#endif -void computeChi2MPlex(const MPlexLS &psErr, const MPlexLV& psPar, const MPlexQI &inChg, - const MPlexHS &msErr, const MPlexHV& msPar, - MPlexQF& outChi2, - const int N_proc); - -void updateParametersEndcapMPlex(const MPlexLS &psErr, const MPlexLV& psPar, const MPlexQI &inChg, - const MPlexHS &msErr, const MPlexHV& msPar, - MPlexLS &outErr, MPlexLV& outPar, - const int N_proc); - -void computeChi2EndcapMPlex(const MPlexLS &psErr, const MPlexLV& psPar, const MPlexQI &inChg, - const MPlexHS &msErr, const MPlexHV& msPar, - MPlexQF& outChi2, - const int N_proc); #endif diff --git a/mkFit/MkBase.h b/mkFit/MkBase.h index f4375aa07f346..acfedee8c38c2 100644 --- a/mkFit/MkBase.h +++ b/mkFit/MkBase.h @@ -5,6 +5,10 @@ #include "PropagationMPlex.h" +//============================================================================== +// MkBase +//============================================================================== + struct MkBase { MPlexLS Err[2]; @@ -21,16 +25,60 @@ struct MkBase MkBase() {} - void PropagateTracksToR(float R, const int N_proc) + //---------------------------------------------------------------------------- + + void PropagateTracksToR(float r, const int N_proc, const PropagationFlags pf) { - propagateHelixToRMPlex(Err[iC], Par[iC], Chg, R, - Err[iP], Par[iP], N_proc); + MPlexQF msRad; +#pragma simd + for (int n = 0; n < NN; ++n) + { + msRad.At(n, 0, 0) = r; + } + + propagateHelixToRMPlex(Err[iC], Par[iC], Chg, msRad, + Err[iP], Par[iP], N_proc, pf); } - void PropagateTracksToZ(float Z, const int N_proc) + void PropagateTracksToHitR(const MPlexHV& par, const int N_proc, const PropagationFlags pf) { - propagateHelixToZMPlex(Err[iC], Par[iC], Chg, Z, - Err[iP], Par[iP], N_proc); + MPlexQF msRad; +#pragma simd + for (int n = 0; n < NN; ++n) + { + msRad.At(n, 0, 0) = std::hypot(par.ConstAt(n, 0, 0), par.ConstAt(n, 1, 0)); + } + + propagateHelixToRMPlex(Err[iC], Par[iC], Chg, msRad, + Err[iP], Par[iP], N_proc, pf); + } + + //---------------------------------------------------------------------------- + + void PropagateTracksToZ(float z, const int N_proc, const PropagationFlags pf) + { + MPlexQF msZ; +#pragma simd + for (int n = 0; n < NN; ++n) + { + msZ.At(n, 0, 0) = z; + } + + propagateHelixToZMPlex(Err[iC], Par[iC], Chg, msZ, + Err[iP], Par[iP], N_proc, pf); + } + + void PropagateTracksToHitZ(const MPlexHV& par, const int N_proc, const PropagationFlags pf) + { + MPlexQF msZ; +#pragma simd + for (int n = 0; n < NN; ++n) + { + msZ.At(n, 0, 0) = par.ConstAt(n, 2, 0); + } + + propagateHelixToZMPlex(Err[iC], Par[iC], Chg, msZ, + Err[iP], Par[iP], N_proc, pf); } }; diff --git a/mkFit/MkBuilder.cc b/mkFit/MkBuilder.cc index 1e971c5c30f5d..23aa1076dd930 100644 --- a/mkFit/MkBuilder.cc +++ b/mkFit/MkBuilder.cc @@ -176,8 +176,8 @@ MkBuilder::MkBuilder() : { const TrackerInfo &ti = Config::TrkInfo; - m_fndfoos_brl = { computeChi2MPlex, updateParametersMPlex, &MkBase::PropagateTracksToR }; - m_fndfoos_ec = { computeChi2EndcapMPlex, updateParametersEndcapMPlex, &MkBase::PropagateTracksToZ }; + m_fndfoos_brl = { kalmanPropagateAndComputeChi2, kalmanPropagateAndUpdate, &MkBase::PropagateTracksToR }; + m_fndfoos_ec = { kalmanPropagateAndComputeChi2Endcap, kalmanPropagateAndUpdateEndcap, &MkBase::PropagateTracksToZ }; { SteeringParams &sp = m_steering_params[TrackerInfo::Reg_Endcap_Neg]; sp.reserve_plan(3 + 3 + 6 + 18); @@ -1250,7 +1250,9 @@ void MkBuilder::FindTracksBestHit() dcall(pre_prop_print(curr_layer, mkfndr.get())); - (mkfndr.get()->*fnd_foos.m_propagate_foo)(layer_info.m_propagate_to, curr_tridx); + (mkfndr.get()->*fnd_foos.m_propagate_foo)(layer_info.m_propagate_to, curr_tridx, + PF_use_param_b_field | PF_apply_material); + // XXXX-MFMAT-review dcall(post_prop_print(curr_layer, mkfndr.get())); @@ -1415,7 +1417,10 @@ void MkBuilder::FindTracksStandard() //propagate to layer dcall(pre_prop_print(curr_layer, mkfndr.get())); - (mkfndr.get()->*fnd_foos.m_propagate_foo)(layer_info.m_propagate_to, end - itrack); + (mkfndr.get()->*fnd_foos.m_propagate_foo)(layer_info.m_propagate_to, end - itrack, + PF_use_param_b_field | PF_apply_material); + // XXXX-MFMAT-review + dcall(post_prop_print(curr_layer, mkfndr.get())); @@ -1613,7 +1618,10 @@ void MkBuilder::find_tracks_in_layers(CandCloner &cloner, MkFinder *mkfndr, #endif // propagate to current layer - (mkfndr->*fnd_foos.m_propagate_foo)(layer_info.m_propagate_to, end - itrack); + (mkfndr->*fnd_foos.m_propagate_foo)(layer_info.m_propagate_to, end - itrack, + PF_use_param_b_field | PF_apply_material); + // XXXX-MFMAT-review + // copy_out the propagated track params, errors only (hit-idcs and chi2 already updated) mkfndr->CopyOutParErr(eoccs.m_candidates, end - itrack, true); diff --git a/mkFit/MkFinder.cc b/mkFit/MkFinder.cc index 219826ec82607..8a433b93b35d0 100644 --- a/mkFit/MkFinder.cc +++ b/mkFit/MkFinder.cc @@ -434,7 +434,7 @@ void MkFinder::AddBestHit(const LayerOfHits &layer_of_hits, const int N_proc, //now compute the chi2 of track state vs hit MPlexQF outChi2; (*fnd_foos.m_compute_chi2_foo)(Err[iP], Par[iP], Chg, msErr, msPar, - outChi2, N_proc); + outChi2, N_proc, PF_use_param_b_field); #ifndef NO_PREFETCH // Prefetch to L1 the hits we'll process in the next loop iteration. @@ -536,7 +536,7 @@ void MkFinder::AddBestHit(const LayerOfHits &layer_of_hits, const int N_proc, dprint("update parameters"); (*fnd_foos.m_update_param_foo)(Err[iP], Par[iP], Chg, msErr, msPar, - Err[iC], Par[iC], N_proc); + Err[iC], Par[iC], N_proc, PF_use_param_b_field); //std::cout << "Par[iP](0,0,0)=" << Par[iP](0,0,0) << " Par[iC](0,0,0)=" << Par[iC](0,0,0)<< std::endl; } @@ -616,7 +616,7 @@ void MkFinder::FindCandidates(const LayerOfHits &layer_of_hits, //now compute the chi2 of track state vs hit MPlexQF outChi2; (*fnd_foos.m_compute_chi2_foo)(Err[iP], Par[iP], Chg, msErr, msPar, - outChi2, N_proc); + outChi2, N_proc, PF_use_param_b_field); // Prefetch to L1 the hits we'll (probably) process in the next loop iteration. for (int itrack = 0; itrack < N_proc; ++itrack) @@ -648,7 +648,7 @@ void MkFinder::FindCandidates(const LayerOfHits &layer_of_hits, if (oneCandPassCut) { (*fnd_foos.m_update_param_foo)(Err[iP], Par[iP], Chg, msErr, msPar, - Err[iC], Par[iC], N_proc); + Err[iC], Par[iC], N_proc, PF_use_param_b_field); dprint("update parameters" << std::endl << "propagated track parameters x=" << Par[iP].ConstAt(0, 0, 0) << " y=" << Par[iP].ConstAt(0, 1, 0) << std::endl @@ -788,7 +788,7 @@ void MkFinder::FindCandidatesCloneEngine(const LayerOfHits &layer_of_hits, CandC //now compute the chi2 of track state vs hit MPlexQF outChi2; - (*fnd_foos.m_compute_chi2_foo)(Err[iP], Par[iP], Chg, msErr, msPar, outChi2, N_proc); + (*fnd_foos.m_compute_chi2_foo)(Err[iP], Par[iP], Chg, msErr, msPar, outChi2, N_proc, PF_use_param_b_field); // Prefetch to L1 the hits we'll (probably) process in the next loop iteration. for (int itrack = 0; itrack < N_proc; ++itrack) @@ -882,7 +882,7 @@ void MkFinder::UpdateWithLastHit(const LayerOfHits &layer_of_hits, int N_proc, } (*fnd_foos.m_update_param_foo)(Err[iP], Par[iP], Chg, msErr, msPar, - Err[iC], Par[iC], N_proc); + Err[iC], Par[iC], N_proc, PF_use_param_b_field); //now that we have moved propagation at the end of the sequence we lost the handle of //using the propagated parameters instead of the updated for the missing hit case. @@ -995,7 +995,9 @@ void MkFinder::BkFitFitTracks(const EventOfHits & eventofhits, { // Prototyping final backward fit. // This works with track-finding indices, before remapping. - // layers should be collected during track finding and list all layers that have actual hits. + // + // Layers should be collected during track finding and list all layers that have actual hits. + // Then we could avoid checking which layers actually do have hits. MPlexQF tmp_chi2; float tmp_err[6] = { 666, 0, 666, 0, 0, 666 }; @@ -1037,27 +1039,17 @@ void MkFinder::BkFitFitTracks(const EventOfHits & eventofhits, if (LI.is_barrel()) { - // ZZZ for cmsGeom propagation done in also done in update! fix - propagateHelixToRMPlex(Err[iC], Par[iC], Chg, msPar, - Err[iP], Par[iP], N_proc, useParamBfield); + PropagateTracksToHitR(msPar, N_proc, PF_use_param_b_field | PF_apply_material); - // ZZZ should have common chi2 / update function - computeChi2MPlex(Err[iP], Par[iP], Chg, msErr, msPar, - tmp_chi2, N_proc); - - updateParametersMPlex(Err[iP], Par[iP], Chg, msErr, msPar, - Err[iC], Par[iC], N_proc); + kalmanOperation(KFO_Calculate_Chi2 | KFO_Update_Params, + Err[iP], Par[iP], msErr, msPar, Err[iC], Par[iC], tmp_chi2, N_proc); } else { - propagateHelixToZMPlex(Err[iC], Par[iC], Chg, msPar, - Err[iP], Par[iP], N_proc, useParamBfield); - - computeChi2EndcapMPlex(Err[iP], Par[iP], Chg, msErr, msPar, - tmp_chi2, N_proc); + PropagateTracksToHitZ(msPar, N_proc, PF_use_param_b_field | PF_apply_material); - updateParametersEndcapMPlex(Err[iP], Par[iP], Chg, msErr, msPar, - Err[iC], Par[iC], N_proc); + kalmanOperationEndcap(KFO_Calculate_Chi2 | KFO_Update_Params, + Err[iP], Par[iP], msErr, msPar, Err[iC], Par[iC], tmp_chi2, N_proc); } #ifdef DEBUG_BACKWARD_FIT diff --git a/mkFit/MkFinder.h b/mkFit/MkFinder.h index 023d4a451c0ee..0b83189b9c9c2 100644 --- a/mkFit/MkFinder.h +++ b/mkFit/MkFinder.h @@ -77,10 +77,17 @@ class MkFinder : public MkBase MPlexQI XHitSize; MPlexHitIdx XHitArr; - // Hit error / parameters for hit matching, update. + // Hit errors / parameters for hit matching, update. MPlexHS msErr; MPlexHV msPar; + // An idea: Do propagation to hit in FindTracksXYZZ functions. + // Have some state / functions here that make this short to write. + // This would simplify KalmanUtils (remove the propagate functions). + // Track errors / parameters propagated to current hit. + // MPlexLS candErrAtCurrHit; + // MPlexLV candParAtCurrHit; + //============================================================================ MkFinder() {} diff --git a/mkFit/MkFitter.cc b/mkFit/MkFitter.cc index d6f58e3f9703b..ea036d65fd8a6 100644 --- a/mkFit/MkFitter.cc +++ b/mkFit/MkFitter.cc @@ -444,11 +444,10 @@ void MkFitter::FitTracksWithInterSlurp(const std::vector<HitVec>& layersohits, msErr[0].SlurpIn(varr + off_error, idx, N_proc); #endif - propagateHelixToRMPlex(Err[iC], Par[iC], Chg, msPar[0], - Err[iP], Par[iP], N_proc); + PropagateTracksToHitR(msPar[0], N_proc, PF_use_param_b_field | PF_apply_material); - updateParametersMPlex(Err[iP], Par[iP], Chg, msErr[0], msPar[0], - Err[iC], Par[iC], N_proc); + kalmanUpdate(Err[iP], Par[iP], msErr[0], msPar[0], + Err[iC], Par[iC], N_proc); } } @@ -521,17 +520,18 @@ void MkFitter::FitTracks(const int N_proc, const Event * ev, const bool useParam { // Fitting loop. + PropagationFlags pflags((useParamBfield ? PF_use_param_b_field : 0) | PF_apply_material); + for (int hi = 0; hi < Nhits; ++hi) { // Note, charge is not passed (line propagation). // propagateLineToRMPlex(Err[iC], Par[iC], msErr[hi], msPar[hi], // Err[iP], Par[iP]); - propagateHelixToRMPlex(Err[iC], Par[iC], Chg, msPar[hi], - Err[iP], Par[iP], N_proc, useParamBfield); + PropagateTracksToHitR(msPar[hi], N_proc, pflags); - updateParametersMPlex(Err[iP], Par[iP], Chg, msErr[hi], msPar[hi], - Err[iC], Par[iC], N_proc); + kalmanUpdate(Err[iP], Par[iP], msErr[hi], msPar[hi], + Err[iC], Par[iC], N_proc); if (Config::fit_val) MkFitter::CollectFitValidation(hi,N_proc,ev); } @@ -567,6 +567,8 @@ void MkFitter::FitTracksSteered(const bool is_barrel[], const int N_proc, const dprintf("MkFitter::FitTracksSteered %d %d %d\n", is_barrel[0], is_barrel[1], is_barrel[2]); + PropagationFlags pflags((useParamBfield ? PF_use_param_b_field : 0) | PF_apply_material); + for (int hi = 0; hi < Nhits; ++hi) { // Note, charge is not passed (line propagation). @@ -575,19 +577,17 @@ void MkFitter::FitTracksSteered(const bool is_barrel[], const int N_proc, const if (is_barrel[hi]) { - propagateHelixToRMPlex(Err[iC], Par[iC], Chg, msPar[hi], - Err[iP], Par[iP], N_proc, useParamBfield); + PropagateTracksToHitR(msPar[hi], N_proc, pflags); - updateParametersMPlex(Err[iP], Par[iP], Chg, msErr[hi], msPar[hi], - Err[iC], Par[iC], N_proc); + kalmanUpdate(Err[iP], Par[iP], msErr[hi], msPar[hi], + Err[iC], Par[iC], N_proc); } else { - propagateHelixToZMPlex(Err[iC], Par[iC], Chg, msPar[hi], - Err[iP], Par[iP], N_proc, useParamBfield); + PropagateTracksToHitZ(msPar[hi], N_proc, pflags); - updateParametersEndcapMPlex(Err[iP], Par[iP], Chg, msErr[hi], msPar[hi], - Err[iC], Par[iC], N_proc); + kalmanUpdateEndcap(Err[iP], Par[iP], msErr[hi], msPar[hi], + Err[iC], Par[iC], N_proc); } if (Config::fit_val) MkFitter::CollectFitValidation(hi,N_proc,ev); diff --git a/mkFit/PropagationMPlex.cc b/mkFit/PropagationMPlex.cc index 2fbf0b599fc9d..87062ea966bef 100644 --- a/mkFit/PropagationMPlex.cc +++ b/mkFit/PropagationMPlex.cc @@ -246,8 +246,8 @@ void MultHelixPropTranspFull(const MPlexLL& A, const MPlexLL& B, MPlexLL& C) #endif } // end unnamed namespace -void helixAtRFromIterativeCCSFullJac(const MPlexLV& inPar, const MPlexQI& inChg, MPlexLV& outPar, - const MPlexQF &msRad, MPlexLL& errorProp, +void helixAtRFromIterativeCCSFullJac(const MPlexLV& inPar, const MPlexQI& inChg, const MPlexQF &msRad, + MPlexLV& outPar, MPlexLL& errorProp, const int N_proc) { errorProp.SetVal(0.f); @@ -377,18 +377,18 @@ void helixAtRFromIterativeCCSFullJac(const MPlexLV& inPar, const MPlexQI& inChg, //#pragma omp declare simd simdlen(NN) notinbranch linear(n) #include "PropagationMPlex.icc" -void helixAtRFromIterativeCCS(const MPlexLV& inPar, const MPlexQI& inChg, MPlexLV& outPar, - const MPlexQF &msRad, MPlexLL& errorProp, - const int N_proc, const bool useParamBfield) +void helixAtRFromIterativeCCS(const MPlexLV& inPar, const MPlexQI& inChg, const MPlexQF &msRad, + MPlexLV& outPar, MPlexLL& errorProp, + const int N_proc, const PropagationFlags pf) { errorProp.SetVal(0.f); - - helixAtRFromIterativeCCS_impl(inPar, inChg, outPar, msRad, errorProp, 0, NN, N_proc, useParamBfield); + + helixAtRFromIterativeCCS_impl(inPar, inChg, msRad, outPar, errorProp, 0, NN, N_proc, pf); } void applyMaterialEffects(const MPlexQF &hitsRl, const MPlexQF& hitsXi, MPlexLS &outErr, MPlexLV& outPar, - const int N_proc) + const int N_proc) { #pragma simd for (int n = 0; n < NN; ++n) @@ -433,22 +433,17 @@ void applyMaterialEffects(const MPlexQF &hitsRl, const MPlexQF& hitsXi, } void propagateHelixToRMPlex(const MPlexLS &inErr, const MPlexLV& inPar, - const MPlexQI &inChg, const MPlexHV& msPar, + const MPlexQI &inChg, const MPlexQF& msRad, MPlexLS &outErr, MPlexLV& outPar, - const int N_proc, const bool useParamBfield) + const int N_proc, const PropagationFlags pf) { outErr = inErr; - outPar = inPar; + outPar = inPar; // XXXX-4G is this requirement for helixAtRFromIterativeCCS_XXX ??? We should document it there. + // Or even better, do it there as it is somewhat counterintuitive :) MPlexLL errorProp; - MPlexQF msRad; -#pragma simd - for (int n = 0; n < NN; ++n) { - msRad.At(n, 0, 0) = hipo(msPar.ConstAt(n, 0, 0), msPar.ConstAt(n, 1, 0)); - } - - helixAtRFromIterativeCCS(inPar, inChg, outPar, msRad, errorProp, N_proc, useParamBfield); + helixAtRFromIterativeCCS(inPar, inChg, msRad, outPar, errorProp, N_proc, pf.use_param_b_field); #ifdef DEBUG { @@ -468,26 +463,29 @@ void propagateHelixToRMPlex(const MPlexLS &inErr, const MPlexLV& inPar, } #endif - if (Config::useCMSGeom && useParamBfield) // useParamBfield is proxy for fittest + if (Config::useCMSGeom && pf.apply_material) { MPlexQF hitsRl; MPlexQF hitsXi; #pragma simd for (int n = 0; n < NN; ++n) { - const int zbin = getZbinME(msPar(n, 2, 0)); - const int rbin = getRbinME(msRad(n, 0, 0)); + const int zbin = getZbinME(outPar(n, 2, 0)); // XXXX-4K changed msPar_z to outPar_z + const int rbin = getRbinME(msRad (n, 0, 0)); hitsRl(n, 0, 0) = (zbin>=0 && zbin<Config::nBinsZME) ? getRlVal(zbin,rbin) : 0.f; // protect against crazy propagations hitsXi(n, 0, 0) = (zbin>=0 && zbin<Config::nBinsZME) ? getXiVal(zbin,rbin) : 0.f; // protect against crazy propagations } - applyMaterialEffects(hitsRl, hitsXi, outErr, outPar, N_proc); + applyMaterialEffects(hitsRl, hitsXi, outErr, outPar, N_proc); // XXXX-4K - why are we doing this on error before propagation ??? + // Isn't the material calculated for this/current layer? } squashPhiMPlex(outPar,N_proc); // ensure phi is between |pi| // Matriplex version of: // result.errors = ROOT::Math::Similarity(errorProp, outErr); + + // MultHelixProp can be optimized for CCS coordinates, see GenMPlexOps.pl MPlexLL temp; MultHelixProp (errorProp, outErr, temp); MultHelixPropTransp(errorProp, temp, outErr); @@ -514,11 +512,21 @@ void propagateHelixToRMPlex(const MPlexLS &inErr, const MPlexLV& inPar, } } #endif + + /* + // To be used with: MPT_DIM = 1 + if (fabs(sqrt(outPar[0]*outPar[0]+outPar[1]*outPar[1])-r)>0.0001) { + std::cout << "DID NOT GET TO R, dR=" << fabs(sqrt(outPar[0]*outPar[0]+outPar[1]*outPar[1])-r) + << " r=" << r << " r0in=" << r0in << " rout=" << sqrt(outPar[0]*outPar[0]+outPar[1]*outPar[1]) << std::endl; + std::cout << "pt=" << pt << " pz=" << inPar.At(n, 2) << std::endl; + } + */ } +/* void propagateHelixToRMPlex(const MPlexLS& inErr, const MPlexLV& inPar, const MPlexQI& inChg, const float r, - MPlexLS& outErr, MPlexLV& outPar, + MPlexLS& outErr, MPlexLV& outPar, const int N_proc) { outErr = inErr; @@ -532,7 +540,7 @@ void propagateHelixToRMPlex(const MPlexLS& inErr, const MPlexLV& inPar, msRad.At(n, 0, 0) = r; } - helixAtRFromIterativeCCS(inPar, inChg, outPar, msRad, errorProp, N_proc); + helixAtRFromIterativeCCS(inPar, inChg, msRad, outPar, errorProp, N_proc); //add multiple scattering uncertainty and energy loss (FIXME: in this way it is not applied in track fit) if (Config::useCMSGeom) @@ -579,33 +587,20 @@ void propagateHelixToRMPlex(const MPlexLS& inErr, const MPlexLV& inPar, } } #endif - - /* - if (fabs(sqrt(outPar[0]*outPar[0]+outPar[1]*outPar[1])-r)>0.0001) { - std::cout << "DID NOT GET TO R, dR=" << fabs(sqrt(outPar[0]*outPar[0]+outPar[1]*outPar[1])-r) - << " r=" << r << " r0in=" << r0in << " rout=" << sqrt(outPar[0]*outPar[0]+outPar[1]*outPar[1]) << std::endl; - std::cout << "pt=" << pt << " pz=" << inPar.At(n, 2) << std::endl; - } - */ } +*/ void propagateHelixToZMPlex(const MPlexLS &inErr, const MPlexLV& inPar, - const MPlexQI &inChg, const MPlexHV& msPar, + const MPlexQI &inChg, const MPlexQF& msZ, MPlexLS &outErr, MPlexLV& outPar, - const int N_proc, const bool useParamBfield) + const int N_proc, const PropagationFlags pf) { outErr = inErr; outPar = inPar; MPlexLL errorProp; - MPlexQF msZ; -#pragma simd - for (int n = 0; n < NN; ++n) { - msZ.At(n, 0, 0) = msPar.ConstAt(n, 2, 0); - } - - helixAtZ(inPar, inChg, outPar, msZ, errorProp, N_proc, useParamBfield); + helixAtZ(inPar, inChg, msZ, outPar, errorProp, N_proc, pf); #ifdef DEBUG { @@ -625,7 +620,7 @@ void propagateHelixToZMPlex(const MPlexLS &inErr, const MPlexLV& inPar, } #endif - if (Config::useCMSGeom && useParamBfield) // param bfield only used in fitting + if (Config::useCMSGeom && pf.apply_material) { MPlexQF hitsRl; MPlexQF hitsXi; @@ -633,7 +628,7 @@ void propagateHelixToZMPlex(const MPlexLS &inErr, const MPlexLV& inPar, for (int n = 0; n < NN; ++n) { const int zbin = getZbinME(msZ(n, 0, 0)); - const int rbin = getRbinME(hipo(msPar(n, 0, 0), msPar(n, 1, 0))); + const int rbin = getRbinME(std::hypot(outPar(n, 0, 0), outPar(n, 1, 0))); // XXXX-4K changed msPar_xy to outPar_xy hitsRl(n, 0, 0) = (rbin>=0 && rbin<Config::nBinsRME) ? getRlVal(zbin,rbin) : 0.f; // protect against crazy propagations hitsXi(n, 0, 0) = (rbin>=0 && rbin<Config::nBinsRME) ? getXiVal(zbin,rbin) : 0.f; // protect against crazy propagations @@ -673,7 +668,7 @@ void propagateHelixToZMPlex(const MPlexLS &inErr, const MPlexLV& inPar, #endif } - +/* void propagateHelixToZMPlex(const MPlexLS &inErr, const MPlexLV& inPar, const MPlexQI &inChg, const float z, MPlexLS &outErr, MPlexLV& outPar, @@ -684,14 +679,6 @@ void propagateHelixToZMPlex(const MPlexLS &inErr, const MPlexLV& inPar, MPlexLL errorProp; - MPlexQF msZ; -#pragma simd - for (int n = 0; n < NN; ++n) { - // XXXXMT4G - // msZ.At(n, 0, 0) = (inPar.ConstAt(n, 2, 0) > 0) ? z : -z; - msZ.At(n, 0, 0) = z; - } - helixAtZ(inPar, inChg, outPar, msZ, errorProp, N_proc); #ifdef DEBUG @@ -762,11 +749,11 @@ void propagateHelixToZMPlex(const MPlexLS &inErr, const MPlexLV& inPar, } #endif } +*/ - -void helixAtZ(const MPlexLV& inPar, const MPlexQI& inChg, MPlexLV& outPar, - const MPlexQF &msZ, MPlexLL& errorProp, - const int N_proc, const bool useParamBfield) +void helixAtZ(const MPlexLV& inPar, const MPlexQI& inChg, const MPlexQF &msZ, + MPlexLV& outPar, MPlexLL& errorProp, + const int N_proc, const PropagationFlags pf) { errorProp.SetVal(0.f); @@ -787,7 +774,7 @@ void helixAtZ(const MPlexLV& inPar, const MPlexQI& inChg, MPlexLV& outPar, const float phiin = inPar.ConstAt(n, 4, 0); const float theta = inPar.ConstAt(n, 5, 0); - const float k = inChg.ConstAt(n, 0, 0) * 100.f / (-Config::sol*(useParamBfield?Config::BfieldFromZR(zin,hipo(inPar.ConstAt(n,0,0),inPar.ConstAt(n,1,0))):Config::Bfield)); + const float k = inChg.ConstAt(n, 0, 0) * 100.f / (-Config::sol*(pf.use_param_b_field?Config::BfieldFromZR(zin,hipo(inPar.ConstAt(n,0,0),inPar.ConstAt(n,1,0))):Config::Bfield)); dprint_np(n, std::endl << "input parameters" << " inPar.ConstAt(n, 0, 0)=" << std::setprecision(9) << inPar.ConstAt(n, 0, 0) diff --git a/mkFit/PropagationMPlex.h b/mkFit/PropagationMPlex.h index cc36fd008f21c..f246db82a71bb 100644 --- a/mkFit/PropagationMPlex.h +++ b/mkFit/PropagationMPlex.h @@ -26,42 +26,41 @@ void propagateLineToRMPlex(const MPlexLS &psErr, const MPlexLV& psPar, const int N_proc); void propagateHelixToRMPlex(const MPlexLS &inErr, const MPlexLV& inPar, - const MPlexQI &inChg, const MPlexHV& msPar, + const MPlexQI &inChg, const MPlexQF& msRad, MPlexLS &outErr, MPlexLV& outPar, - const int N_proc, const bool useParamBfield = false); - + const int N_proc, const PropagationFlags pf); +/* void propagateHelixToRMPlex(const MPlexLS& inErr, const MPlexLV& inPar, const MPlexQI& inChg, const float r, MPlexLS& outErr, MPlexLV& outPar, - const int N_proc); - -void helixAtRFromIterativeCCSFullJac(const MPlexLV& inPar, const MPlexQI& inChg, - MPlexLV& outPar, const MPlexQF &msRad, - MPlexLL& errorProp, + const int N_proc, const PropagationFlags pf); +*/ +void helixAtRFromIterativeCCSFullJac(const MPlexLV& inPar, const MPlexQI& inChg, const MPlexQF &msRad, + MPlexLV& outPar, MPlexLL& errorProp, const int N_proc); -void helixAtRFromIterativeCCS(const MPlexLV& inPar, const MPlexQI& inChg, - MPlexLV& outPar, const MPlexQF &msRad, - MPlexLL& errorProp, - const int N_proc, const bool useParamBfield = false); +void helixAtRFromIterativeCCS(const MPlexLV& inPar, const MPlexQI& inChg, const MPlexQF &msRad, + MPlexLV& outPar, MPlexLL& errorProp, + const int N_proc, const PropagationFlags pf); void propagateHelixToZMPlex(const MPlexLS &inErr, const MPlexLV& inPar, - const MPlexQI &inChg, const MPlexHV& msPar, + const MPlexQI &inChg, const MPlexQF& msZ, MPlexLS &outErr, MPlexLV& outPar, - const int N_proc, const bool useParamBfield = false); + const int N_proc, const PropagationFlags pf); +/* void propagateHelixToZMPlex(const MPlexLS &inErr, const MPlexLV& inPar, const MPlexQI &inChg, const float z, MPlexLS &outErr, MPlexLV& outPar, const int N_proc); +*/ -void helixAtZ(const MPlexLV& inPar, const MPlexQI& inChg, - MPlexLV& outPar, const MPlexQF &msZ, - MPlexLL& errorProp, - const int N_proc, const bool useParamBfield = false); +void helixAtZ(const MPlexLV& inPar, const MPlexQI& inChg, const MPlexQF &msZ, + MPlexLV& outPar, MPlexLL& errorProp, + const int N_proc, const PropagationFlags pf); void applyMaterialEffects(const MPlexQF &hitsRl, const MPlexQF& hitsXi, - MPlexLS &outErr, MPlexLV& outPar, + MPlexLS &outErr, MPlexLV& outPar, const int N_proc); #endif diff --git a/mkFit/PropagationMPlex.icc b/mkFit/PropagationMPlex.icc index 946dae9527d0d..eefcc8c90b2ab 100644 --- a/mkFit/PropagationMPlex.icc +++ b/mkFit/PropagationMPlex.icc @@ -6,14 +6,14 @@ template<typename Tf, typename Ti, typename TfLL1, typename Tf11, typename TfLLL #ifdef __CUDACC__ __device__ #endif -static inline void helixAtRFromIterativeCCS_impl(const Tf& __restrict__ inPar, - const Ti& __restrict__ inChg, - TfLL1& __restrict__ outPar, - const Tf11& __restrict__ msRad, - TfLLL& __restrict__ errorProp, +static inline void helixAtRFromIterativeCCS_impl(const Tf& __restrict__ inPar, + const Ti& __restrict__ inChg, + const Tf11& __restrict__ msRad, + TfLL1& __restrict__ outPar, + TfLLL& __restrict__ errorProp, const int nmin, const int nmax, const int N_proc, - const bool useParamBfield) + const PropagationFlags pf) { #pragma simd for (int n = nmin; n < nmax; ++n) @@ -27,7 +27,7 @@ static inline void helixAtRFromIterativeCCS_impl(const Tf& __restrict__ inPar, errorProp(n,5,5) = 1.f; float r0 = hipo(inPar(n, 0, 0), inPar(n, 1, 0)); - const float k = inChg(n, 0, 0) * 100.f / (-Config::sol*(useParamBfield?Config::BfieldFromZR(inPar(n,2,0),r0):Config::Bfield)); + const float k = inChg(n, 0, 0) * 100.f / (-Config::sol*(pf.use_param_b_field ? Config::BfieldFromZR(inPar(n,2,0),r0) : Config::Bfield)); const float r = msRad(n, 0, 0); // if (std::abs(r-r0)<0.0001f) { diff --git a/mkFit/SteeringParams.h b/mkFit/SteeringParams.h index d46f7e9c05ff4..a296d79add10d 100644 --- a/mkFit/SteeringParams.h +++ b/mkFit/SteeringParams.h @@ -11,23 +11,23 @@ class MkFinder; #define COMPUTE_CHI2_ARGS const MPlexLS &, const MPlexLV &, const MPlexQI &, \ const MPlexHS &, const MPlexHV &, \ - MPlexQF &, const int + MPlexQF &, const int, const PropagationFlags #define UPDATE_PARAM_ARGS const MPlexLS &, const MPlexLV &, const MPlexQI &, \ const MPlexHS &, const MPlexHV &, \ - MPlexLS &, MPlexLV &, const int + MPlexLS &, MPlexLV &, const int, const PropagationFlags struct FindingFoos { void (*m_compute_chi2_foo) (COMPUTE_CHI2_ARGS); void (*m_update_param_foo) (UPDATE_PARAM_ARGS); - void (MkBase::*m_propagate_foo) (float, const int); + void (MkBase::*m_propagate_foo) (float, const int, const PropagationFlags); FindingFoos() {} FindingFoos(void (*cch2_f) (COMPUTE_CHI2_ARGS), void (*updp_f) (UPDATE_PARAM_ARGS), - void (MkBase::*p_f) (float, const int)) : + void (MkBase::*p_f) (float, const int, const PropagationFlags)) : m_compute_chi2_foo(cch2_f), m_update_param_foo(updp_f), m_propagate_foo(p_f) From c40e9d610dfeb760f257e4652c2a0662f73579c6 Mon Sep 17 00:00:00 2001 From: Matevz Tadel <mtadel@ucsd.edu> Date: Thu, 21 Dec 2017 16:36:05 -0800 Subject: [PATCH 02/11] Remove old dead code, move applyMaterialEffects() to the bottom (as is in .h). --- mkFit/PropagationMPlex.cc | 243 ++++++++------------------------------ mkFit/PropagationMPlex.h | 14 +-- 2 files changed, 52 insertions(+), 205 deletions(-) diff --git a/mkFit/PropagationMPlex.cc b/mkFit/PropagationMPlex.cc index 87062ea966bef..c93529ed6ce52 100644 --- a/mkFit/PropagationMPlex.cc +++ b/mkFit/PropagationMPlex.cc @@ -246,6 +246,9 @@ void MultHelixPropTranspFull(const MPlexLL& A, const MPlexLL& B, MPlexLL& C) #endif } // end unnamed namespace + +//============================================================================== + void helixAtRFromIterativeCCSFullJac(const MPlexLV& inPar, const MPlexQI& inChg, const MPlexQF &msRad, MPlexLV& outPar, MPlexLL& errorProp, const int N_proc) @@ -386,51 +389,6 @@ void helixAtRFromIterativeCCS(const MPlexLV& inPar, const MPlexQI& inChg, co helixAtRFromIterativeCCS_impl(inPar, inChg, msRad, outPar, errorProp, 0, NN, N_proc, pf); } -void applyMaterialEffects(const MPlexQF &hitsRl, const MPlexQF& hitsXi, - MPlexLS &outErr, MPlexLV& outPar, - const int N_proc) -{ -#pragma simd - for (int n = 0; n < NN; ++n) - { - float radL = hitsRl.ConstAt(n,0,0); - if (radL<0.0000000000001f) continue;//ugly, please fixme - const float theta = outPar.ConstAt(n,0,5); - const float pt = 1.f/outPar.ConstAt(n,0,3); - const float p = pt/std::sin(theta); - const float p2 = p*p; - constexpr float mpi = 0.140; // m=140 MeV, pion - constexpr float mpi2 = mpi*mpi; // m=140 MeV, pion - const float beta2 = p2/(p2+mpi2); - const float beta = std::sqrt(beta2); - //radiation lenght, corrected for the crossing angle (cos alpha from dot product of radius vector and momentum) - const float invCos = p/pt; - radL = radL * invCos; //fixme works only for barrel geom - // multiple scattering - //vary independently phi and theta by the rms of the planar multiple scattering angle - const float thetaMSC = 0.0136f*std::sqrt(radL)*(1.f+0.038f*std::log(radL))/(beta*p);// eq 32.15 - const float thetaMSC2 = thetaMSC*thetaMSC; - outErr.At(n, 4, 4) += thetaMSC2; - outErr.At(n, 5, 5) += thetaMSC2; - //std::cout << "beta=" << beta << " p=" << p << std::endl; - //std::cout << "multiple scattering thetaMSC=" << thetaMSC << " thetaMSC2=" << thetaMSC2 << " radL=" << radL << std::endl; - // energy loss - const float gamma = 1.f/std::sqrt(1.f - beta2); - const float gamma2 = gamma*gamma; - constexpr float me = 0.0005; // m=0.5 MeV, electron - const float wmax = 2.f*me*beta2*gamma2 / ( 1.f + 2.f*gamma*me/mpi + me*me/(mpi*mpi) ); - constexpr float I = 16.0e-9 * 10.75; - const float deltahalf = std::log(28.816e-9f * std::sqrt(2.33f*0.498f)/I) + std::log(beta*gamma) - 0.5f; - const float dEdx = beta<1.f ? (2.f*(hitsXi.ConstAt(n,0,0) * invCos * (0.5f*std::log(2.f*me*beta2*gamma2*wmax/(I*I)) - beta2 - deltahalf) / beta2)) : 0.f;//protect against infs and nans - // dEdx = dEdx*2.;//xi in cmssw is defined with an extra factor 0.5 with respect to formula 27.1 in pdg - //std::cout << "dEdx=" << dEdx << " delta=" << deltahalf << " wmax=" << wmax << " Xi=" << hitsXi.ConstAt(n,0,0) << std::endl; - const float dP = dEdx/beta; - outPar.At(n, 0, 3) = p/((p+dP)*pt); - //assume 100% uncertainty - outErr.At(n, 3, 3) += dP*dP/(p2*pt*pt); - } - -} void propagateHelixToRMPlex(const MPlexLS &inErr, const MPlexLV& inPar, const MPlexQI &inChg, const MPlexQF& msRad, @@ -523,72 +481,8 @@ void propagateHelixToRMPlex(const MPlexLS &inErr, const MPlexLV& inPar, */ } -/* -void propagateHelixToRMPlex(const MPlexLS& inErr, const MPlexLV& inPar, - const MPlexQI& inChg, const float r, - MPlexLS& outErr, MPlexLV& outPar, - const int N_proc) -{ - outErr = inErr; - outPar = inPar; - - MPlexLL errorProp; - - MPlexQF msRad; -#pragma simd - for (int n = 0; n < NN; ++n) { - msRad.At(n, 0, 0) = r; - } - - helixAtRFromIterativeCCS(inPar, inChg, msRad, outPar, errorProp, N_proc); - - //add multiple scattering uncertainty and energy loss (FIXME: in this way it is not applied in track fit) - if (Config::useCMSGeom) - { - MPlexQF hitsRl; - MPlexQF hitsXi; - -#pragma simd - for (int n = 0; n < N_proc; ++n) - { - const int zbin = getZbinME(outPar(n, 2, 0)); - const int rbin = getRbinME(r); - - hitsRl(n, 0, 0) = (zbin>=0 && zbin<Config::nBinsZME) ? getRlVal(zbin,rbin) : 0.f; // protect against crazy propagations - hitsXi(n, 0, 0) = (zbin>=0 && zbin<Config::nBinsZME) ? getXiVal(zbin,rbin) : 0.f; // protect against crazy propagations - } - applyMaterialEffects(hitsRl, hitsXi, outErr, outPar, N_proc); - } - - squashPhiMPlex(outPar,N_proc); // ensure phi is between |pi| - - // Matriplex version of: - // result.errors = ROOT::Math::Similarity(errorProp, outErr); - - //MultHelixProp can be optimized for CCS coordinates, see GenMPlexOps.pl - MPlexLL temp; - MultHelixProp (errorProp, outErr, temp); - MultHelixPropTransp(errorProp, temp, outErr); - - // This dump is now out of its place as similarity is done with matriplex ops. -#ifdef DEBUG - if (debug) { - for (int kk = 0; kk < N_proc; ++kk) - { - dprintf("outErr %d\n", kk); - for (int i = 0; i < 6; ++i) { for (int j = 0; j < 6; ++j) - dprintf("%8f ", outErr.At(kk,i,j)); printf("\n"); - } dprintf("\n"); - dprintf("outPar %d\n", kk); - for (int i = 0; i < 6; ++i) { - dprintf("%8f ", outPar.At(kk,i,0)); printf("\n"); - } dprintf("\n"); - } - } -#endif -} -*/ +//============================================================================== void propagateHelixToZMPlex(const MPlexLS &inErr, const MPlexLV& inPar, const MPlexQI &inChg, const MPlexQF& msZ, @@ -668,88 +562,6 @@ void propagateHelixToZMPlex(const MPlexLS &inErr, const MPlexLV& inPar, #endif } -/* -void propagateHelixToZMPlex(const MPlexLS &inErr, const MPlexLV& inPar, - const MPlexQI &inChg, const float z, - MPlexLS &outErr, MPlexLV& outPar, - const int N_proc) -{ - outErr = inErr; - outPar = inPar; - - MPlexLL errorProp; - - helixAtZ(inPar, inChg, outPar, msZ, errorProp, N_proc); - -#ifdef DEBUG - { - for (int kk = 0; kk < N_proc; ++kk) - { - dprintf("inErr %d\n", kk); - for (int i = 0; i < 6; ++i) { for (int j = 0; j < 6; ++j) - dprintf("%8f ", inErr.ConstAt(kk,i,j)); printf("\n"); - } dprintf("\n"); - - dprintf("errorProp %d\n", kk); - for (int i = 0; i < 6; ++i) { for (int j = 0; j < 6; ++j) - dprintf("%8f ", errorProp.At(kk,i,j)); printf("\n"); - } dprintf("\n"); - - } - } -#endif - - if (Config::useCMSGeom) - { - MPlexQF hitsRl; - MPlexQF hitsXi; -#pragma simd - for (int n = 0; n < N_proc; ++n) - { - const int zbin = getZbinME(z); - const int rbin = getRbinME(hipo(outPar(n, 0, 0), outPar(n, 1, 0))); - - hitsRl(n, 0, 0) = (rbin>=0 && rbin<Config::nBinsRME) ? getRlVal(zbin,rbin) : 0.f; // protect against crazy propagations - hitsXi(n, 0, 0) = (rbin>=0 && rbin<Config::nBinsRME) ? getXiVal(zbin,rbin) : 0.f; // protect against crazy propagations - } - applyMaterialEffects(hitsRl, hitsXi, outErr, outPar, N_proc); - } - - squashPhiMPlex(outPar,N_proc); // ensure phi is between |pi| - - // Matriplex version of: - // result.errors = ROOT::Math::Similarity(errorProp, outErr); - MPlexLL temp; - MultHelixPropEndcap (errorProp, outErr, temp); - MultHelixPropTranspEndcap(errorProp, temp, outErr); - - // This dump is now out of its place as similarity is done with matriplex ops. -#ifdef DEBUG - { - for (int kk = 0; kk < N_proc; ++kk) - { - dprintf("outErr %d\n", kk); - for (int i = 0; i < 6; ++i) { for (int j = 0; j < 6; ++j) - dprintf("%8f ", outErr.At(kk,i,j)); printf("\n"); - } dprintf("\n"); - - dprintf("outPar %d\n", kk); - for (int i = 0; i < 6; ++i) { - dprintf("%8f ", outPar.At(kk,i,0)); printf("\n"); - } dprintf("\n"); - // XXXXMT4G - // if (std::abs(outPar.At(kk,2,0) - inPar.ConstAt(kk, 2, 0))>0.0001) { - // dprint_np(kk, "DID NOT GET TO Z, dZ=" << std::abs(outPar.At(kk,2,0)-inPar.ConstAt(kk, 2, 0)) - if (std::abs(outPar.At(kk,2,0) - z) > 0.0001) { - dprint_np(kk, "DID NOT GET TO Z, dZ=" << std::abs(outPar.At(kk,2,0) - z) - << " z=" << msZ.ConstAt(kk, 2, 0) << " zin=" << inPar.ConstAt(kk,2,0) << " zout=" << outPar.At(kk,2,0) - << "\n pt=" << hipo(inPar.ConstAt(kk,3,0), inPar.ConstAt(kk,4,0)) << " pz=" << inPar.ConstAt(kk,5,0)); - } - } - } -#endif -} -*/ void helixAtZ(const MPlexLV& inPar, const MPlexQI& inChg, const MPlexQF &msZ, MPlexLV& outPar, MPlexLL& errorProp, @@ -851,3 +663,50 @@ void helixAtZ(const MPlexLV& inPar, const MPlexQI& inChg, const MPlexQF &msZ, #endif } } + +//============================================================================== + +void applyMaterialEffects(const MPlexQF &hitsRl, const MPlexQF& hitsXi, + MPlexLS &outErr, MPlexLV& outPar, + const int N_proc) +{ +#pragma simd + for (int n = 0; n < NN; ++n) + { + float radL = hitsRl.ConstAt(n,0,0); + if (radL<0.0000000000001f) continue;//ugly, please fixme + const float theta = outPar.ConstAt(n,0,5); + const float pt = 1.f/outPar.ConstAt(n,0,3); + const float p = pt/std::sin(theta); + const float p2 = p*p; + constexpr float mpi = 0.140; // m=140 MeV, pion + constexpr float mpi2 = mpi*mpi; // m=140 MeV, pion + const float beta2 = p2/(p2+mpi2); + const float beta = std::sqrt(beta2); + //radiation lenght, corrected for the crossing angle (cos alpha from dot product of radius vector and momentum) + const float invCos = p/pt; + radL = radL * invCos; //fixme works only for barrel geom + // multiple scattering + //vary independently phi and theta by the rms of the planar multiple scattering angle + const float thetaMSC = 0.0136f*std::sqrt(radL)*(1.f+0.038f*std::log(radL))/(beta*p);// eq 32.15 + const float thetaMSC2 = thetaMSC*thetaMSC; + outErr.At(n, 4, 4) += thetaMSC2; + outErr.At(n, 5, 5) += thetaMSC2; + //std::cout << "beta=" << beta << " p=" << p << std::endl; + //std::cout << "multiple scattering thetaMSC=" << thetaMSC << " thetaMSC2=" << thetaMSC2 << " radL=" << radL << std::endl; + // energy loss + const float gamma = 1.f/std::sqrt(1.f - beta2); + const float gamma2 = gamma*gamma; + constexpr float me = 0.0005; // m=0.5 MeV, electron + const float wmax = 2.f*me*beta2*gamma2 / ( 1.f + 2.f*gamma*me/mpi + me*me/(mpi*mpi) ); + constexpr float I = 16.0e-9 * 10.75; + const float deltahalf = std::log(28.816e-9f * std::sqrt(2.33f*0.498f)/I) + std::log(beta*gamma) - 0.5f; + const float dEdx = beta<1.f ? (2.f*(hitsXi.ConstAt(n,0,0) * invCos * (0.5f*std::log(2.f*me*beta2*gamma2*wmax/(I*I)) - beta2 - deltahalf) / beta2)) : 0.f;//protect against infs and nans + // dEdx = dEdx*2.;//xi in cmssw is defined with an extra factor 0.5 with respect to formula 27.1 in pdg + //std::cout << "dEdx=" << dEdx << " delta=" << deltahalf << " wmax=" << wmax << " Xi=" << hitsXi.ConstAt(n,0,0) << std::endl; + const float dP = dEdx/beta; + outPar.At(n, 0, 3) = p/((p+dP)*pt); + //assume 100% uncertainty + outErr.At(n, 3, 3) += dP*dP/(p2*pt*pt); + } +} diff --git a/mkFit/PropagationMPlex.h b/mkFit/PropagationMPlex.h index f246db82a71bb..394fd1a9f986c 100644 --- a/mkFit/PropagationMPlex.h +++ b/mkFit/PropagationMPlex.h @@ -29,12 +29,7 @@ void propagateHelixToRMPlex(const MPlexLS &inErr, const MPlexLV& inPar, const MPlexQI &inChg, const MPlexQF& msRad, MPlexLS &outErr, MPlexLV& outPar, const int N_proc, const PropagationFlags pf); -/* -void propagateHelixToRMPlex(const MPlexLS& inErr, const MPlexLV& inPar, - const MPlexQI& inChg, const float r, - MPlexLS& outErr, MPlexLV& outPar, - const int N_proc, const PropagationFlags pf); -*/ + void helixAtRFromIterativeCCSFullJac(const MPlexLV& inPar, const MPlexQI& inChg, const MPlexQF &msRad, MPlexLV& outPar, MPlexLL& errorProp, const int N_proc); @@ -48,13 +43,6 @@ void propagateHelixToZMPlex(const MPlexLS &inErr, const MPlexLV& inPar, MPlexLS &outErr, MPlexLV& outPar, const int N_proc, const PropagationFlags pf); -/* -void propagateHelixToZMPlex(const MPlexLS &inErr, const MPlexLV& inPar, - const MPlexQI &inChg, const float z, - MPlexLS &outErr, MPlexLV& outPar, - const int N_proc); -*/ - void helixAtZ(const MPlexLV& inPar, const MPlexQI& inChg, const MPlexQF &msZ, MPlexLV& outPar, MPlexLL& errorProp, const int N_proc, const PropagationFlags pf); From a2743177534a5cba292aafbcbbd11a09e0ccbbef Mon Sep 17 00:00:00 2001 From: Matevz Tadel <mtadel@ucsd.edu> Date: Fri, 22 Dec 2017 12:01:40 -0800 Subject: [PATCH 03/11] Make Finder-FV use PropagationFlags. --- mkFit/MkBuilder.cc | 2 +- mkFit/MkFinderFV.cc | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/mkFit/MkBuilder.cc b/mkFit/MkBuilder.cc index 6d16a87da4fa4..445431a69e551 100644 --- a/mkFit/MkBuilder.cc +++ b/mkFit/MkBuilder.cc @@ -1812,7 +1812,7 @@ void MkBuilder::find_tracks_in_layersFV(int start_seed, int end_seed, int region auto& mkfndr = finders[index]; // propagate to current layer - (mkfndr.*fnd_foos.m_propagate_foo)(layer_info.m_propagate_to, mkfndr.nnfv()); + (mkfndr.*fnd_foos.m_propagate_foo)(layer_info.m_propagate_to, mkfndr.nnfv(), PF_use_param_b_field); mkfndr.SelectHitIndices(layer_of_hits); diff --git a/mkFit/MkFinderFV.cc b/mkFit/MkFinderFV.cc index b30a126572f7d..348e91895c059 100644 --- a/mkFit/MkFinderFV.cc +++ b/mkFit/MkFinderFV.cc @@ -318,7 +318,7 @@ void MkFinderFV<nseeds, ncands>::FindCandidates(const LayerOfHits &layer_of_hits #endif //now compute the chi2 of track state vs hit - (*fnd_foos.m_compute_chi2_foo)(Err[iP], Par[iP], Chg, msErr, msPar, XHitChi2[hit_cnt], NNFV); + (*fnd_foos.m_compute_chi2_foo)(Err[iP], Par[iP], Chg, msErr, msPar, XHitChi2[hit_cnt], NNFV, PF_use_param_b_field); // Prefetch to L1 the hits we'll (probably) process in the next loop iteration. for (int itrack = 0; itrack < NNFV; ++itrack) @@ -357,7 +357,7 @@ void MkFinderFV<nseeds, ncands>::UpdateWithLastHit(const LayerOfHits &layer_of_h } (*fnd_foos.m_update_param_foo)(Err[iP], Par[iP], Chg, msErr, msPar, - Err[iC], Par[iC], NNFV); + Err[iC], Par[iC], NNFV, PF_use_param_b_field); //now that we have moved propagation at the end of the sequence we lost the handle of //using the propagated parameters instead of the updated for the missing hit case. From 23fd2a47fe4fb1709ad415ee183ea4d269d024ef Mon Sep 17 00:00:00 2001 From: Matevz Tadel <mtadel@ucsd.edu> Date: Thu, 4 Jan 2018 13:28:38 -0800 Subject: [PATCH 04/11] Use PropagationFlags also in Propagation and Simulation. --- Propagation.cc | 30 ++++++++++++++++-------------- Propagation.h | 8 ++++---- Simulation.cc | 4 +++- buildtest.cc | 8 +++++--- seedtest.cc | 22 ++++++++++++++++------ 5 files changed, 44 insertions(+), 28 deletions(-) diff --git a/Propagation.cc b/Propagation.cc index 126ee699cf150..62d2a86952ef3 100644 --- a/Propagation.cc +++ b/Propagation.cc @@ -234,13 +234,13 @@ void HelixState::propagateErrors(const HelixState& in, float totalDistance, bool // each step travels for a path length equal to the safe step between the current position and the nearest object. TrackState propagateHelixToNextSolid(TrackState inputState, const Geometry& geom, - const bool useParamBfield) + const PropagationFlags pflags) { bool debug = true; - const HelixState hsin(inputState,useParamBfield); + const HelixState hsin(inputState, pflags.use_param_b_field); TrackState result(inputState); - HelixState hsout(result,useParamBfield); + HelixState hsout(result, pflags.use_param_b_field); #ifdef CHECKSTATEVALID if (!hsout.state.valid) { @@ -341,17 +341,19 @@ TrackState propagateHelixToNextSolid(TrackState inputState, const Geometry& geom // XXMT4K The following is only used in buildtest.cc, maybe obsolete? -// Propagate to the next obj -// each step travels for a path length equal to the safe step between the current position and the nearest object. -TrackState propagateHelixToLayer(TrackState inputState, int layer, const Geometry& geom, const bool useParamBfield) +// Propagate to the next obj. Eeach step travels for a path length equal to +// the safe step between the current position and the nearest object. + +TrackState propagateHelixToLayer(TrackState inputState, int layer, const Geometry& geom, + const PropagationFlags pflags) { bool debug = true; const VUSolid* target = geom.Layer(layer); - const HelixState hsin(inputState,useParamBfield); + const HelixState hsin(inputState, pflags.use_param_b_field); TrackState result(inputState); - HelixState hsout(result,useParamBfield); + HelixState hsout(result, pflags.use_param_b_field); #ifdef CHECKSTATEVALID if (!hsout.state.valid) { @@ -405,13 +407,13 @@ TrackState propagateHelixToLayer(TrackState inputState, int layer, const Geometr // each step travels for a path lenght equal to delta r between the current position and the target radius. // for track with pT>=1 GeV this converges to the correct path lenght in <5 iterations // derivatives need to be updated at each iteration -TrackState propagateHelixToR(TrackState inputState, float r, const bool useParamBfield) +TrackState propagateHelixToR(TrackState inputState, float r, const PropagationFlags pflags) { bool debug = false; - const HelixState hsin(inputState,useParamBfield); + const HelixState hsin(inputState, pflags.use_param_b_field); TrackState result(inputState); - HelixState hsout(result,useParamBfield); + HelixState hsout(result, pflags.use_param_b_field); #ifdef CHECKSTATEVALID if (!hsout.state.valid) { @@ -470,8 +472,8 @@ TrackState propagateHelixToR(TrackState inputState, float r, const bool useParam return hsout.state; } -TrackState propagateHelixToZ(TrackState inputState, float zout, const bool useParamBfield) { - +TrackState propagateHelixToZ(TrackState inputState, float zout, const PropagationFlags pflags) +{ TrackState result = inputState; const float z = inputState.z(); @@ -485,7 +487,7 @@ TrackState propagateHelixToZ(TrackState inputState, float zout, const bool usePa const float cosP = std::cos(phi); const float sinP = std::sin(phi); - const float k = inputState.charge*100.f/(-Config::sol*(useParamBfield?Config::BfieldFromZR(z,inputState.posR()):Config::Bfield)); + const float k = inputState.charge*100.f/(-Config::sol*(pflags.use_param_b_field ? Config::BfieldFromZR(z,inputState.posR()) : Config::Bfield)); const float s = (zout - z)/cosT; const float angPath = s*sinT*ipT/k; diff --git a/Propagation.h b/Propagation.h index d0ba55f0b41d8..1d35feccb9949 100644 --- a/Propagation.h +++ b/Propagation.h @@ -8,16 +8,16 @@ // assuming radial direction (i.e. origin at (0,0)) TrackState propagateLineToR(const TrackState& inputState, float r); -TrackState propagateHelixToNextSolid(TrackState inputState, const Geometry& geom, const bool useParamBfield = false); -TrackState propagateHelixToLayer(TrackState inputState, int layer, const Geometry& geom, const bool useParamBfield = false); +TrackState propagateHelixToNextSolid(TrackState inputState, const Geometry& geom, const PropagationFlags pflags); +TrackState propagateHelixToLayer(TrackState inputState, int layer, const Geometry& geom, const PropagationFlags pflags); // helix propagation in steps along helix trajectory. // each step travels for a path lenght equal to delta r between the current position and the target radius. // for track with pT>=1 GeV this converges to the correct path lenght in <5 iterations // derivatives need to be updated at each iteration -TrackState propagateHelixToR(TrackState inputState, float r, const bool useParamBfield = false); +TrackState propagateHelixToR(TrackState inputState, float r, const PropagationFlags pflags); // -TrackState propagateHelixToZ(TrackState inputState, float z, const bool useParamBfield = false); +TrackState propagateHelixToZ(TrackState inputState, float z, const PropagationFlags pflags); #endif diff --git a/Simulation.cc b/Simulation.cc index 2f052b0b56ca2..8ed192efd98d8 100644 --- a/Simulation.cc +++ b/Simulation.cc @@ -89,6 +89,8 @@ void setupTrackByToyMC(SVector3& pos, SVector3& mom, SMatrixSym66& covtrk, // XXMT4M - This should become while not in outer layer (or at least propState.state == invalid) // Try the invalid thingy first ... but would be good to also know what layers are final. + const PropagationFlags pflags(PF_use_param_b_field); + for (int ihit = 0; ihit < Config::nTotHit; ++ihit) { dprintf("\n================================================================================\n"); @@ -96,7 +98,7 @@ void setupTrackByToyMC(SVector3& pos, SVector3& mom, SMatrixSym66& covtrk, ihit, simLayer, tmpState.x(), tmpState.y(), tmpState.posR(), tmpState.z(), tmpState.posPhi()); dprintf("--------------------------------------------------------------------------------\n"); - auto propState = propagateHelixToNextSolid(tmpState,geom,true); + auto propState = propagateHelixToNextSolid(tmpState, geom, pflags); float initX = propState.parameters.At(0); float initY = propState.parameters.At(1); diff --git a/buildtest.cc b/buildtest.cc index 9d07839d644c4..ad2a9a878992a 100644 --- a/buildtest.cc +++ b/buildtest.cc @@ -176,14 +176,16 @@ void extendCandidate(const BinInfoMap & segmentMap, const Event& ev, const cand_ const auto& evt_lay_hits(ev.layerHits_); const auto& segLayMap(segmentMap[ilayer]); + const PropagationFlags pflags(PF_none); + dprint("processing candidate with nHits=" << tkcand.nFoundHits()); #ifdef LINEARINTERP - TrackState propState = propagateHelixToR(updatedState,ev.geom_.Radius(ilayer)); + TrackState propState = propagateHelixToR(updatedState, ev.geom_.Radius(ilayer), pflags); #else #ifdef TBB #error "Invalid combination of options (thread safety)" #endif - TrackState propState = propagateHelixToLayer(updatedState,ilayer,ev.geom_); + TrackState propState = propagateHelixToLayer(updatedState, ilayer,ev.geom_, pflags); #endif // LINEARINTERP #ifdef CHECKSTATEVALID if (!propState.valid) { @@ -221,7 +223,7 @@ void extendCandidate(const BinInfoMap & segmentMap, const Event& ev, const cand_ const float deltaR = maxR - minR; dprint("min, max, delta: " << minR << ", " << maxR << ", " << deltaR); const TrackState propStateMin = propState; - const TrackState propStateMax = propagateHelixToR(updatedState,maxR); + const TrackState propStateMax = propagateHelixToR(updatedState, maxR, pflags); #ifdef CHECKSTATEVALID if (!propStateMax.valid) { return; diff --git a/seedtest.cc b/seedtest.cc index d9b0d30dfc5c0..ad327993e44cd 100644 --- a/seedtest.cc +++ b/seedtest.cc @@ -20,7 +20,10 @@ void buildSeedsByMC(const TrackVec& evt_sim_tracks, TrackVec& evt_seed_tracks, T bool debug(true); #endif - for (int itrack=0;itrack<evt_sim_tracks.size();++itrack) { + const PropagationFlags pflags(PF_none); + + for (int itrack=0;itrack<evt_sim_tracks.size();++itrack) + { const Track& trk = evt_sim_tracks[itrack]; int seedhits[Config::nLayers]; //float sumchi2 = 0; @@ -37,10 +40,12 @@ void buildSeedsByMC(const TrackVec& evt_sim_tracks, TrackVec& evt_seed_tracks, T dprint("processing sim track # " << itrack << " par=" << trk.parameters()); TSLayerPairVec updatedStates; // validation for position pulls - for (auto ilayer=0;ilayer<Config::nlayers_per_seed;++ilayer) {//seeds have first three layers as seeds + // use first layers as seed hits + for (auto ilayer=0;ilayer<Config::nlayers_per_seed;++ilayer) + { auto hitidx = trk.getHitIdx(ilayer); const Hit& seed_hit = ev.layerHits_[ilayer][hitidx]; - TrackState propState = propagateHelixToR(updatedState,seed_hit.r()); + TrackState propState = propagateHelixToR(updatedState, seed_hit.r(), pflags); #ifdef CHECKSTATEVALID if (!propState.valid) { std::cout << "Seeding failed to propagate to layer: " << ilayer << " for sim track: " << itrack << std::endl; @@ -527,6 +532,8 @@ void buildSeedsFromTriplets(const std::vector<HitVec>& evt_lay_hits, const Tripl fprintf(stderr, "__FILE__::__LINE__ Needs fixing for B/E support, search for XXMT4K\n"); exit(1); + const PropagationFlags pflags(PF_none); + // now perform kalman fit on seeds --> first need initial parameters --> get from Conformal fitter! unsigned int seedID = 0; for(auto&& hit_triplet : filtered_triplets){ @@ -548,7 +555,7 @@ void buildSeedsFromTriplets(const std::vector<HitVec>& evt_lay_hits, const Tripl for (auto ilayer = 0; ilayer < Config::nlayers_per_seed; ++ilayer) { const Hit& seed_hit = evt_lay_hits[ilayer][hit_triplet[ilayer]]; - TrackState propState = propagateHelixToR(updatedState,seed_hit.r()); + TrackState propState = propagateHelixToR(updatedState, seed_hit.r(), pflags); #ifdef CHECKSTATEVALID if (!propState.valid) { break; @@ -573,7 +580,10 @@ void buildSeedsFromTriplets(const std::vector<HitVec>& evt_lay_hits, const Tripl void fitSeeds(const std::vector<HitVec>& evt_lay_hits, TrackVec& evt_seed_tracks, Event& ev){ // copy+paste code to fit needs here... - for(auto&& seed : evt_seed_tracks) { + const PropagationFlags pflags(PF_none); + + for (auto&& seed : evt_seed_tracks) + { // state to save to track TrackState updatedState; //const int seedID = seed.label(); // label == seedID! @@ -586,7 +596,7 @@ void fitSeeds(const std::vector<HitVec>& evt_lay_hits, TrackVec& evt_seed_tracks float sumchi2 = 0.0f; // chi2 to save for (auto ilayer = 0; ilayer < Config::nlayers_per_seed; ++ilayer) { const Hit& seed_hit = evt_lay_hits[ilayer][seed.getHitIdx(ilayer)]; - TrackState propState = propagateHelixToR(updatedState,seed_hit.r()); + TrackState propState = propagateHelixToR(updatedState, seed_hit.r(), pflags); #ifdef CHECKSTATEVALID if (!propState.valid) { break; From 74acf9fd50279c9c932b6876d0b693710517dbfd Mon Sep 17 00:00:00 2001 From: Matevz Tadel <mtadel@ucsd.edu> Date: Thu, 4 Jan 2018 15:33:33 -0800 Subject: [PATCH 05/11] Add 'functional' PropagationFlags into Config, define them for each geometry in the plugin, and use them in code. --- Config.cc | 7 +++ Config.h | 14 ++++-- Geoms/CMS-2017.cc | 8 ++++ Geoms/CylCowWLids.cc | 7 +++ mkFit/KalmanUtilsMPlex.cc | 96 +++++++++++++++++++-------------------- mkFit/MkBuilder.cc | 11 ++--- mkFit/MkFinder.cc | 26 ++++++----- mkFit/MkFinderFV.cc | 4 +- mkFit/MkFitter.cc | 10 ++-- mkFit/MkFitter.h | 4 +- mkFit/PropagationMPlex.cc | 20 ++++---- mkFit/PropagationMPlex.h | 8 ++-- 12 files changed, 122 insertions(+), 93 deletions(-) diff --git a/Config.cc b/Config.cc index 3d52122a622b9..1ba4918c3570e 100644 --- a/Config.cc +++ b/Config.cc @@ -48,6 +48,13 @@ namespace Config seedOpts seedInput = simSeeds; cleanOpts seedCleaning = noCleaning; + bool finding_requires_propagation_to_hit_pos; + PropagationFlags finding_inter_layer_pflags; + PropagationFlags finding_intra_layer_pflags; + PropagationFlags backward_fit_pflags; + PropagationFlags forward_fit_pflags; + PropagationFlags seed_fit_pflags; + bool useCMSGeom = false; bool readCmsswTracks = false; diff --git a/Config.h b/Config.h index 0443b67e9a054..61f8dacf355ad 100644 --- a/Config.h +++ b/Config.h @@ -221,18 +221,26 @@ namespace Config constexpr float seed_d0cut = 0.5f; // 5mm extern bool cf_seeding; - // Config for propagation + // Config for propagation - could/should enter into PropagationFlags?! constexpr int Niter = 5; constexpr int NiterSim = 10; // Can make more steps due to near volume misses. constexpr bool useTrigApprox = true; - // Config for Bfield + // PropagationFlags as used during finding and fitting. Defined for each Geom in its plugin. + extern bool finding_requires_propagation_to_hit_pos; + extern PropagationFlags finding_inter_layer_pflags; + extern PropagationFlags finding_intra_layer_pflags; + extern PropagationFlags backward_fit_pflags; + extern PropagationFlags forward_fit_pflags; + extern PropagationFlags seed_fit_pflags; + + // Config for Bfield. Note: for now the same for CMS-2017 and CylCowWLids. constexpr float Bfield = 3.8112; constexpr float mag_c1 = 3.8114; constexpr float mag_b0 = -3.94991e-06; constexpr float mag_b1 = 7.53701e-06; constexpr float mag_a = 2.43878e-11; - + // Config for seeding as well... needed bfield constexpr float maxCurvR = (100 * minSimPt) / (sol * Bfield); // in cm diff --git a/Geoms/CMS-2017.cc b/Geoms/CMS-2017.cc index 1ff2683ec62d6..d676b884734b0 100644 --- a/Geoms/CMS-2017.cc +++ b/Geoms/CMS-2017.cc @@ -12,12 +12,20 @@ namespace void Create_CMS_2017(TrackerInfo& ti, bool verbose) { Config::nTotalLayers = 18 + 2 * 27; + Config::useCMSGeom = true; Config::nlayers_per_seed = 4; Config::maxCandsPerSeed = 6; // GC said 3 is enough ??? Config::maxHolesPerCand = 12; // should be reduced Config::chi2Cut = 30.0; + Config::finding_requires_propagation_to_hit_pos = true; + Config::finding_inter_layer_pflags = PropagationFlags(PF_use_param_b_field | PF_apply_material); + Config::finding_intra_layer_pflags = PropagationFlags(PF_none); + Config::backward_fit_pflags = PropagationFlags(PF_use_param_b_field | PF_apply_material); + Config::forward_fit_pflags = PropagationFlags(PF_use_param_b_field | PF_apply_material); + Config::seed_fit_pflags = PropagationFlags(PF_none); + ti.set_eta_regions(0.9, 1.7, 2.45, false); ti.create_layers(18, 27, 27); diff --git a/Geoms/CylCowWLids.cc b/Geoms/CylCowWLids.cc index 5df0500033e68..3209664458fb1 100644 --- a/Geoms/CylCowWLids.cc +++ b/Geoms/CylCowWLids.cc @@ -211,6 +211,13 @@ namespace { Config::nTotalLayers = 10 + 2 * 9; + Config::finding_requires_propagation_to_hit_pos = false; + Config::finding_inter_layer_pflags = PropagationFlags(PF_use_param_b_field); + Config::finding_intra_layer_pflags = PropagationFlags(PF_none); + Config::backward_fit_pflags = PropagationFlags(PF_use_param_b_field); + Config::forward_fit_pflags = PropagationFlags(PF_use_param_b_field); + Config::seed_fit_pflags = PropagationFlags(PF_none); + CylCowWLidsCreator creator(ti); creator.FillTrackerInfo(); diff --git a/mkFit/KalmanUtilsMPlex.cc b/mkFit/KalmanUtilsMPlex.cc index ced1c8b51e32f..97372f528ff28 100644 --- a/mkFit/KalmanUtilsMPlex.cc +++ b/mkFit/KalmanUtilsMPlex.cc @@ -470,14 +470,10 @@ void kalmanPropagateAndUpdate(const MPlexLS &psErr, const MPlexLV& psPar, const MPlexLS &outErr, MPlexLV& outPar, const int N_proc, const PropagationFlags propFlags) { - MPlexLS propErr; - MPlexLV propPar; - - // XXXX The following if could go away if we use proper update_foos in - // steering and set them from geom. For now it is still here (I'm going - // crazy already as things are). - - if (Config::useCMSGeom) { + if (Config::finding_requires_propagation_to_hit_pos) + { + MPlexLS propErr; + MPlexLV propPar; MPlexQF msRad; #pragma simd for (int n = 0; n < NN; ++n) @@ -486,13 +482,15 @@ void kalmanPropagateAndUpdate(const MPlexLS &psErr, const MPlexLV& psPar, const } propagateHelixToRMPlex(psErr, psPar, inChg, msRad, propErr, propPar, N_proc, propFlags); - } else { - propErr = psErr; - propPar = psPar; - } - kalmanOperation(KFO_Update_Params, propErr, propPar, msErr, msPar, - outErr, outPar, dummy_chi2, N_proc); + kalmanOperation(KFO_Update_Params, propErr, propPar, msErr, msPar, + outErr, outPar, dummy_chi2, N_proc); + } + else + { + kalmanOperation(KFO_Update_Params, psErr, psPar, msErr, msPar, + outErr, outPar, dummy_chi2, N_proc); + } } //------------------------------------------------------------------------------ @@ -511,10 +509,10 @@ void kalmanPropagateAndComputeChi2(const MPlexLS &psErr, const MPlexLV& psPar, MPlexQF& outChi2, const int N_proc, const PropagationFlags propFlags) { - MPlexLS propErr; - MPlexLV propPar; - - if (Config::useCMSGeom) { + if (Config::finding_requires_propagation_to_hit_pos) + { + MPlexLS propErr; + MPlexLV propPar; MPlexQF msRad; #pragma simd for (int n = 0; n < NN; ++n) @@ -523,13 +521,15 @@ void kalmanPropagateAndComputeChi2(const MPlexLS &psErr, const MPlexLV& psPar, } propagateHelixToRMPlex(psErr, psPar, inChg, msRad, propErr, propPar, N_proc, propFlags); - } else { - propErr = psErr; - propPar = psPar; - } - kalmanOperation(KFO_Calculate_Chi2, propErr, propPar, msErr, msPar, - dummy_err, dummy_par, outChi2, N_proc); + kalmanOperation(KFO_Calculate_Chi2, propErr, propPar, msErr, msPar, + dummy_err, dummy_par, outChi2, N_proc); + } + else + { + kalmanOperation(KFO_Calculate_Chi2, psErr, psPar, msErr, msPar, + dummy_err, dummy_par, outChi2, N_proc); + } } //------------------------------------------------------------------------------ @@ -688,14 +688,10 @@ void kalmanPropagateAndUpdateEndcap(const MPlexLS &psErr, const MPlexLV& psPar, MPlexLS &outErr, MPlexLV& outPar, const int N_proc, const PropagationFlags propFlags) { - MPlexLS propErr; - MPlexLV propPar; - - // XXXX The following if could go away if we use proper update_foos in - // steering and set them from geom. For now it is still here (I'm going - // crazy already as things are). - - if (Config::useCMSGeom) { + if (Config::finding_requires_propagation_to_hit_pos) + { + MPlexLS propErr; + MPlexLV propPar; MPlexQF msZ; #pragma simd for (int n = 0; n < NN; ++n) @@ -704,13 +700,15 @@ void kalmanPropagateAndUpdateEndcap(const MPlexLS &psErr, const MPlexLV& psPar, } propagateHelixToZMPlex(psErr, psPar, inChg, msZ, propErr, propPar, N_proc, propFlags); - } else { - propErr = psErr; - propPar = psPar; - } - kalmanOperationEndcap(KFO_Update_Params, propErr, propPar, msErr, msPar, - outErr, outPar, dummy_chi2, N_proc); + kalmanOperationEndcap(KFO_Update_Params, propErr, propPar, msErr, msPar, + outErr, outPar, dummy_chi2, N_proc); + } + else + { + kalmanOperationEndcap(KFO_Update_Params, psErr, psPar, msErr, msPar, + outErr, outPar, dummy_chi2, N_proc); + } } //------------------------------------------------------------------------------ @@ -729,10 +727,10 @@ void kalmanPropagateAndComputeChi2Endcap(const MPlexLS &psErr, const MPlexLV& p MPlexQF& outChi2, const int N_proc, const PropagationFlags propFlags) { - MPlexLS propErr; - MPlexLV propPar; - - if (Config::useCMSGeom) { + if (Config::finding_requires_propagation_to_hit_pos) + { + MPlexLS propErr; + MPlexLV propPar; MPlexQF msZ; #pragma simd for (int n = 0; n < NN; ++n) @@ -741,13 +739,15 @@ void kalmanPropagateAndComputeChi2Endcap(const MPlexLS &psErr, const MPlexLV& p } propagateHelixToZMPlex(psErr, psPar, inChg, msZ, propErr, propPar, N_proc, propFlags); - } else { - propErr = psErr; - propPar = psPar; - } - kalmanOperationEndcap(KFO_Calculate_Chi2, propErr, propPar, msErr, msPar, - dummy_err, dummy_par, outChi2, N_proc); + kalmanOperationEndcap(KFO_Calculate_Chi2, propErr, propPar, msErr, msPar, + dummy_err, dummy_par, outChi2, N_proc); + } + else + { + kalmanOperationEndcap(KFO_Calculate_Chi2, psErr, psPar, msErr, msPar, + dummy_err, dummy_par, outChi2, N_proc); + } } //------------------------------------------------------------------------------ diff --git a/mkFit/MkBuilder.cc b/mkFit/MkBuilder.cc index 445431a69e551..05a42d8effe54 100644 --- a/mkFit/MkBuilder.cc +++ b/mkFit/MkBuilder.cc @@ -735,7 +735,7 @@ inline void MkBuilder::fit_one_seed_set(TrackVec& seedtracks, int itrack, int en if (Config::seedInput != cmsswSeeds) { - mkfttr->FitTracksSteered(is_brl, end - itrack, m_event); + mkfttr->FitTracksSteered(is_brl, end - itrack, m_event, Config::seed_fit_pflags); } mkfttr->OutputFittedTracksAndHitIdx(m_event->seedTracks_, itrack, end, false); @@ -1244,8 +1244,7 @@ void MkBuilder::FindTracksBestHit() dcall(pre_prop_print(curr_layer, mkfndr.get())); (mkfndr.get()->*fnd_foos.m_propagate_foo)(layer_info.m_propagate_to, curr_tridx, - PF_use_param_b_field | PF_apply_material); - // XXXX-MFMAT-review + Config::finding_inter_layer_pflags); dcall(post_prop_print(curr_layer, mkfndr.get())); @@ -1410,7 +1409,7 @@ void MkBuilder::FindTracksStandard() dcall(pre_prop_print(curr_layer, mkfndr.get())); (mkfndr.get()->*fnd_foos.m_propagate_foo)(layer_info.m_propagate_to, end - itrack, - PF_use_param_b_field | PF_apply_material); + Config::finding_inter_layer_pflags); // XXXX-MFMAT-review @@ -1610,7 +1609,7 @@ void MkBuilder::find_tracks_in_layers(CandCloner &cloner, MkFinder *mkfndr, // propagate to current layer (mkfndr->*fnd_foos.m_propagate_foo)(layer_info.m_propagate_to, end - itrack, - PF_use_param_b_field | PF_apply_material); + Config::finding_inter_layer_pflags); // XXXX-MFMAT-review @@ -1812,7 +1811,7 @@ void MkBuilder::find_tracks_in_layersFV(int start_seed, int end_seed, int region auto& mkfndr = finders[index]; // propagate to current layer - (mkfndr.*fnd_foos.m_propagate_foo)(layer_info.m_propagate_to, mkfndr.nnfv(), PF_use_param_b_field); + (mkfndr.*fnd_foos.m_propagate_foo)(layer_info.m_propagate_to, mkfndr.nnfv(), Config::finding_inter_layer_pflags); mkfndr.SelectHitIndices(layer_of_hits); diff --git a/mkFit/MkFinder.cc b/mkFit/MkFinder.cc index 157d80a0cc84e..81f8cd92a609f 100644 --- a/mkFit/MkFinder.cc +++ b/mkFit/MkFinder.cc @@ -172,13 +172,16 @@ void MkFinder::SelectHitIndices(const LayerOfHits &layer_of_hits, float z = Par[iI].ConstAt(itrack, 2, 0); float dz = nSigmaZ * std::sqrt(Err[iI].ConstAt(itrack, 2, 2)); - dz = std::max(std::abs(dz), L.min_dq()); + dz = std::max(std::abs(dz), L.min_dq()); - if (Config::useCMSGeom) + // NOTE -- once issues in this block are resolved the changes should also be + // ported to MkFinderFV. + if (Config::useCMSGeom) // should be Config::finding_requires_propagation_to_hit_pos { //now correct for bending and for layer thickness unsing linear approximation //fixme! using constant value, to be taken from layer properties //XXXXMT4GC should we also increase dz? + //XXXXMT4GC an we just take half od layer dR? const float deltaR = Config::cmsDeltaRad; const float r = std::sqrt(r2); //here alpha is the difference between posPhi and momPhi @@ -207,11 +210,12 @@ void MkFinder::SelectHitIndices(const LayerOfHits &layer_of_hits, dr = std::max(std::abs(dr), L.min_dq()); - if (Config::useCMSGeom) + if (Config::useCMSGeom) // should be Config::finding_requires_propagation_to_hit_pos { //now correct for bending and for layer thickness unsing linear approximation //fixme! using constant value, to be taken from layer properties //XXXXMT4GC should we also increase dr? + //XXXXMT4GC can we just take half of layer dz? const float deltaZ = 5; float cosT = std::cos(Par[iI].ConstAt(itrack, 5, 0)); float sinT = std::sin(Par[iI].ConstAt(itrack, 5, 0)); @@ -434,7 +438,7 @@ void MkFinder::AddBestHit(const LayerOfHits &layer_of_hits, const int N_proc, //now compute the chi2 of track state vs hit MPlexQF outChi2; (*fnd_foos.m_compute_chi2_foo)(Err[iP], Par[iP], Chg, msErr, msPar, - outChi2, N_proc, PF_use_param_b_field); + outChi2, N_proc, Config::finding_intra_layer_pflags); #ifndef NO_PREFETCH // Prefetch to L1 the hits we'll process in the next loop iteration. @@ -536,7 +540,7 @@ void MkFinder::AddBestHit(const LayerOfHits &layer_of_hits, const int N_proc, dprint("update parameters"); (*fnd_foos.m_update_param_foo)(Err[iP], Par[iP], Chg, msErr, msPar, - Err[iC], Par[iC], N_proc, PF_use_param_b_field); + Err[iC], Par[iC], N_proc, Config::finding_intra_layer_pflags); //std::cout << "Par[iP](0,0,0)=" << Par[iP](0,0,0) << " Par[iC](0,0,0)=" << Par[iC](0,0,0)<< std::endl; } @@ -616,7 +620,7 @@ void MkFinder::FindCandidates(const LayerOfHits &layer_of_hits, //now compute the chi2 of track state vs hit MPlexQF outChi2; (*fnd_foos.m_compute_chi2_foo)(Err[iP], Par[iP], Chg, msErr, msPar, - outChi2, N_proc, PF_use_param_b_field); + outChi2, N_proc, Config::finding_intra_layer_pflags); // Prefetch to L1 the hits we'll (probably) process in the next loop iteration. for (int itrack = 0; itrack < N_proc; ++itrack) @@ -648,7 +652,7 @@ void MkFinder::FindCandidates(const LayerOfHits &layer_of_hits, if (oneCandPassCut) { (*fnd_foos.m_update_param_foo)(Err[iP], Par[iP], Chg, msErr, msPar, - Err[iC], Par[iC], N_proc, PF_use_param_b_field); + Err[iC], Par[iC], N_proc, Config::finding_intra_layer_pflags); dprint("update parameters" << std::endl << "propagated track parameters x=" << Par[iP].ConstAt(0, 0, 0) << " y=" << Par[iP].ConstAt(0, 1, 0) << std::endl @@ -788,7 +792,7 @@ void MkFinder::FindCandidatesCloneEngine(const LayerOfHits &layer_of_hits, CandC //now compute the chi2 of track state vs hit MPlexQF outChi2; - (*fnd_foos.m_compute_chi2_foo)(Err[iP], Par[iP], Chg, msErr, msPar, outChi2, N_proc, PF_use_param_b_field); + (*fnd_foos.m_compute_chi2_foo)(Err[iP], Par[iP], Chg, msErr, msPar, outChi2, N_proc, Config::finding_intra_layer_pflags); // Prefetch to L1 the hits we'll (probably) process in the next loop iteration. for (int itrack = 0; itrack < N_proc; ++itrack) @@ -880,7 +884,7 @@ void MkFinder::UpdateWithLastHit(const LayerOfHits &layer_of_hits, int N_proc, } (*fnd_foos.m_update_param_foo)(Err[iP], Par[iP], Chg, msErr, msPar, - Err[iC], Par[iC], N_proc, PF_use_param_b_field); + Err[iC], Par[iC], N_proc, Config::finding_intra_layer_pflags); //now that we have moved propagation at the end of the sequence we lost the handle of //using the propagated parameters instead of the updated for the missing hit case. @@ -1037,14 +1041,14 @@ void MkFinder::BkFitFitTracks(const EventOfHits & eventofhits, if (LI.is_barrel()) { - PropagateTracksToHitR(msPar, N_proc, PF_use_param_b_field | PF_apply_material); + PropagateTracksToHitR(msPar, N_proc, Config::backward_fit_pflags); kalmanOperation(KFO_Calculate_Chi2 | KFO_Update_Params, Err[iP], Par[iP], msErr, msPar, Err[iC], Par[iC], tmp_chi2, N_proc); } else { - PropagateTracksToHitZ(msPar, N_proc, PF_use_param_b_field | PF_apply_material); + PropagateTracksToHitZ(msPar, N_proc, Config::backward_fit_pflags); kalmanOperationEndcap(KFO_Calculate_Chi2 | KFO_Update_Params, Err[iP], Par[iP], msErr, msPar, Err[iC], Par[iC], tmp_chi2, N_proc); diff --git a/mkFit/MkFinderFV.cc b/mkFit/MkFinderFV.cc index 348e91895c059..83e8abe0b1204 100644 --- a/mkFit/MkFinderFV.cc +++ b/mkFit/MkFinderFV.cc @@ -318,7 +318,7 @@ void MkFinderFV<nseeds, ncands>::FindCandidates(const LayerOfHits &layer_of_hits #endif //now compute the chi2 of track state vs hit - (*fnd_foos.m_compute_chi2_foo)(Err[iP], Par[iP], Chg, msErr, msPar, XHitChi2[hit_cnt], NNFV, PF_use_param_b_field); + (*fnd_foos.m_compute_chi2_foo)(Err[iP], Par[iP], Chg, msErr, msPar, XHitChi2[hit_cnt], NNFV, Config::finding_intra_layer_pflags); // Prefetch to L1 the hits we'll (probably) process in the next loop iteration. for (int itrack = 0; itrack < NNFV; ++itrack) @@ -357,7 +357,7 @@ void MkFinderFV<nseeds, ncands>::UpdateWithLastHit(const LayerOfHits &layer_of_h } (*fnd_foos.m_update_param_foo)(Err[iP], Par[iP], Chg, msErr, msPar, - Err[iC], Par[iC], NNFV, PF_use_param_b_field); + Err[iC], Par[iC], NNFV, Config::finding_intra_layer_pflags); //now that we have moved propagation at the end of the sequence we lost the handle of //using the propagated parameters instead of the updated for the missing hit case. diff --git a/mkFit/MkFitter.cc b/mkFit/MkFitter.cc index ea036d65fd8a6..56bd7135e770c 100644 --- a/mkFit/MkFitter.cc +++ b/mkFit/MkFitter.cc @@ -444,7 +444,7 @@ void MkFitter::FitTracksWithInterSlurp(const std::vector<HitVec>& layersohits, msErr[0].SlurpIn(varr + off_error, idx, N_proc); #endif - PropagateTracksToHitR(msPar[0], N_proc, PF_use_param_b_field | PF_apply_material); + PropagateTracksToHitR(msPar[0], N_proc, Config::forward_fit_pflags); kalmanUpdate(Err[iP], Par[iP], msErr[0], msPar[0], Err[iC], Par[iC], N_proc); @@ -516,12 +516,10 @@ void MkFitter::ConformalFitTracks(bool fitting, int beg, int end) } } -void MkFitter::FitTracks(const int N_proc, const Event * ev, const bool useParamBfield) +void MkFitter::FitTracks(const int N_proc, const Event * ev, const PropagationFlags pflags) { // Fitting loop. - PropagationFlags pflags((useParamBfield ? PF_use_param_b_field : 0) | PF_apply_material); - for (int hi = 0; hi < Nhits; ++hi) { // Note, charge is not passed (line propagation). @@ -561,14 +559,12 @@ void MkFitter::CollectFitValidation(const int hi, const int N_proc, const Event } } -void MkFitter::FitTracksSteered(const bool is_barrel[], const int N_proc, const Event * ev, const bool useParamBfield) +void MkFitter::FitTracksSteered(const bool is_barrel[], const int N_proc, const Event * ev, const PropagationFlags pflags) { // Fitting loop. dprintf("MkFitter::FitTracksSteered %d %d %d\n", is_barrel[0], is_barrel[1], is_barrel[2]); - PropagationFlags pflags((useParamBfield ? PF_use_param_b_field : 0) | PF_apply_material); - for (int hi = 0; hi < Nhits; ++hi) { // Note, charge is not passed (line propagation). diff --git a/mkFit/MkFitter.h b/mkFit/MkFitter.h index dcc69a053d54f..86085aa0a92ef 100644 --- a/mkFit/MkFitter.h +++ b/mkFit/MkFitter.h @@ -75,8 +75,8 @@ class MkFitter : public MkBase void FitTracksWithInterSlurp(const std::vector<HitVec>& layersohits, int N_proc); void ConformalFitTracks(bool fitting, int beg, int end); - void FitTracks(const int N_proc, const Event * ev, const bool useParamBfield = false); - void FitTracksSteered(const bool is_barrel[], const int N_proc, const Event * ev, const bool useParamBfield = false); + void FitTracks(const int N_proc, const Event * ev, const PropagationFlags pflags); + void FitTracksSteered(const bool is_barrel[], const int N_proc, const Event * ev, const PropagationFlags pflags); void CollectFitValidation(const int hi, const int N_proc, const Event * ev) const; diff --git a/mkFit/PropagationMPlex.cc b/mkFit/PropagationMPlex.cc index c93529ed6ce52..0c09eca74e691 100644 --- a/mkFit/PropagationMPlex.cc +++ b/mkFit/PropagationMPlex.cc @@ -382,18 +382,18 @@ void helixAtRFromIterativeCCSFullJac(const MPlexLV& inPar, const MPlexQI& inChg, void helixAtRFromIterativeCCS(const MPlexLV& inPar, const MPlexQI& inChg, const MPlexQF &msRad, MPlexLV& outPar, MPlexLL& errorProp, - const int N_proc, const PropagationFlags pf) + const int N_proc, const PropagationFlags pflags) { errorProp.SetVal(0.f); - helixAtRFromIterativeCCS_impl(inPar, inChg, msRad, outPar, errorProp, 0, NN, N_proc, pf); + helixAtRFromIterativeCCS_impl(inPar, inChg, msRad, outPar, errorProp, 0, NN, N_proc, pflags); } void propagateHelixToRMPlex(const MPlexLS &inErr, const MPlexLV& inPar, const MPlexQI &inChg, const MPlexQF& msRad, MPlexLS &outErr, MPlexLV& outPar, - const int N_proc, const PropagationFlags pf) + const int N_proc, const PropagationFlags pflags) { outErr = inErr; outPar = inPar; // XXXX-4G is this requirement for helixAtRFromIterativeCCS_XXX ??? We should document it there. @@ -401,7 +401,7 @@ void propagateHelixToRMPlex(const MPlexLS &inErr, const MPlexLV& inPar, MPlexLL errorProp; - helixAtRFromIterativeCCS(inPar, inChg, msRad, outPar, errorProp, N_proc, pf.use_param_b_field); + helixAtRFromIterativeCCS(inPar, inChg, msRad, outPar, errorProp, N_proc, pflags.use_param_b_field); #ifdef DEBUG { @@ -421,7 +421,7 @@ void propagateHelixToRMPlex(const MPlexLS &inErr, const MPlexLV& inPar, } #endif - if (Config::useCMSGeom && pf.apply_material) + if (pflags.apply_material) { MPlexQF hitsRl; MPlexQF hitsXi; @@ -487,14 +487,14 @@ void propagateHelixToRMPlex(const MPlexLS &inErr, const MPlexLV& inPar, void propagateHelixToZMPlex(const MPlexLS &inErr, const MPlexLV& inPar, const MPlexQI &inChg, const MPlexQF& msZ, MPlexLS &outErr, MPlexLV& outPar, - const int N_proc, const PropagationFlags pf) + const int N_proc, const PropagationFlags pflags) { outErr = inErr; outPar = inPar; MPlexLL errorProp; - helixAtZ(inPar, inChg, msZ, outPar, errorProp, N_proc, pf); + helixAtZ(inPar, inChg, msZ, outPar, errorProp, N_proc, pflags); #ifdef DEBUG { @@ -514,7 +514,7 @@ void propagateHelixToZMPlex(const MPlexLS &inErr, const MPlexLV& inPar, } #endif - if (Config::useCMSGeom && pf.apply_material) + if (pflags.apply_material) { MPlexQF hitsRl; MPlexQF hitsXi; @@ -565,7 +565,7 @@ void propagateHelixToZMPlex(const MPlexLS &inErr, const MPlexLV& inPar, void helixAtZ(const MPlexLV& inPar, const MPlexQI& inChg, const MPlexQF &msZ, MPlexLV& outPar, MPlexLL& errorProp, - const int N_proc, const PropagationFlags pf) + const int N_proc, const PropagationFlags pflags) { errorProp.SetVal(0.f); @@ -586,7 +586,7 @@ void helixAtZ(const MPlexLV& inPar, const MPlexQI& inChg, const MPlexQF &msZ, const float phiin = inPar.ConstAt(n, 4, 0); const float theta = inPar.ConstAt(n, 5, 0); - const float k = inChg.ConstAt(n, 0, 0) * 100.f / (-Config::sol*(pf.use_param_b_field?Config::BfieldFromZR(zin,hipo(inPar.ConstAt(n,0,0),inPar.ConstAt(n,1,0))):Config::Bfield)); + const float k = inChg.ConstAt(n, 0, 0) * 100.f / (-Config::sol*(pflags.use_param_b_field?Config::BfieldFromZR(zin,hipo(inPar.ConstAt(n,0,0),inPar.ConstAt(n,1,0))):Config::Bfield)); dprint_np(n, std::endl << "input parameters" << " inPar.ConstAt(n, 0, 0)=" << std::setprecision(9) << inPar.ConstAt(n, 0, 0) diff --git a/mkFit/PropagationMPlex.h b/mkFit/PropagationMPlex.h index 394fd1a9f986c..b30fd3100e913 100644 --- a/mkFit/PropagationMPlex.h +++ b/mkFit/PropagationMPlex.h @@ -28,7 +28,7 @@ void propagateLineToRMPlex(const MPlexLS &psErr, const MPlexLV& psPar, void propagateHelixToRMPlex(const MPlexLS &inErr, const MPlexLV& inPar, const MPlexQI &inChg, const MPlexQF& msRad, MPlexLS &outErr, MPlexLV& outPar, - const int N_proc, const PropagationFlags pf); + const int N_proc, const PropagationFlags pflags); void helixAtRFromIterativeCCSFullJac(const MPlexLV& inPar, const MPlexQI& inChg, const MPlexQF &msRad, MPlexLV& outPar, MPlexLL& errorProp, @@ -36,16 +36,16 @@ void helixAtRFromIterativeCCSFullJac(const MPlexLV& inPar, const MPlexQI& inChg, void helixAtRFromIterativeCCS(const MPlexLV& inPar, const MPlexQI& inChg, const MPlexQF &msRad, MPlexLV& outPar, MPlexLL& errorProp, - const int N_proc, const PropagationFlags pf); + const int N_proc, const PropagationFlags pflags); void propagateHelixToZMPlex(const MPlexLS &inErr, const MPlexLV& inPar, const MPlexQI &inChg, const MPlexQF& msZ, MPlexLS &outErr, MPlexLV& outPar, - const int N_proc, const PropagationFlags pf); + const int N_proc, const PropagationFlags pflags); void helixAtZ(const MPlexLV& inPar, const MPlexQI& inChg, const MPlexQF &msZ, MPlexLV& outPar, MPlexLL& errorProp, - const int N_proc, const PropagationFlags pf); + const int N_proc, const PropagationFlags pflags); void applyMaterialEffects(const MPlexQF &hitsRl, const MPlexQF& hitsXi, MPlexLS &outErr, MPlexLV& outPar, From 46901dfcd3b6fb307b262e28dbfd23e0f4750a13 Mon Sep 17 00:00:00 2001 From: Matevz Tadel <mtadel@ucsd.edu> Date: Thu, 4 Jan 2018 15:41:46 -0800 Subject: [PATCH 06/11] Remove XXXX-MFMAT-review comments. --- mkFit/MkBuilder.cc | 2 -- 1 file changed, 2 deletions(-) diff --git a/mkFit/MkBuilder.cc b/mkFit/MkBuilder.cc index 05a42d8effe54..4670bc42dbe22 100644 --- a/mkFit/MkBuilder.cc +++ b/mkFit/MkBuilder.cc @@ -1410,7 +1410,6 @@ void MkBuilder::FindTracksStandard() (mkfndr.get()->*fnd_foos.m_propagate_foo)(layer_info.m_propagate_to, end - itrack, Config::finding_inter_layer_pflags); - // XXXX-MFMAT-review dcall(post_prop_print(curr_layer, mkfndr.get())); @@ -1610,7 +1609,6 @@ void MkBuilder::find_tracks_in_layers(CandCloner &cloner, MkFinder *mkfndr, // propagate to current layer (mkfndr->*fnd_foos.m_propagate_foo)(layer_info.m_propagate_to, end - itrack, Config::finding_inter_layer_pflags); - // XXXX-MFMAT-review // copy_out the propagated track params, errors only (hit-idcs and chi2 already updated) From bc7909c627276628f12981b066e4ab1a4c814c5b Mon Sep 17 00:00:00 2001 From: Matevz Tadel <mtadel@ucsd.edu> Date: Fri, 5 Jan 2018 11:05:31 -0800 Subject: [PATCH 07/11] Use non-squared disk hole limits, we know radius where we call the function now. --- TrackerInfo.cc | 4 ++-- TrackerInfo.h | 12 +++++------- mkFit/HitStructures.h | 2 -- 3 files changed, 7 insertions(+), 11 deletions(-) diff --git a/TrackerInfo.cc b/TrackerInfo.cc index ce7e1b6151ec8..92116755345ef 100644 --- a/TrackerInfo.cc +++ b/TrackerInfo.cc @@ -21,8 +21,8 @@ void LayerInfo::set_selection_limits(float p1, float p2, float q1, float q2) void LayerInfo::set_r_hole_range(float rh1, float rh2) { m_has_r_range_hole = true; - m_hole_r2_min = rh1 * rh1; - m_hole_r2_max = rh2 * rh2; + m_hole_r_min = rh1; + m_hole_r_max = rh2; } diff --git a/TrackerInfo.h b/TrackerInfo.h index 54305d211ae60..2f6cbe79d46de 100644 --- a/TrackerInfo.h +++ b/TrackerInfo.h @@ -33,7 +33,7 @@ struct WSR_Result class LayerInfo { private: - bool is_in_r2_hole(float r2) const { return r2 > m_hole_r2_min && r2 < m_hole_r2_max; } + bool is_in_r_hole_no_check(float r) const { return r > m_hole_r_min && r < m_hole_r_max; } public: enum LayerType_e { Undef = -1, Barrel = 0, EndCapPos = 1, EndCapNeg = 2 }; @@ -49,7 +49,7 @@ class LayerInfo bool m_is_outer = false; bool m_has_r_range_hole = false; - float m_hole_r2_min, m_hole_r2_max; // This could be turned into std::function when needed. + float m_hole_r_min, m_hole_r_max; // This could be turned into std::function when needed. // Selection limits float m_q_bin; // > 0 - bin width, < 0 - number of bins @@ -83,8 +83,7 @@ class LayerInfo bool is_within_r_limits(float r) const { return r > m_rin && r < m_rout; } bool is_within_q_limits(float q) const { return is_barrel() ? is_within_z_limits(q) : is_within_r_limits(q); } - bool is_in_rsqr_hole(float r2) const { return m_has_r_range_hole ? is_in_r2_hole(r2) : false; } - bool is_in_xy_hole(float x, float y) const { return m_has_r_range_hole ? is_in_r2_hole(x*x + y*y): false; } + bool is_in_r_hole (float r) const { return m_has_r_range_hole ? is_in_r_hole_no_check(r) : false; } WSR_Result is_within_z_sensitive_region(float z, float dz) const { @@ -100,9 +99,8 @@ class LayerInfo { if (m_has_r_range_hole) { - const float r2 = r*r; - if (r2 < m_hole_r2_max - dr && r2 > m_hole_r2_min + dr) return WSR_Result(WSR_Outside, true); - if (r2 < m_hole_r2_max + dr && r2 > m_hole_r2_min - dr ) return WSR_Result(WSR_Edge, true); + if (r < m_hole_r_max - dr && r > m_hole_r_min + dr) return WSR_Result(WSR_Outside, true); + if (r < m_hole_r_max + dr && r > m_hole_r_min - dr ) return WSR_Result(WSR_Edge, true); } return WSR_Result(WSR_Inside, false); } diff --git a/mkFit/HitStructures.h b/mkFit/HitStructures.h index 465691b6febb0..e4b4e786e23d1 100644 --- a/mkFit/HitStructures.h +++ b/mkFit/HitStructures.h @@ -107,8 +107,6 @@ class LayerOfHits float min_dq() const { return m_layer_info->m_select_min_dq; } float max_dq() const { return m_layer_info->m_select_max_dq; } - bool is_in_xy_hole(float x, float y) const { return m_layer_info->is_in_xy_hole(x, y); } - // Testing bin filling static constexpr float m_fphi = Config::m_nphi / Config::TwoPI; static constexpr int m_phi_mask = 0x3ff; From 4d917e3aebe1c151722e1002ef5796ecbbb3aee2 Mon Sep 17 00:00:00 2001 From: Matevz Tadel <mtadel@ucsd.edu> Date: Thu, 11 Jan 2018 14:15:32 -0800 Subject: [PATCH 08/11] Remove questions/comments for Giuseppe and Kevin. --- mkFit/PropagationMPlex.cc | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/mkFit/PropagationMPlex.cc b/mkFit/PropagationMPlex.cc index 0c09eca74e691..c6109c5406611 100644 --- a/mkFit/PropagationMPlex.cc +++ b/mkFit/PropagationMPlex.cc @@ -395,9 +395,12 @@ void propagateHelixToRMPlex(const MPlexLS &inErr, const MPlexLV& inPar, MPlexLS &outErr, MPlexLV& outPar, const int N_proc, const PropagationFlags pflags) { + // This is used further down when calculating similarity with errorProp (and before in DEBUG). + // MT: I don't think this really needed if we use inErr where required. outErr = inErr; - outPar = inPar; // XXXX-4G is this requirement for helixAtRFromIterativeCCS_XXX ??? We should document it there. - // Or even better, do it there as it is somewhat counterintuitive :) + // This requirement for helixAtRFromIterativeCCS_impl() and for helixAtRFromIterativeCCSFullJac(). + // MT: This should be properly handled in both functions (expecting input in output parameters sucks). + outPar = inPar; MPlexLL errorProp; @@ -428,14 +431,13 @@ void propagateHelixToRMPlex(const MPlexLS &inErr, const MPlexLV& inPar, #pragma simd for (int n = 0; n < NN; ++n) { - const int zbin = getZbinME(outPar(n, 2, 0)); // XXXX-4K changed msPar_z to outPar_z + const int zbin = getZbinME(outPar(n, 2, 0)); const int rbin = getRbinME(msRad (n, 0, 0)); hitsRl(n, 0, 0) = (zbin>=0 && zbin<Config::nBinsZME) ? getRlVal(zbin,rbin) : 0.f; // protect against crazy propagations hitsXi(n, 0, 0) = (zbin>=0 && zbin<Config::nBinsZME) ? getXiVal(zbin,rbin) : 0.f; // protect against crazy propagations } - applyMaterialEffects(hitsRl, hitsXi, outErr, outPar, N_proc); // XXXX-4K - why are we doing this on error before propagation ??? - // Isn't the material calculated for this/current layer? + applyMaterialEffects(hitsRl, hitsXi, outErr, outPar, N_proc); } squashPhiMPlex(outPar,N_proc); // ensure phi is between |pi| @@ -522,8 +524,8 @@ void propagateHelixToZMPlex(const MPlexLS &inErr, const MPlexLV& inPar, for (int n = 0; n < NN; ++n) { const int zbin = getZbinME(msZ(n, 0, 0)); - const int rbin = getRbinME(std::hypot(outPar(n, 0, 0), outPar(n, 1, 0))); // XXXX-4K changed msPar_xy to outPar_xy - + const int rbin = getRbinME(std::hypot(outPar(n, 0, 0), outPar(n, 1, 0))); + hitsRl(n, 0, 0) = (rbin>=0 && rbin<Config::nBinsRME) ? getRlVal(zbin,rbin) : 0.f; // protect against crazy propagations hitsXi(n, 0, 0) = (rbin>=0 && rbin<Config::nBinsRME) ? getXiVal(zbin,rbin) : 0.f; // protect against crazy propagations } From 865c46db3834e20c776ee6a81315c5998c58c20d Mon Sep 17 00:00:00 2001 From: Matevz Tadel <mtadel@ucsd.edu> Date: Thu, 11 Jan 2018 14:57:51 -0800 Subject: [PATCH 09/11] Backward-fit after finding for BestHit, emove obsolete BH classes in HitStructures.h. --- mkFit/HitStructures.h | 84 +---------------------- mkFit/MkBuilder.cc | 153 +++++++++++++++++++++++++++++------------- mkFit/MkBuilder.h | 1 + mkFit/MkFinder.cc | 92 ++++++++++++++++++++----- mkFit/MkFinder.h | 8 ++- 5 files changed, 191 insertions(+), 147 deletions(-) diff --git a/mkFit/HitStructures.h b/mkFit/HitStructures.h index e4b4e786e23d1..3aeef68fddef3 100644 --- a/mkFit/HitStructures.h +++ b/mkFit/HitStructures.h @@ -210,90 +210,11 @@ class EventOfHits }; + //============================================================================== +// CombinedCandidate and EventOfCombinedCandidates //============================================================================== -// This should actually be a BunchOfCandidates that share common hit vector. -// At the moment this is an EtaBin ... - -class EtaBinOfCandidates -{ -public: - std::vector<Track> m_candidates; - - int m_capacity; - int m_size; - -public: - EtaBinOfCandidates() : - m_candidates(Config::maxCandsPerEtaBin), - m_capacity (Config::maxCandsPerEtaBin), - m_size (0) - {} - - void Reset() - { - m_size = 0; - } - - void InsertTrack(const Track& track) - { - assert (m_size < m_capacity); // or something - - m_candidates[m_size] = track; - ++m_size; - } - - void SortByPhi() - { - std::sort(m_candidates.begin(), m_candidates.begin() + m_size, sortTrksByPhiMT); - } -}; - -class EventOfCandidates -{ -public: - std::vector<EtaBinOfCandidates> m_etabins_of_candidates; - -public: - EventOfCandidates() : - m_etabins_of_candidates(Config::nEtaBin) - {} - - void Reset() - { - for (auto &i : m_etabins_of_candidates) - { - i.Reset(); - } - } - - void InsertCandidate(const Track& track) - { - int bin = getEtaBinExtendedEdge(track.posEta()); - // XXXX MT Had to add back this conditional for best-hit (bad seeds) - // Note also the ExtendedEdge above, this practically removes bin = -1 - // occurence and improves efficiency. - if (bin != -1) - { - m_etabins_of_candidates[bin].InsertTrack(track); - } - } - - void SortByPhi() - { - for (auto &i : m_etabins_of_candidates) - { - i.SortByPhi(); - } - } -}; - - -//------------------------------------------------------- -// for combinatorial version, switch to vector of vectors -//------------------------------------------------------- - // This inheritance sucks but not doing it will require more changes. class CombCandidate : public std::vector<Track> @@ -305,6 +226,7 @@ class CombCandidate : public std::vector<Track> int m_last_seed_layer = -1; }; + class EventOfCombCandidates { public: diff --git a/mkFit/MkBuilder.cc b/mkFit/MkBuilder.cc index 4670bc42dbe22..e484e2201c488 100644 --- a/mkFit/MkBuilder.cc +++ b/mkFit/MkBuilder.cc @@ -1263,6 +1263,12 @@ void MkBuilder::FindTracksBestHit() mkfndr->OutputTracksAndHitIdx(cands, trk_idcs, 0, curr_tridx, false); + // final backward fit + if (Config::backwardFit) + { + BackwardFitBH(mkfndr.get(), rng.m_beg, rng.m_end, region); + } + ++rng; } // end of loop over candidates in a tbb chunk }); // end parallel_for over candidates in a region @@ -1671,54 +1677,6 @@ void MkBuilder::find_tracks_in_layers(CandCloner &cloner, MkFinder *mkfndr, } } -//============================================================================== - -void MkBuilder::BackwardFit(MkFinder *mkfndr, int start_seed, int end_seed, int region) -{ - EventOfCombCandidates &eoccs = m_event_of_comb_cands; - const SteeringParams &st_par = m_steering_params[region]; - - for (int iseed = start_seed; iseed < end_seed; iseed += NN) - { - const int end = std::min(iseed + NN, end_seed); - - // printf("Pre Final fit for %d - %d\n", iseed, end); - // for (int i = iseed; i < end; ++i) { const Track &t = eoccs[i][0]; - // printf(" %4d with q=%+d chi2=%7.3f pT=%7.3f eta=% 7.3f x=%.3f y=%.3f z=%.3f nHits=%2d label=%4d findable=%d\n", - // i, t.charge(), t.chi2(), t.pT(), t.momEta(), t.x(), t.y(), t.z(), t.nFoundHits(), t.label(), t.isFindable()); - // } - - bool chi_debug = false; - redo_fit: - - mkfndr->BkFitInputTracks(eoccs, iseed, end); - mkfndr->BkFitFitTracks(m_event_of_hits, st_par, end - iseed, false, chi_debug); - -#ifdef DEBUG_BACKWARD_FIT - // Dump tracks with pT > 2 and chi2/dof > 20. Assumes MPT_SIZE=1. - if (! chi_debug && 1.0f/mkfndr->Par[MkBase::iC].At(0,3,0) > 2.0f && - mkfndr->Chi2(0,0,0) / (eoccs[iseed][0].nFoundHits() * 3 - 6) > 20.0f) - { - chi_debug = true; - printf("CHIHDR Event %d, Seed %3d, pT %f, chipdof %f ### NOTE x,y,z in cm, sigmas, deltas in mum ### !!!\n", - m_event->evtID(), iseed, 1.0f/mkfndr->Par[MkBase::iC].At(0,3,0), - mkfndr->Chi2(0,0,0) / (eoccs[iseed][0].nFoundHits() * 3 - 6)); - printf("CHIHDR %3s %10s %10s %10s %10s %10s %11s %11s %11s %10s %10s %10s %10s %11s %11s %11s %10s %10s %10s %10s %10s %11s %11s\n", - "lyr","chi2","x_h","y_h","z_h","r_h","sx_h","sy_h","sz_h","x_t","y_t","z_t","r_t","sx_t","sy_t","sz_t","pt","phi","theta","phi_h","phi_t","d_xy","d_z"); - goto redo_fit; - } -#endif - - mkfndr->BkFitOutputTracks(eoccs, iseed, end); - - // printf("Post Final fit for %d - %d\n", iseed, end); - // for (int i = iseed; i < end; ++i) { const Track &t = eoccs[i][0]; - // printf(" %4d with q=%+d chi2=%7.3f pT=%7.3f eta=% 7.3f x=%.3f y=%.3f z=%.3f nHits=%2d label=%4d findable=%d\n", - // i, t.charge(), t.chi2(), t.pT(), t.momEta(), t.x(), t.y(), t.z(), t.nFoundHits(), t.label(), t.isFindable()); - // } - } - -} //------------------------------------------------------------------------------ // FindTracksCombinatorial: FullVector TBB @@ -1844,3 +1802,102 @@ void MkBuilder::find_tracks_in_layersFV(int start_seed, int end_seed, int region } #endif } + + +//============================================================================== +// BackwardFit +//============================================================================== + +void MkBuilder::BackwardFitBH(MkFinder *mkfndr, int start_seed, int end_seed, int region) +{ + TrackVec &cands = m_event->candidateTracks_; + + const SteeringParams &st_par = m_steering_params[region]; + + for (int iseed = start_seed; iseed < end_seed; iseed += NN) + { + const int end = std::min(iseed + NN, end_seed); + + // printf("Pre Final fit for %d - %d\n", iseed, end); + // for (int i = iseed; i < end; ++i) { const Track &t = eoccs[i][0]; + // printf(" %4d with q=%+d chi2=%7.3f pT=%7.3f eta=% 7.3f x=%.3f y=%.3f z=%.3f nHits=%2d label=%4d findable=%d\n", + // i, t.charge(), t.chi2(), t.pT(), t.momEta(), t.x(), t.y(), t.z(), t.nFoundHits(), t.label(), t.isFindable()); + // } + + bool chi_debug = false; + redo_fit: + + mkfndr->BkFitInputTracks(cands, iseed, end); + mkfndr->BkFitFitTracks(m_event_of_hits, st_par, end - iseed, false, chi_debug); + +#ifdef DEBUG_BACKWARD_FIT + // Dump tracks with pT > 2 and chi2/dof > 20. Assumes MPT_SIZE=1. + if (! chi_debug && 1.0f/mkfndr->Par[MkBase::iC].At(0,3,0) > 2.0f && + mkfndr->Chi2(0,0,0) / (eoccs[iseed][0].nFoundHits() * 3 - 6) > 20.0f) + { + chi_debug = true; + printf("CHIHDR Event %d, Seed %3d, pT %f, chipdof %f ### NOTE x,y,z in cm, sigmas, deltas in mum ### !!!\n", + m_event->evtID(), iseed, 1.0f/mkfndr->Par[MkBase::iC].At(0,3,0), + mkfndr->Chi2(0,0,0) / (eoccs[iseed][0].nFoundHits() * 3 - 6)); + printf("CHIHDR %3s %10s %10s %10s %10s %10s %11s %11s %11s %10s %10s %10s %10s %11s %11s %11s %10s %10s %10s %10s %10s %11s %11s\n", + "lyr","chi2","x_h","y_h","z_h","r_h","sx_h","sy_h","sz_h","x_t","y_t","z_t","r_t","sx_t","sy_t","sz_t","pt","phi","theta","phi_h","phi_t","d_xy","d_z"); + goto redo_fit; + } +#endif + + mkfndr->BkFitOutputTracks(cands, iseed, end); + + // printf("Post Final fit for %d - %d\n", iseed, end); + // for (int i = iseed; i < end; ++i) { const Track &t = eoccs[i][0]; + // printf(" %4d with q=%+d chi2=%7.3f pT=%7.3f eta=% 7.3f x=%.3f y=%.3f z=%.3f nHits=%2d label=%4d findable=%d\n", + // i, t.charge(), t.chi2(), t.pT(), t.momEta(), t.x(), t.y(), t.z(), t.nFoundHits(), t.label(), t.isFindable()); + // } + } + +} + +void MkBuilder::BackwardFit(MkFinder *mkfndr, int start_seed, int end_seed, int region) +{ + EventOfCombCandidates &eoccs = m_event_of_comb_cands; + const SteeringParams &st_par = m_steering_params[region]; + + for (int iseed = start_seed; iseed < end_seed; iseed += NN) + { + const int end = std::min(iseed + NN, end_seed); + + // printf("Pre Final fit for %d - %d\n", iseed, end); + // for (int i = iseed; i < end; ++i) { const Track &t = eoccs[i][0]; + // printf(" %4d with q=%+d chi2=%7.3f pT=%7.3f eta=% 7.3f x=%.3f y=%.3f z=%.3f nHits=%2d label=%4d findable=%d\n", + // i, t.charge(), t.chi2(), t.pT(), t.momEta(), t.x(), t.y(), t.z(), t.nFoundHits(), t.label(), t.isFindable()); + // } + + bool chi_debug = false; + redo_fit: + + mkfndr->BkFitInputTracks(eoccs, iseed, end); + mkfndr->BkFitFitTracks(m_event_of_hits, st_par, end - iseed, false, chi_debug); + +#ifdef DEBUG_BACKWARD_FIT + // Dump tracks with pT > 2 and chi2/dof > 20. Assumes MPT_SIZE=1. + if (! chi_debug && 1.0f/mkfndr->Par[MkBase::iC].At(0,3,0) > 2.0f && + mkfndr->Chi2(0,0,0) / (eoccs[iseed][0].nFoundHits() * 3 - 6) > 20.0f) + { + chi_debug = true; + printf("CHIHDR Event %d, Seed %3d, pT %f, chipdof %f ### NOTE x,y,z in cm, sigmas, deltas in mum ### !!!\n", + m_event->evtID(), iseed, 1.0f/mkfndr->Par[MkBase::iC].At(0,3,0), + mkfndr->Chi2(0,0,0) / (eoccs[iseed][0].nFoundHits() * 3 - 6)); + printf("CHIHDR %3s %10s %10s %10s %10s %10s %11s %11s %11s %10s %10s %10s %10s %11s %11s %11s %10s %10s %10s %10s %10s %11s %11s\n", + "lyr","chi2","x_h","y_h","z_h","r_h","sx_h","sy_h","sz_h","x_t","y_t","z_t","r_t","sx_t","sy_t","sz_t","pt","phi","theta","phi_h","phi_t","d_xy","d_z"); + goto redo_fit; + } +#endif + + mkfndr->BkFitOutputTracks(eoccs, iseed, end); + + // printf("Post Final fit for %d - %d\n", iseed, end); + // for (int i = iseed; i < end; ++i) { const Track &t = eoccs[i][0]; + // printf(" %4d with q=%+d chi2=%7.3f pT=%7.3f eta=% 7.3f x=%.3f y=%.3f z=%.3f nHits=%2d label=%4d findable=%d\n", + // i, t.charge(), t.chi2(), t.pT(), t.momEta(), t.x(), t.y(), t.z(), t.nFoundHits(), t.label(), t.isFindable()); + // } + } +} diff --git a/mkFit/MkBuilder.h b/mkFit/MkBuilder.h index 47542e224066f..62fff55ac163d 100644 --- a/mkFit/MkBuilder.h +++ b/mkFit/MkBuilder.h @@ -168,6 +168,7 @@ class MkBuilder void FindTracksCloneEngine(); void FindTracksFV(); + void BackwardFitBH(MkFinder *mkfndr, int start_seed, int end_seed, int region); void BackwardFit(MkFinder *mkfndr, int start_seed, int end_seed, int region); #ifdef USE_CUDA diff --git a/mkFit/MkFinder.cc b/mkFit/MkFinder.cc index 81f8cd92a609f..609e4737a02b0 100644 --- a/mkFit/MkFinder.cc +++ b/mkFit/MkFinder.cc @@ -949,6 +949,64 @@ void MkFinder::CopyOutParErr(std::vector<CombCandidate>& seed_cand_vec, // Backward Fit hack //============================================================================== +void MkFinder::BkFitInputTracks(TrackVec& cands, int beg, int end) +{ + // SlurpIn based on XHit array - so Nhits is irrelevant. + // Could as well use HotArrays from tracks directly + a local cursor array to last hit. + + const int N_proc = end - beg; + const Track &trk0 = cands[beg]; + const char *varr = (char*) &trk0; + const int off_error = (char*) trk0.errors().Array() - varr; + const int off_param = (char*) trk0.parameters().Array() - varr; + + int idx[NN] __attribute__((aligned(64))); + + int itrack = 0; + + for (int i = beg; i < end; ++i, ++itrack) + { + const Track &trk = cands[i]; + + Chg(itrack, 0, 0) = trk.charge(); + CurHit[itrack] = trk.nTotalHits() - 1; + HoTArr[itrack] = trk.getHitsOnTrackArray(); + + idx[itrack] = (char*) &trk - varr; + } + + Chi2.SetVal(0); + +#ifdef MIC_INTRINSICS + __m512i vi = _mm512_load_epi32(idx); + Err[iC].SlurpIn(varr + off_error, vi, N_proc); + Par[iC].SlurpIn(varr + off_param, vi, N_proc); +#else + Err[iC].SlurpIn(varr + off_error, idx, N_proc); + Par[iC].SlurpIn(varr + off_param, idx, N_proc); +#endif + + Err[iC].Scale(100.0f); +} + +void MkFinder::BkFitOutputTracks(TrackVec& cands, int beg, int end) +{ + // Only copy out track params / errors / chi2, all the rest is ok. + + int itrack = 0; + for (int i = beg; i < end; ++i, ++itrack) + { + Track &trk = cands[i]; + + Err[iC].CopyOut(itrack, trk.errors_nc().Array()); + Par[iC].CopyOut(itrack, trk.parameters_nc().Array()); + + trk.setChi2(Chi2(itrack, 0, 0)); + } +} + +//------------------------------------------------------------------------------ + void MkFinder::BkFitInputTracks(EventOfCombCandidates& eocss, int beg, int end) { // SlurpIn based on XHit array - so Nhits is irrelevant. @@ -989,6 +1047,24 @@ void MkFinder::BkFitInputTracks(EventOfCombCandidates& eocss, int beg, int end) Err[iC].Scale(100.0f); } +void MkFinder::BkFitOutputTracks(EventOfCombCandidates& eocss, int beg, int end) +{ + // Only copy out track params / errors / chi2, all the rest is ok. + + int itrack = 0; + for (int i = beg; i < end; ++i, ++itrack) + { + Track &trk = eocss[i][0]; + + Err[iC].CopyOut(itrack, trk.errors_nc().Array()); + Par[iC].CopyOut(itrack, trk.parameters_nc().Array()); + + trk.setChi2(Chi2(itrack, 0, 0)); + } +} + +//------------------------------------------------------------------------------ + namespace { float e2s(float x) { return 1e4 * std::sqrt(x); } } void MkFinder::BkFitFitTracks(const EventOfHits & eventofhits, @@ -1089,19 +1165,3 @@ void MkFinder::BkFitFitTracks(const EventOfHits & eventofhits, Chi2.Add(tmp_chi2); } } - -void MkFinder::BkFitOutputTracks(EventOfCombCandidates& eocss, int beg, int end) -{ - // Only copy out track params / errors / chi2, all the rest is ok. - - int itrack = 0; - for (int i = beg; i < end; ++i, ++itrack) - { - Track &trk = eocss[i][0]; - - Err[iC].CopyOut(itrack, trk.errors_nc().Array()); - Par[iC].CopyOut(itrack, trk.parameters_nc().Array()); - - trk.setChi2(Chi2(itrack, 0, 0)); - } -} diff --git a/mkFit/MkFinder.h b/mkFit/MkFinder.h index 36b1dbde14443..ae4f4bd1ae3a1 100644 --- a/mkFit/MkFinder.h +++ b/mkFit/MkFinder.h @@ -148,10 +148,14 @@ class MkFinder : public MkBase int CurHit[NN]; const HitOnTrack *HoTArr[NN]; - void BkFitInputTracks(EventOfCombCandidates& eocss, int beg, int end); + void BkFitInputTracks (TrackVec& cands, int beg, int end); + void BkFitOutputTracks(TrackVec& cands, int beg, int end); + + void BkFitInputTracks (EventOfCombCandidates& eocss, int beg, int end); + void BkFitOutputTracks(EventOfCombCandidates& eocss, int beg, int end); + void BkFitFitTracks(const EventOfHits& eventofhits, const SteeringParams& st_par, int N_proc, bool useParamBfield = false, bool chiDebug = false); - void BkFitOutputTracks(EventOfCombCandidates& eocss, int beg, int end); //---------------------------------------------------------------------------- From 371f171fb28b859e14be7deb590f6b3570da5dd9 Mon Sep 17 00:00:00 2001 From: Kevin McDermott <kpm82@cornell.edu> Date: Mon, 15 Jan 2018 10:58:28 -0500 Subject: [PATCH 10/11] remove param b-field from inter-layer propagation during finding --- Geoms/CMS-2017.cc | 2 +- Geoms/CylCowWLids.cc | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Geoms/CMS-2017.cc b/Geoms/CMS-2017.cc index d676b884734b0..7555c0134df92 100644 --- a/Geoms/CMS-2017.cc +++ b/Geoms/CMS-2017.cc @@ -20,7 +20,7 @@ namespace Config::chi2Cut = 30.0; Config::finding_requires_propagation_to_hit_pos = true; - Config::finding_inter_layer_pflags = PropagationFlags(PF_use_param_b_field | PF_apply_material); + Config::finding_inter_layer_pflags = PropagationFlags(PF_apply_material); Config::finding_intra_layer_pflags = PropagationFlags(PF_none); Config::backward_fit_pflags = PropagationFlags(PF_use_param_b_field | PF_apply_material); Config::forward_fit_pflags = PropagationFlags(PF_use_param_b_field | PF_apply_material); diff --git a/Geoms/CylCowWLids.cc b/Geoms/CylCowWLids.cc index 3209664458fb1..665622f90e175 100644 --- a/Geoms/CylCowWLids.cc +++ b/Geoms/CylCowWLids.cc @@ -212,7 +212,7 @@ namespace Config::nTotalLayers = 10 + 2 * 9; Config::finding_requires_propagation_to_hit_pos = false; - Config::finding_inter_layer_pflags = PropagationFlags(PF_use_param_b_field); + Config::finding_inter_layer_pflags = PropagationFlags(PF_none); Config::finding_intra_layer_pflags = PropagationFlags(PF_none); Config::backward_fit_pflags = PropagationFlags(PF_use_param_b_field); Config::forward_fit_pflags = PropagationFlags(PF_use_param_b_field); From 51ef65fcc9a6b5b4b5da08b6083a14cfeb20af7f Mon Sep 17 00:00:00 2001 From: Kevin McDermott <kpm82@cornell.edu> Date: Mon, 15 Jan 2018 12:43:04 -0800 Subject: [PATCH 11/11] protect against insane propagations which return bins out of range of binned material constants --- mkFit/PropagationMPlex.cc | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/mkFit/PropagationMPlex.cc b/mkFit/PropagationMPlex.cc index c6109c5406611..1bd16b2771ba2 100644 --- a/mkFit/PropagationMPlex.cc +++ b/mkFit/PropagationMPlex.cc @@ -434,8 +434,8 @@ void propagateHelixToRMPlex(const MPlexLS &inErr, const MPlexLV& inPar, const int zbin = getZbinME(outPar(n, 2, 0)); const int rbin = getRbinME(msRad (n, 0, 0)); - hitsRl(n, 0, 0) = (zbin>=0 && zbin<Config::nBinsZME) ? getRlVal(zbin,rbin) : 0.f; // protect against crazy propagations - hitsXi(n, 0, 0) = (zbin>=0 && zbin<Config::nBinsZME) ? getXiVal(zbin,rbin) : 0.f; // protect against crazy propagations + hitsRl(n, 0, 0) = (zbin>=0 && zbin<Config::nBinsZME && rbin>=0 && rbin<Config::nBinsRME) ? getRlVal(zbin,rbin) : 0.f; // protect against crazy propagations + hitsXi(n, 0, 0) = (zbin>=0 && zbin<Config::nBinsZME && rbin>=0 && rbin<Config::nBinsRME) ? getXiVal(zbin,rbin) : 0.f; // protect against crazy propagations } applyMaterialEffects(hitsRl, hitsXi, outErr, outPar, N_proc); } @@ -526,8 +526,8 @@ void propagateHelixToZMPlex(const MPlexLS &inErr, const MPlexLV& inPar, const int zbin = getZbinME(msZ(n, 0, 0)); const int rbin = getRbinME(std::hypot(outPar(n, 0, 0), outPar(n, 1, 0))); - hitsRl(n, 0, 0) = (rbin>=0 && rbin<Config::nBinsRME) ? getRlVal(zbin,rbin) : 0.f; // protect against crazy propagations - hitsXi(n, 0, 0) = (rbin>=0 && rbin<Config::nBinsRME) ? getXiVal(zbin,rbin) : 0.f; // protect against crazy propagations + hitsRl(n, 0, 0) = (zbin>=0 && zbin<Config::nBinsZME && rbin>=0 && rbin<Config::nBinsRME) ? getRlVal(zbin,rbin) : 0.f; // protect against crazy propagations + hitsXi(n, 0, 0) = (zbin>=0 && zbin<Config::nBinsZME && rbin>=0 && rbin<Config::nBinsRME) ? getXiVal(zbin,rbin) : 0.f; // protect against crazy propagations } applyMaterialEffects(hitsRl, hitsXi, outErr, outPar, N_proc); }