Skip to content

Commit

Permalink
Merge pull request cms-sw#334 from slava77/from-devel-pr332/ptSign-cosah
Browse files Browse the repository at this point in the history 
invpt sign fix and numerical precision improvements for high pt
  • Loading branch information
@cerati
cerati committed Jul 23, 2021
2 parents 2c817f5 + 77a7f37 commit 795127e
Show file tree
Hide file tree
Showing 6 changed files with 82 additions and 43 deletions.
2 changes: 1 addition & 1 deletion mkFit/CandCloner.cc
View file
Expand Up @@ -62,7 +62,7 @@ void CandCloner::ProcessSeedRange(int is_beg, int is_end)
dprint("trkIdx=" << hitsForSeed[ih].trkIdx << " hitIdx=" << hitsForSeed[ih].hitIdx << " chi2=" << hitsForSeed[ih].chi2 << std::endl
<< " "
<< "original pt=" << ccand[hitsForSeed[ih].trkIdx].pT() << " "
<< "nTotalHits=" << cdand[hitsForSeed[ih].trkIdx].nTotalHits() << " "
<< "nTotalHits=" << ccand[hitsForSeed[ih].trkIdx].nTotalHits() << " "
<< "nFoundHits=" << ccand[hitsForSeed[ih].trkIdx].nFoundHits() << " "
<< "chi2=" << ccand[hitsForSeed[ih].trkIdx].chi2());
}
Expand Down
11 changes: 11 additions & 0 deletions mkFit/MkFinder.cc
View file
Expand Up @@ -1309,6 +1309,17 @@ void MkFinder::BkFitFitTracksBH(const EventOfHits & eventofhits,
Err[iP], Par[iP], msErr, msPar, Err[iC], Par[iC], tmp_chi2, N_proc);
}

//fixup invpt sign and charge
for (int n = 0; n < N_proc; ++n)
{
if (Par[iC].At(n,3,0) < 0)
{
Chg.At(n, 0, 0) = -Chg.At(n, 0, 0);
Par[iC].At(n,3,0) = -Par[iC].At(n,3,0);
}
}


