192 LogTrace(
"TestOutliers") <<
"tOld size=" << tracksOld->size() <<
" tOut size=" << tracksOut->size()
193 <<
" aOld size=" << recSimCollOld.
size() <<
" aOut size=" << recSimCollOut.
size();
196 LogTrace(
"TestOutliers") <<
"recSimCollOld.size()=" << recSimCollOld.
size() ;
200 LogTrace(
"TestOutliers") <<
"trackOld->pt()=" << trackOld->pt() <<
" trackOld->numberOfValidHits()=" << trackOld->numberOfValidHits();
201 std::vector<pair<TrackingParticleRef, double> > tpOld;
202 if(recSimCollOld.
find(trackOld) != recSimCollOld.
end()){
203 tpOld = recSimCollOld[trackOld];
204 if (tpOld.size()!=0)
LogTrace(
"TestOutliers") <<
" associated" ;
205 else LogTrace(
"TestOutliers") <<
" NOT associated" ;
206 }
else LogTrace(
"TestOutliers") <<
" NOT associated" ;
208 LogTrace(
"TestOutliers") <<
"recSimCollOut.size()=" << recSimCollOut.
size() ;
212 LogTrace(
"TestOutliers") <<
"trackOut->pt()=" << trackOut->pt() <<
" trackOut->numberOfValidHits()=" << trackOut->numberOfValidHits();
213 std::vector<pair<TrackingParticleRef, double> > tpOut;
214 if(recSimCollOut.
find(trackOut) != recSimCollOut.
end()){
215 tpOut = recSimCollOut[trackOut];
216 if (tpOut.size()!=0)
LogTrace(
"TestOutliers") <<
" associated" ;
217 else LogTrace(
"TestOutliers") <<
" NOT associated" ;
218 }
else LogTrace(
"TestOutliers") <<
" NOT associated" ;
222 LogTrace(
"TestOutliers") <<
"tracksOut->size()=" << tracksOut->size();
223 LogTrace(
"TestOutliers") <<
"tracksOld->size()=" << tracksOld->size();
225 std::vector<unsigned int> outused;
226 for (
unsigned int j = 0;
j < tracksOld->size(); ++
j) {
229 LogTrace(
"TestOutliers") <<
"now track old with id=" <<
j <<
" seed ref=" << trackOld->seedRef().get()
230 <<
" pt=" << trackOld->pt() <<
" eta=" << trackOld->eta()
231 <<
" chi2=" << trackOld->normalizedChi2()
232 <<
" tip= " << fabs(trackOld->dxy(
beamSpot->position()))
233 <<
" lip= " << fabs(trackOld->dsz(
beamSpot->position()));
248 std::vector<unsigned int> outtest;
249 for (
unsigned int k = 0;
k < tracksOut->size(); ++
k) {
251 if (tmpOut->
seedRef() == trackOld->seedRef()) {
252 outtest.push_back(
k);
257 if (outtest.size() == 1) {
259 LogTrace(
"TestOutliers") <<
"now track out with id=" << outtest[0] <<
" seed ref=" << trackOut->
seedRef().
get()
260 <<
" pt=" << trackOut->
pt();
261 outused.push_back(outtest[0]);
262 }
else if (outtest.size() > 1) {
263 for (
unsigned int p = 0;
p < outtest.size(); ++
p) {
267 if ((*itOut)->isValid()) {
268 bool thishit =
false;
270 if ((*itOld)->isValid()) {
284 LogTrace(
"TestOutliers") <<
"now track out with id=" << outtest[
p]
285 <<
" seed ref=" << trackOut->
seedRef().
get() <<
" pt=" << trackOut->
pt();
286 outused.push_back(outtest[
p]);
291 if (outtest.empty() || trackOut.
get() ==
nullptr) {
292 if (recSimCollOld.
find(trackOld) != recSimCollOld.
end()) {
293 vector<pair<TrackingParticleRef, double> > tpOld;
294 tpOld = recSimCollOld[trackOld];
295 if (!tpOld.empty()) {
296 LogTrace(
"TestOutliers") <<
"no match: old associated and out lost! old has #hits="
297 << trackOld->numberOfValidHits() <<
" and fraction=" << tpOld.begin()->second;
298 if (tpOld.begin()->second > 0.5)
302 LogTrace(
"TestOutliers") <<
"...skip to next old track";
308 <<
" trackOld->seedRef()=" << trackOld->seedRef().get();
309 bool oldAssoc = recSimCollOld.
find(trackOld) != recSimCollOld.
end();
310 bool outAssoc = recSimCollOut.
find(trackOut) != recSimCollOut.
end();
311 LogTrace(
"TestOutliers") <<
"outAssoc=" << outAssoc <<
" oldAssoc=" << oldAssoc;
339 std::vector<unsigned int> tpids;
340 std::vector<std::pair<TrackingParticleRef, double> > tpOut;
341 std::vector<pair<TrackingParticleRef, double> > tpOld;
345 tpOut = recSimCollOut[trackOut];
346 if (!tpOut.empty()) {
348 tprOut = tpOut.begin()->first;
349 fracOut = tpOut.begin()->second;
351 tpids.push_back(g4T->trackId());
357 tpOld = recSimCollOld[trackOld];
358 if (!tpOld.empty()) {
359 tprOld = tpOld.begin()->first;
365 tpids.push_back(g4T->trackId());
371 if (tprOut.
get() !=
nullptr || tprOld.
get() !=
nullptr) {
373 tpr = tprOut.
get() !=
nullptr ? tprOut : tprOld;
375 const SimTrack *assocTrack = &(*tpr->g4Track_begin());
382 (*(trackOld->recHitsEnd() - 1)),
388 <<
"Out->pt=" << trackOut->
pt() <<
" Old->pt=" << trackOld->pt() <<
" tp->pt="
389 <<
sqrt(tpr->momentum().perp2())
391 <<
" Old->validHits=" << trackOld->numberOfValidHits() <<
" Out->validHits="
394 <<
" fracOut=" << fracOut <<
" deltaHits=" << trackOld->numberOfValidHits() - trackOut->
numberOfValidHits();
397 double PtPullOut = (trackOut->
pt() -
sqrt(tpr->momentum().perp2())) / trackOut->
ptError();
398 double PtPullOld = (trackOld->pt() -
sqrt(tpr->momentum().perp2())) / trackOld->ptError();
406 GlobalPoint(tpr->vertex().x(), tpr->vertex().y(), tpr->vertex().z()),
412 tscpBuilder(ftsAtProduction,
GlobalPoint(0, 0, 0));
417 double qoverpSim = tsAtClosestApproach.
charge() /
p.mag();
418 double lambdaSim =
M_PI / 2 -
p.theta();
419 double phiSim =
p.phi();
420 double dxySim = (-
v.x() *
sin(
p.phi()) +
v.y() *
cos(
p.phi()));
421 double dszSim =
v.z() *
p.perp() /
p.mag() - (
v.x() *
p.x() +
v.y() *
p.y()) /
p.perp() *
p.z() /
p.mag();
422 double d0Sim = -dxySim;
423 double dzSim = dszSim *
p.mag() /
p.perp();
426 double qoverpPullOut = (trackOut->
qoverp() - qoverpSim) / trackOut->
qoverpError();
427 double qoverpPullOld = (trackOld->qoverp() - qoverpSim) / trackOld->qoverpError();
428 double lambdaPullOut = (trackOut->
lambda() - lambdaSim) / trackOut->
thetaError();
429 double lambdaPullOld = (trackOld->lambda() - lambdaSim) / trackOld->thetaError();
430 double phi0PullOut = (trackOut->
phi() - phiSim) / trackOut->
phiError();
431 double phi0PullOld = (trackOld->phi() - phiSim) / trackOld->phiError();
432 double d0PullOut = (trackOut->
d0() - d0Sim) / trackOut->
d0Error();
433 double d0PullOld = (trackOld->d0() - d0Sim) / trackOld->d0Error();
434 double dzPullOut = (trackOut->
dz() - dzSim) / trackOut->
dzError();
435 double dzPullOld = (trackOld->dz() - dzSim) / trackOld->dzError();
453 if (tprOut.
get() !=
nullptr && tprOld.
get() ==
nullptr) {
454 if (!tpOld.empty() && tpOld.begin()->second <= 0.5) {
458 <<
" old #hits=" << trackOld->numberOfValidHits();
463 <<
" old #hits=" << trackOld->numberOfValidHits();
465 }
else if (tprOut.
get() ==
nullptr && tprOld.
get() !=
nullptr) {
466 LogTrace(
"TestOutliers") <<
"old associated and out not! old has #hits=" << trackOld->numberOfValidHits()
467 <<
" and fraction=" << tpOld.begin()->second;
468 if (tpOld.begin()->second > 0.5) {
474 if (fabs(PtPullOut) < fabs(PtPullOld))
483 LogTrace(
"TestOutliers") <<
"track old";
485 LogTrace(
"TestOutliers") << (*itOld)->isValid() <<
" " << (*itOld)->geographicalId().rawId();
487 LogTrace(
"TestOutliers") <<
"track out";
489 LogTrace(
"TestOutliers") << (*itOut)->isValid() <<
" " << (*itOut)->geographicalId().rawId();
493 vector<pair<int, trackingRecHit_iterator> > gainedlostoutliers;
497 if ((*itOut)->isValid()) {
499 if ((*itOld)->geographicalId().rawId() == (*itOut)->geographicalId().rawId())
503 gainedlostoutliers.push_back(pair<int, trackingRecHit_iterator>(1, itOut));
504 LogTrace(
"TestOutliers") <<
"broken trajectory during old fit... gained hit "
505 << (*itOut)->geographicalId().rawId();
513 bool outlier =
false;
517 if ((*itOld)->geographicalId().rawId() == (*itOut)->geographicalId().rawId()) {
519 if ((*itOld)->isValid() && !(*itOut)->isValid() &&
520 (*itOld)->geographicalId().rawId() == (*itOut)->geographicalId().rawId()) {
521 LogTrace(
"TestOutliers") << (*itOld)->isValid() <<
" " << (*itOut)->isValid() <<
" "
522 << (*itOld)->geographicalId().rawId() <<
" "
523 << (*itOut)->geographicalId().rawId();
529 gainedlostoutliers.push_back(pair<int, trackingRecHit_iterator>(2, itOld));
533 gainedlostoutliers.push_back(pair<int, trackingRecHit_iterator>(3, itOld));
536 for (
std::vector<pair<int, trackingRecHit_iterator> >::iterator it = gainedlostoutliers.begin();
537 it != gainedlostoutliers.end();
539 LogTrace(
"TestOutliers") <<
"type of processed hit:" << it->first;
543 bool outlier =
false;
546 else if (it->first == 2)
548 else if (it->first == 3)
551 if (outlier || lost || gained) {
554 if (lost && (*itHit)->isValid() ==
false) {
556 LogTrace(
"TestOutliers") <<
"lost invalid";
558 }
else if (gained && (*itHit)->isValid() ==
false) {
560 LogTrace(
"TestOutliers") <<
"gained invalid";
566 std::vector<SimHitIdpr> simTrackIds =
hitAssociator.associateHitId(**itHit);
567 bool goodhit =
false;
568 for (
size_t j = 0;
j < simTrackIds.size();
j++) {
569 for (
size_t jj = 0;
jj < tpids.size();
jj++) {
570 if (simTrackIds[
j].
first == tpids[
jj])
580 if (dynamic_cast<const SiPixelRecHit *>(&**itHit)) {
581 LogTrace(
"TestOutliers") <<
"SiPixelRecHit";
582 clustersize = ((
const SiPixelRecHit *)(&**itHit))->cluster()->size();
584 }
else if (dynamic_cast<const SiStripRecHit2D *>(&**itHit)) {
585 LogTrace(
"TestOutliers") <<
"SiStripRecHit2D";
586 clustersize = ((
const SiStripRecHit2D *)(&**itHit))->cluster()->amplitudes().size();
588 }
else if (dynamic_cast<const SiStripMatchedRecHit2D *>(&**itHit)) {
589 LogTrace(
"TestOutliers") <<
"SiStripMatchedRecHit2D";
592 if (clsize1 > clsize2)
593 clustersize = clsize1;
595 clustersize = clsize2;
597 }
else if (dynamic_cast<const ProjectedSiStripRecHit2D *>(&**itHit)) {
598 LogTrace(
"TestOutliers") <<
"ProjectedSiStripRecHit2D";
605 int subdetId = (*itHit)->geographicalId().subdetId();
606 DetId id = (*itHit)->geographicalId();
607 int layerId = tTopo->
layer(
id);
608 layerval = subdetId * 10 + layerId;
626 double energyLoss_ = 0.;
627 unsigned int monoId = 0;
628 std::vector<double> energyLossM;
629 std::vector<double> energyLossS;
630 std::vector<PSimHit> assSimHits =
hitAssociator.associateHit(**itHit);
631 if (assSimHits.empty())
634 std::vector<unsigned int> trackIds;
638 for (std::vector<PSimHit>::const_iterator
m = assSimHits.begin();
m < assSimHits.end();
m++) {
641 unsigned int tId =
m->trackId();
642 if (
find(trackIds.begin(), trackIds.end(), tId) == trackIds.end())
643 trackIds.push_back(tId);
644 LogTrace(
"TestOutliers") <<
"id=" << tId;
645 if (
m->energyLoss() > energyLoss_) {
647 energyLoss_ =
m->energyLoss();
649 if (hittypeval == 3) {
651 monoId =
m->detUnitId();
652 if (monoId ==
m->detUnitId()) {
653 energyLossM.push_back(
m->energyLoss());
655 energyLossS.push_back(
m->energyLoss());
659 energyLossM.push_back(
m->energyLoss());
662 unsigned int nIds = trackIds.size();
666 posxy->Fill(fabs(gpos.
x()), fabs(gpos.
y()));
667 poszr->Fill(fabs(gpos.
z()),
sqrt(gpos.
x() * gpos.
x() + gpos.
y() * gpos.
y()));
674 bool ioniOnly =
true;
675 unsigned int idc = 0;
676 for (
size_t jj = 0;
jj < tpids.size();
jj++) {
677 idc +=
std::count(trackIds.begin(), trackIds.end(), tpids[
jj]);
679 if (idc == trackIds.size()) {
682 for (std::vector<PSimHit>::const_iterator
m = assSimHits.begin() + 1;
m < assSimHits.end();
m++) {
683 if ((
m->processType() != 7 &&
m->processType() != 8 &&
m->processType() != 9) &&
684 abs(
m->particleType()) != 11) {
689 if (ioniOnly && !shared) {
690 LogTrace(
"TestOutliers") <<
"delta";
696 if (clustersize >= 4) {
704 if (hittypeval == 1 && clustersize >= 4)
706 if (hittypeval == 1 && clustersize < 4)
708 if (hittypeval == 2 && clustersize >= 4)
710 if (hittypeval == 2 && clustersize < 4)
712 if (hittypeval == 3 && clustersize >= 4)
714 if (hittypeval == 3 && clustersize < 4)
716 if (hittypeval == 4 && clustersize >= 4)
718 if (hittypeval == 4 && clustersize < 4)
741 if (hittypeval != 3) {
742 if (energyLossM.size() > 1) {
743 sort(energyLossM.begin(), energyLossM.end(), greater<double>());
748 if (energyLossM.size() > 1 && energyLossS.size() <= 1) {
749 sort(energyLossM.begin(), energyLossM.end(), greater<double>());
751 }
else if (energyLossS.size() > 1 && energyLossM.size() <= 1) {
752 sort(energyLossS.begin(), energyLossS.end(), greater<double>());
754 }
else if (energyLossS.size() > 1 && energyLossM.size() > 1) {
755 sort(energyLossM.begin(), energyLossM.end(), greater<double>());
756 sort(energyLossS.begin(), energyLossS.end(), greater<double>());
757 energyLossM[1] / energyLossM[0] > energyLossS[1] / energyLossS[0]
763 LogTrace(
"TestOutliers") <<
"before merged";
766 LogTrace(
"TestOutliers") <<
"assSimHits.size()=" << assSimHits.size();
767 if ((assSimHits.size() > 1 &&
tmp ==
nullptr) || (assSimHits.size() > 2 &&
tmp !=
nullptr)) {
774 for (std::vector<PSimHit>::const_iterator
m = assSimHits.begin();
m < assSimHits.end();
m++) {
775 unsigned int tId =
m->trackId();
776 LogTrace(
"TestOutliers") <<
"component with id=" << tId <<
" eLoss=" <<
m->energyLoss()
777 <<
" proc=" <<
m->processType() <<
" part=" <<
m->particleType();
778 if (
find(tpids.begin(), tpids.end(), tId) == tpids.end())
780 if (
m->processType() == 2) {
784 .
transform((*itHit)->localPositionError(),
785 theG->
idToDet((*itHit)->geographicalId())->surface())
791 ROOT::Math::SVector<double, 3>
delta;
798 LogTrace(
"TestOutliers") <<
"hit pull=" << mpull;
802 LogTrace(
"TestOutliers") <<
"end merged";
811 const SiPixelRecHit *pix = dynamic_cast<const SiPixelRecHit *>(&**itHit);
812 if ((hittypeval != 3 && assSimHits.size() < 2) || (hittypeval == 3 && assSimHits.size() < 3)) {
821 LogTrace(
"TestOutliers") <<
"out good";
824 LogTrace(
"TestOutliers") <<
"lost good";
827 LogTrace(
"TestOutliers") <<
"gained good";
834 if (ioniOnly && !shared) {
840 }
else if (!ioniOnly && !shared) {
853 LogTrace(
"TestOutliers") <<
"out merged, ioniOnly=" << ioniOnly <<
" shared=" << shared;
855 if (ioniOnly && !shared)
857 else if (!ioniOnly && !shared)
861 LogTrace(
"TestOutliers") <<
"lost merged, ioniOnly=" << ioniOnly <<
" shared=" << shared;
863 if (ioniOnly && !shared)
865 else if (!ioniOnly && !shared)
869 LogTrace(
"TestOutliers") <<
"gained merged, ioniOnly=" << ioniOnly <<
" shared=" << shared;
871 LogTrace(
"TestOutliers") <<
"merged, ioniOnly=" << ioniOnly <<
" shared=" << shared;
882 const SiPixelRecHit *pix = dynamic_cast<const SiPixelRecHit *>(&**itHit);
887 LogTrace(
"TestOutliers") <<
"out bad";
890 LogTrace(
"TestOutliers") <<
"lost bad";
893 LogTrace(
"TestOutliers") <<
"gained bad";
905 <<
"Out->pt=" << trackOut->
pt() <<
" Old->pt=" << trackOld->pt() <<
" tp->pt="
906 <<
sqrt(tpr->momentum().perp2())
908 <<
" Old->validHits=" << trackOld->numberOfValidHits() <<
" Out->validHits="
911 <<
" fracOut=" << fracOut <<
" deltaHits=" << trackOld->numberOfValidHits() - trackOut->
numberOfValidHits();
912 LogTrace(
"TestOutliers") <<
"track with gained hits";
918 LogTrace(
"TestOutliers") <<
"end track old #" <<
j;
921 for (
unsigned int k = 0;
k < tracksOut->size(); ++
k) {
922 if (
find(outused.begin(), outused.end(),
k) == outused.end()) {
924 bool outAssoc = recSimCollOut.
find(trackOut) != recSimCollOut.
end();