#ifdef DEBUG_BACKWARD_FIT_BH
// Dump per hit chi2
for (int i = 0; i < N_proc; ++i)
Expand Down
12 changes: 10 additions & 2 deletions mkFit/MkStdSeqs.cc
View file
Expand Up @@ -631,6 +631,9 @@ void find_duplicates_sharedhits_pixelseed(TrackVec &tracks, const float fraction

void find_and_remove_duplicates(TrackVec &tracks, const IterationConfig &itconf)
{
#ifdef DEBUG
std::cout<<" find_and_remove_duplicates: input track size " <<tracks.size()<<std::endl;
#endif
if (itconf.m_require_quality_filter)
{
find_duplicates_sharedhits(tracks, itconf.m_params.fracSharedHits);
Expand All @@ -645,6 +648,12 @@ void find_and_remove_duplicates(TrackVec &tracks, const IterationConfig &itconf)
remove_duplicates(tracks);
}

#ifdef DEBUG
std::cout<<" find_and_remove_duplicates: output track size " <<tracks.size()<<std::endl;
for (auto const& tk : tracks) {
std::cout<<tk.parameters()<<std::endl;
}
#endif
}

//=========================================================================
Expand Down Expand Up @@ -675,8 +684,7 @@ void dump_simtracks(Event *m_event)
{
dprint("track #" << itrack << " hit #" << ihit
<< " hit pos=" << simtracks[itrack].hitsVector(m_event->layerHits_)[ihit].position()
<< " phi=" << simtracks[itrack].hitsVector(m_event->layerHits_)[ihit].phi()
<< " phiPart=" << getPhiPartition(simtracks[itrack].hitsVector(m_event->layerHits_)[ihit].phi()));
<< " phi=" << simtracks[itrack].hitsVector(m_event->layerHits_)[ihit].phi());
}
}
}
Expand Down
74 changes: 43 additions & 31 deletions mkFit/PropagationMPlex.cc
View file
Expand Up @@ -295,7 +295,7 @@ void helixAtRFromIterativeCCSFullJac(const MPlexLV& inPar, const MPlexQI& inChg,
errorPropTmp(n,4,4) = 1.f;
errorPropTmp(n,5,5) = 1.f;

float cosa = 0., sina = 0.;
float cosah = 0., sinah = 0.;
//no trig approx here, phi and theta can be large
float cosP = std::cos(phiin), sinP = std::sin(phiin);
const float cosT = std::cos(theta), sinT = std::sin(theta);
Expand All @@ -313,19 +313,21 @@ void helixAtRFromIterativeCCSFullJac(const MPlexLV& inPar, const MPlexQI& inChg,
//alpha+=ialpha;

if (Config::useTrigApprox) {
sincos4(ialpha, sina, cosa);
sincos4(ialpha*0.5f, sinah, cosah);
} else {
cosa=std::cos(ialpha);
sina=std::sin(ialpha);
cosah=std::cos(ialpha*0.5f);
sinah=std::sin(ialpha*0.5f);
}
const float cosa=1.f-2.f*sinah*sinah;
const float sina=2.f*sinah*cosah;

//derivatives of alpha
const float dadx = -outPar.At(n, 0, 0)*ipt/(k*r0);
const float dady = -outPar.At(n, 1, 0)*ipt/(k*r0);
const float dadipt = (r-r0)/k;

outPar.At(n, 0, 0) = outPar.ConstAt(n, 0, 0) + k*(pxin*sina-pyin*(1.f-cosa));
outPar.At(n, 1, 0) = outPar.ConstAt(n, 1, 0) + k*(pyin*sina+pxin*(1.f-cosa));
outPar.At(n, 0, 0) = outPar.ConstAt(n, 0, 0) + 2.f*k*sinah*(pxin*cosah-pyin*sinah);
outPar.At(n, 1, 0) = outPar.ConstAt(n, 1, 0) + 2.f*k*sinah*(pyin*cosah+pxin*sinah);
const float pxinold = pxin;//copy before overwriting
pxin = pxin*cosa-pyin*sina;
pyin = pyin*cosa+pxinold*sina;
Expand Down Expand Up @@ -615,6 +617,7 @@ void helixAtZ(const MPlexLV& inPar, const MPlexQI& inChg, const MPlexQF &msZ,
const float theta = inPar.ConstAt(n, 5, 0);

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));
const float kinv = 1.f/k;

dprint_np(n, std::endl << "input parameters"
<< " inPar.ConstAt(n, 0, 0)=" << std::setprecision(9) << inPar.ConstAt(n, 0, 0)
Expand All @@ -627,10 +630,11 @@ void helixAtZ(const MPlexLV& inPar, const MPlexQI& inChg, const MPlexQF &msZ,

const float pt = 1.f/ipt;

float cosaTmp = 0., sinaTmp = 0.;
float cosahTmp = 0., sinahTmp = 0.;
//no trig approx here, phi can be large
const float cosP = std::cos(phiin), sinP = std::sin(phiin);
const float cosT = std::cos(theta), sinT = std::sin(theta);
const float tanT = sinT/cosT;
const float pxin = cosP*pt;
const float pyin = sinP*pt;

Expand All @@ -640,42 +644,50 @@ void helixAtZ(const MPlexLV& inPar, const MPlexQI& inChg, const MPlexQF &msZ,
<< " sinP=" << std::setprecision(9) << sinP << " pt=" << std::setprecision(9) << pt);

const float deltaZ = zout - zin;
const float alpha = deltaZ*sinT*ipt/(cosT*k);
const float alpha = deltaZ*tanT*ipt*kinv;

if (Config::useTrigApprox) {
sincos4(alpha, sinaTmp, cosaTmp);
sincos4(alpha*0.5f, sinahTmp, cosahTmp);
} else {
cosaTmp=std::cos(alpha);
sinaTmp=std::sin(alpha);
cosahTmp=std::cos(alpha*0.5f);
sinahTmp=std::sin(alpha*0.5f);
}
const float cosa = cosaTmp;
const float sina = sinaTmp;
const float cosah = cosahTmp;
const float sinah = sinahTmp;
const float cosa = 1.f - 2.f*sinah*sinah;
const float sina = 2.f*sinah*cosah;

//update parameters
outPar.At(n, 0, 0) = outPar.At(n, 0, 0) + k*(pxin*sina - pyin*(1.f-cosa));
outPar.At(n, 1, 0) = outPar.At(n, 1, 0) + k*(pyin*sina + pxin*(1.f-cosa));
outPar.At(n, 0, 0) = outPar.At(n, 0, 0) + 2.f*k*sinah*(pxin*cosah - pyin*sinah);
outPar.At(n, 1, 0) = outPar.At(n, 1, 0) + 2.f*k*sinah*(pyin*cosah + pxin*sinah);
outPar.At(n, 2, 0) = zout;
outPar.At(n, 4, 0) = phiin+alpha;

dprint_np(n, std::endl << "outPar.At(n, 0, 0)=" << outPar.At(n, 0, 0) << " outPar.At(n, 1, 0)=" << outPar.At(n, 1, 0)
<< " pxin=" << pxin << " pyin=" << pyin);

const float sCosPsina = std::sin(cosP*sina);
const float cCosPsina = std::cos(cosP*sina);

errorProp(n,0,2) = cosP*sinT*(sinP*cosa*sCosPsina - cosa)/cosT;
errorProp(n,0,3) = cosP*sinT*deltaZ*cosa*( 1.f - sinP*sCosPsina )/(cosT*ipt) - k*(cosP*sina - sinP*(1.f-cCosPsina))/(ipt*ipt);
errorProp(n,0,4) = (k/ipt)*( -sinP*sina + sinP*sinP*sina*sCosPsina - cosP*(1.f - cCosPsina ) );
errorProp(n,0,5) = cosP*deltaZ*cosa*( 1.f - sinP*sCosPsina )/(cosT*cosT);

errorProp(n,1,2) = cosa*sinT*(cosP*cosP*sCosPsina - sinP)/cosT;
errorProp(n,1,3) = sinT*deltaZ*cosa*( cosP*cosP*sCosPsina + sinP )/(cosT*ipt) - k*(sinP*sina + cosP*(1.f-cCosPsina))/(ipt*ipt);
errorProp(n,1,4) = (k/ipt)*( -sinP*(1.f - cCosPsina) - sinP*cosP*sina*sCosPsina + cosP*sina );
const float sCosPsinah = std::sin(cosP*sina*0.5f);
const float cCosPsinah = std::cos(cosP*sina*0.5f);
const float sCosPsina = 2.f*sCosPsinah*cCosPsinah;
const float sPsCPless1 = sinP*sCosPsina - 1.f;
const float cAlessFsAA = alpha != 0 ? cosa - sina/alpha : 0.f; // lim-> alpha^2/3
const float sCPx2 = alpha != 0 ? 2.f*sCosPsinah*(cosP*cosa*cCosPsinah - sCosPsinah/alpha) : 0.f;// lim-> ~alpha

errorProp(n,0,2) = cosP*tanT*cosa*sPsCPless1;
// errorProp(n,0,3) = cosP*tanT*deltaZ*cosa*( 1.f - sinP*sCosPsina )*pt - k*(cosP*sina - sinP*(1.f-cCosPsina))*pt*pt;
errorProp(n,0,3) = tanT*deltaZ*pt*( cosP*cAlessFsAA - sinP*sCPx2 );
errorProp(n,0,4) = (k*pt)*( sinP*sina*sPsCPless1 - 2.f*cosP*sCosPsinah*sCosPsinah );
errorProp(n,0,5) = -cosP*deltaZ*cosa*sPsCPless1/(cosT*cosT);

errorProp(n,1,2) = cosa*tanT*(cosP*cosP*sCosPsina - sinP);
// errorProp(n,1,3) = tanT*deltaZ*cosa*( cosP*cosP*sCosPsina + sinP )*pt - k*(sinP*sina + cosP*(1.f-cCosPsina))*pt*pt;
errorProp(n,1,3) = tanT*deltaZ*pt*( sinP*cAlessFsAA + cosP*sCPx2 );
errorProp(n,1,4) = (k*pt)*( -cosP*sina*sPsCPless1 - 2.f*sinP*sCosPsinah*sCosPsinah );
errorProp(n,1,5) = deltaZ*cosa*( cosP*cosP*sCosPsina + sinP )/(cosT*cosT);

errorProp(n,4,2) = -ipt*sinT/(cosT*k);
errorProp(n,4,3) = sinT*deltaZ/(cosT*k);
errorProp(n,4,5) = ipt*deltaZ/(cosT*cosT*k);
errorProp(n,4,2) = -ipt*tanT*kinv;
errorProp(n,4,3) = tanT*deltaZ*kinv;
errorProp(n,4,5) = ipt*deltaZ*kinv/(cosT*cosT);

dprint_np(n, "propagation end, dump parameters" << std::endl
<< "pos = " << outPar.At(n, 0, 0) << " " << outPar.At(n, 1, 0) << " " << outPar.At(n, 2, 0) << std::endl
Expand Down Expand Up @@ -717,7 +729,7 @@ void applyMaterialEffects(const MPlexQF &hitsRl, const MPlexQF& hitsXi, const MP
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 = (isBarrel ? p/pt : 1./std::abs(std::cos(theta)) );
const float invCos = (isBarrel ? p/pt : 1.f/std::abs(std::cos(theta)) );
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
Expand Down Expand Up @@ -746,7 +758,7 @@ void applyMaterialEffects(const MPlexQF &hitsRl, const MPlexQF& hitsXi, const MP
// 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 = propSign.ConstAt(n,0,0)*dEdx/beta;
outPar.At(n, 3, 0) = p/((p+dP)*pt);
outPar.At(n, 3, 0) = p/(std::max(p+dP,0.001f)*pt);//stay above 1MeV
//assume 100% uncertainty
outErr.At(n, 3, 3) += dP*dP/(p2*pt*pt);
}
Expand Down
14 changes: 8 additions & 6 deletions mkFit/PropagationMPlex.icc
View file
Expand Up @@ -54,7 +54,7 @@ static inline void helixAtRFromIterativeCCS_impl(const Tf& __restrict__ inPar
const float kinv = 1.f / k;
const float pt = 1.f / ipt;

float D = 0., cosa = 0., sina = 0., id = 0.;
float D = 0., cosa = 0., sina = 0., cosah = 0., sinah = 0., id = 0.;
//no trig approx here, phi can be large
float cosPorT = std::cos(phiin), sinPorT = std::sin(phiin);
float pxin = cosPorT*pt;
Expand Down Expand Up @@ -97,11 +97,13 @@ static inline void helixAtRFromIterativeCCS_impl(const Tf& __restrict__ inPar
D += id;

if (Config::useTrigApprox) {
sincos4(id*ipt*kinv, sina, cosa);
sincos4(id*ipt*kinv*0.5f, sinah, cosah);
} else {
cosa=std::cos(id*ipt*kinv);
sina=std::sin(id*ipt*kinv);
cosah=std::cos(id*ipt*kinv*0.5f);
sinah=std::sin(id*ipt*kinv*0.5f);
}
cosa=1.f-2.f*sinah*sinah;
sina=2.f*sinah*cosah;

dprint_np(n, "Attempt propagation from r=" << r0 << " to r=" << r << std::endl
<< " x=" << xin << " y=" << yin << " z=" << inPar(n, 2, 0)
Expand Down Expand Up @@ -142,8 +144,8 @@ static inline void helixAtRFromIterativeCCS_impl(const Tf& __restrict__ inPar
}

//update parameters
outPar(n, 0, 0) = outPar(n, 0, 0) + k*(pxin*sina - pyin*(1.f-cosa));
outPar(n, 1, 0) = outPar(n, 1, 0) + k*(pyin*sina + pxin*(1.f-cosa));
outPar(n, 0, 0) = outPar(n, 0, 0) + 2.f*k*sinah*(pxin*cosah - pyin*sinah);
outPar(n, 1, 0) = outPar(n, 1, 0) + 2.f*k*sinah*(pyin*cosah + pxin*sinah);
const float pxinold = pxin;//copy before overwriting
pxin = pxin*cosa - pyin*sina;
pyin = pyin*cosa + pxinold*sina;
Expand Down
12 changes: 9 additions & 3 deletions mkFit/seedtestMPlex.cc
View file
Expand Up @@ -34,7 +34,9 @@ void findSeedsByRoadSearch(TripletIdxConVec & seed_idcs, std::vector<LayerOfHits

// MIMI hack: Config::nlayers_per_seed = 4
// const float seed_z2cut = (Config::nlayers_per_seed * Config::fRadialSpacing) / std::tan(2.0f*std::atan(std::exp(-1.0f*Config::dEtaSeedTrip)));
// const float seed_z2cut = (4 * Config::fRadialSpacing) / std::tan(2.0f*std::atan(std::exp(-1.0f*Config::dEtaSeedTrip)));
#ifdef DEBUG
const float seed_z2cut = (4 * Config::fRadialSpacing) / std::tan(2.0f*std::atan(std::exp(-1.0f*Config::dEtaSeedTrip)));
#endif

// 0 = first layer, 1 = second layer, 2 = third layer
const LayerOfHits& lay1_hits = evt_lay_hits[1];
Expand Down Expand Up @@ -72,7 +74,9 @@ void findSeedsByRoadSearch(TripletIdxConVec & seed_idcs, std::vector<LayerOfHits
// negative points of intersection with third layer
float lay2_negx = 0.0f, lay2_negy = 0.0f;
intersectThirdLayer(aneg,bneg,hit1_x,hit1_y,lay2_negx,lay2_negy);
// MIMI const float lay2_negphi = getPhi(lay2_negx,lay2_negy);
#ifdef DEBUG
const float lay2_negphi = getPhi(lay2_negx,lay2_negy);
#endif

// center of positive curved track
const float apos = 0.5f*((hit0_x+hit1_x)+(hit0_y-hit1_y)*quad);
Expand All @@ -81,7 +85,9 @@ void findSeedsByRoadSearch(TripletIdxConVec & seed_idcs, std::vector<LayerOfHits
// positive points of intersection with third layer
float lay2_posx = 0.0f, lay2_posy = 0.0f;
intersectThirdLayer(apos,bpos,hit1_x,hit1_y,lay2_posx,lay2_posy);
// MIMI const float lay2_posphi = getPhi(lay2_posx,lay2_posy);
#ifdef DEBUG
const float lay2_posphi = getPhi(lay2_posx,lay2_posy);
#endif

std::vector<int> cand_hit2_indices;
// MIMI lay2_hits.SelectHitIndices((2.0f*hit1_z-hit0_z),(lay2_posphi+lay2_negphi)/2.0f,
Expand Down

0 comments on commit 795127e

Please sign in to comment.