7 #include <TDirectory.h>
19 LogTrace(
"TestTrackHits") << iConfig;
45 file =
new TFile(
out.c_str(),
"recreate");
46 for (
int i = 0;
i != 6;
i++)
47 for (
int j = 0;
j != 9;
j++) {
62 title <<
"Chi2Increment_" <<
i + 1 <<
"-" <<
j + 1;
65 title <<
"Chi2IncrementVsEta_" <<
i + 1 <<
"-" <<
j + 1;
67 new TH2F(
title.str().c_str(),
title.str().c_str(), 50, -2.5, 2.5, 1000, 0, 100);
69 title <<
"Chi2GoodHit_" <<
i + 1 <<
"-" <<
j + 1;
72 title <<
"Chi2BadHit_" <<
i + 1 <<
"-" <<
j + 1;
75 title <<
"Chi2DeltaHit_" <<
i + 1 <<
"-" <<
j + 1;
78 title <<
"Chi2NSharedHit_" <<
i + 1 <<
"-" <<
j + 1;
81 title <<
"Chi2SharedHit_" <<
i + 1 <<
"-" <<
j + 1;
85 title <<
"PullGP_X_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts";
88 title <<
"PullGP_Y_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts";
91 title <<
"PullGP_Z_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts";
95 title <<
"PullGM_X_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts";
98 title <<
"PullGM_Y_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts";
101 title <<
"PullGM_Z_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts";
105 title <<
"PullGP_X_" <<
i + 1 <<
"-" <<
j + 1 <<
"_tr";
108 title <<
"PullGP_Y_" <<
i + 1 <<
"-" <<
j + 1 <<
"_tr";
111 title <<
"PullGP_Z_" <<
i + 1 <<
"-" <<
j + 1 <<
"_tr";
115 title <<
"PullGP_X_" <<
i + 1 <<
"-" <<
j + 1 <<
"_rs";
118 title <<
"PullGP_Y_" <<
i + 1 <<
"-" <<
j + 1 <<
"_rs";
121 title <<
"PullGP_Z_" <<
i + 1 <<
"-" <<
j + 1 <<
"_rs";
124 if (((
i == 2 ||
i == 4) && (
j == 0 ||
j == 1)) || (
i == 3 ||
i == 5)) {
127 title <<
"Chi2Increment_mono_" <<
i + 1 <<
"-" <<
j + 1;
131 title <<
"PullGP_X_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts_mono";
134 title <<
"PullGP_Y_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts_mono";
137 title <<
"PullGP_Z_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts_mono";
141 title <<
"PullGM_X_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts_mono";
144 title <<
"PullGM_Y_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts_mono";
147 title <<
"PullGM_Z_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts_mono";
151 title <<
"PullGP_X_" <<
i + 1 <<
"-" <<
j + 1 <<
"_tr_mono";
154 title <<
"PullGP_Y_" <<
i + 1 <<
"-" <<
j + 1 <<
"_tr_mono";
157 title <<
"PullGP_Z_" <<
i + 1 <<
"-" <<
j + 1 <<
"_tr_mono";
161 title <<
"PullGP_X_" <<
i + 1 <<
"-" <<
j + 1 <<
"_rs_mono";
164 title <<
"PullGP_Y_" <<
i + 1 <<
"-" <<
j + 1 <<
"_rs_mono";
167 title <<
"PullGP_Z_" <<
i + 1 <<
"-" <<
j + 1 <<
"_rs_mono";
172 title <<
"Chi2Increment_stereo_" <<
i + 1 <<
"-" <<
j + 1;
176 title <<
"PullGP_X_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts_stereo";
179 title <<
"PullGP_Y_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts_stereo";
182 title <<
"PullGP_Z_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts_stereo";
186 title <<
"PullGM_X_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts_stereo";
189 title <<
"PullGM_Y_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts_stereo";
192 title <<
"PullGM_Z_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts_stereo";
196 title <<
"PullGP_X_" <<
i + 1 <<
"-" <<
j + 1 <<
"_tr_stereo";
199 title <<
"PullGP_Y_" <<
i + 1 <<
"-" <<
j + 1 <<
"_tr_stereo";
202 title <<
"PullGP_Z_" <<
i + 1 <<
"-" <<
j + 1 <<
"_tr_stereo";
206 title <<
"PullGP_X_" <<
i + 1 <<
"-" <<
j + 1 <<
"_rs_stereo";
209 title <<
"PullGP_Y_" <<
i + 1 <<
"-" <<
j + 1 <<
"_rs_stereo";
212 title <<
"PullGP_Z_" <<
i + 1 <<
"-" <<
j + 1 <<
"_rs_stereo";
216 hTotChi2Increment =
new TH1F(
"TotChi2Increment",
"TotChi2Increment", 1000, 0, 100);
217 hTotChi2GoodHit =
new TH1F(
"TotChi2GoodHit",
"TotChi2GoodHit", 1000, 0, 100);
218 hTotChi2BadHit =
new TH1F(
"TotChi2BadHit",
"TotChi2BadHit", 1000, 0, 100);
219 hTotChi2DeltaHit =
new TH1F(
"TotChi2DeltaHit",
"TotChi2DeltaHit", 1000, 0, 100);
220 hTotChi2NSharedHit =
new TH1F(
"TotChi2NSharedHit",
"TotChi2NSharedHit", 1000, 0, 100);
221 hTotChi2SharedHit =
new TH1F(
"TotChi2SharedHit",
"TotChi2SharedHit", 1000, 0, 100);
222 hProcess_vs_Chi2 =
new TH2F(
"Process_vs_Chi2",
"Process_vs_Chi2", 1000, 0, 100, 17, -0.5, 16.5);
223 hClsize_vs_Chi2 =
new TH2F(
"Clsize_vs_Chi2",
"Clsize_vs_Chi2", 1000, 0, 100, 17, -0.5, 16.5);
224 hPixClsize_vs_Chi2 =
new TH2F(
"PixClsize_vs_Chi2",
"PixClsize_vs_Chi2", 1000, 0, 100, 17, -0.5, 16.5);
225 hPrjClsize_vs_Chi2 =
new TH2F(
"PrjClsize_vs_Chi2",
"PrjClsize_vs_Chi2", 1000, 0, 100, 17, -0.5, 16.5);
226 hSt1Clsize_vs_Chi2 =
new TH2F(
"St1Clsize_vs_Chi2",
"St1Clsize_vs_Chi2", 1000, 0, 100, 17, -0.5, 16.5);
227 hSt2Clsize_vs_Chi2 =
new TH2F(
"St2Clsize_vs_Chi2",
"St2Clsize_vs_Chi2", 1000, 0, 100, 17, -0.5, 16.5);
228 hGoodHit_vs_Chi2 =
new TH2F(
"GoodHit_vs_Chi2",
"GoodHit_vs_Chi2", 10000, 0, 1000, 2, -0.5, 1.5);
229 hClusterSize =
new TH1F(
"ClusterSize",
"ClusterSize", 40, -0.5, 39.5);
230 hPixClusterSize =
new TH1F(
"PixClusterSize",
"PixClusterSize", 40, -0.5, 39.5);
231 hPrjClusterSize =
new TH1F(
"PrjClusterSize",
"PrjClusterSize", 40, -0.5, 39.5);
232 hSt1ClusterSize =
new TH1F(
"St1ClusterSize",
"St1ClusterSize", 40, -0.5, 39.5);
233 hSt2ClusterSize =
new TH1F(
"St2ClusterSize",
"St2ClusterSize", 40, -0.5, 39.5);
234 hSimHitVecSize =
new TH1F(
"hSimHitVecSize",
"hSimHitVecSize", 40, -0.5, 39.5);
235 hPixSimHitVecSize =
new TH1F(
"PixSimHitVecSize",
"PixSimHitVecSize", 40, -0.5, 39.5);
236 hPrjSimHitVecSize =
new TH1F(
"PrjSimHitVecSize",
"PrjSimHitVecSize", 40, -0.5, 39.5);
237 hSt1SimHitVecSize =
new TH1F(
"St1SimHitVecSize",
"St1SimHitVecSize", 40, -0.5, 39.5);
238 hSt2SimHitVecSize =
new TH1F(
"St2SimHitVecSize",
"St2SimHitVecSize", 40, -0.5, 39.5);
239 goodbadmerged =
new TH1F(
"goodbadmerged",
"goodbadmerged", 5, 0.5, 5.5);
240 energyLossRatio =
new TH1F(
"energyLossRatio",
"energyLossRatio", 100, 0, 1);
241 mergedPull =
new TH1F(
"mergedPull",
"mergedPull", 200, 0, 20);
242 probXgood =
new TH1F(
"probXgood",
"probXgood", 110, 0, 1.1);
243 probXbad =
new TH1F(
"probXbad",
"probXbad", 110, 0, 1.1);
244 probXdelta =
new TH1F(
"probXdelta",
"probXdelta", 110, 0, 1.1);
245 probXshared =
new TH1F(
"probXshared",
"probXshared", 110, 0, 1.1);
246 probXnoshare =
new TH1F(
"probXnoshare",
"probXnoshare", 110, 0, 1.1);
247 probYgood =
new TH1F(
"probYgood",
"probYgood", 110, 0, 1.1);
248 probYbad =
new TH1F(
"probYbad",
"probYbad", 110, 0, 1.1);
249 probYdelta =
new TH1F(
"probYdelta",
"probYdelta", 110, 0, 1.1);
250 probYshared =
new TH1F(
"probYshared",
"probYshared", 110, 0, 1.1);
251 probYnoshare =
new TH1F(
"probYnoshare",
"probYnoshare", 110, 0, 1.1);
259 LogDebug(
"TestTrackHits") <<
"new event";
285 LogTrace(
"TestTrackHits") <<
"\n*****************new trajectory********************";
288 std::vector<TrajectoryMeasurement> tmColl = it->measurements();
299 std::vector<std::pair<TrackingParticleRef, double> > tP;
301 tP = recSimColl[
track];
304 <<
" associated with quality:" << tP.begin()->second <<
" good track #"
305 << ++
yy <<
" has hits:" <<
track->numberOfValidHits() <<
"\n";
309 <<
" NOT associated to any TrackingParticle"
325 LogTrace(
"TestTrackHits") <<
"a tp is associated with fraction=" << tP.begin()->second;
327 std::vector<unsigned int> tpids;
329 LogTrace(
"TestTrackHits") <<
"tp id=" << g4T->trackId();
330 tpids.push_back(g4T->trackId());
335 for (std::vector<TrajectoryMeasurement>::iterator tm = tmColl.begin(); tm != tmColl.end(); ++tm) {
336 tchi2 += tm->estimate();
338 LogTrace(
"TestTrackHits") <<
"+++++++++++++++++new hit+++++++++++++++++";
339 CTTRHp rhit = tm->recHit();
342 TSOS state =
combiner(tm->backwardPredictedState(), tm->forwardPredictedState());
344 if (rhit->isValid() == 0 && rhit->det() !=
nullptr)
347 LogTrace(
"TestTrackHits") <<
"valid hit #" << ++
pp <<
"of hits=" <<
track->numberOfValidHits();
349 int subdetId = rhit->det()->geographicalId().subdetId();
350 DetId id = rhit->
det()->geographicalId();
351 int layerId = tTopo->
layer(
id);
352 LogTrace(
"TestTrackHits") <<
"subdetId=" << subdetId <<
" layerId=" << layerId;
354 const Surface* surf = rhit->surface();
358 double energyLoss_ = 0.;
359 unsigned int monoId = 0;
360 std::vector<double> energyLossM;
361 std::vector<double> energyLossS;
362 std::vector<PSimHit> assSimHits =
hitAssociator.associateHit(*(rhit)->
hit());
363 unsigned int simhitvecsize = assSimHits.size();
364 if (simhitvecsize == 0)
367 std::vector<unsigned int> trackIds;
370 for (std::vector<PSimHit>::const_iterator
m = assSimHits.begin();
m < assSimHits.end();
m++) {
371 unsigned int tId =
m->trackId();
372 if (
find(trackIds.begin(), trackIds.end(), tId) == trackIds.end())
373 trackIds.push_back(tId);
374 if (
m->energyLoss() > energyLoss_) {
376 energyLoss_ =
m->energyLoss();
378 if (dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())) {
380 monoId =
m->detUnitId();
381 if (monoId ==
m->detUnitId()) {
382 energyLossM.push_back(
m->energyLoss());
384 energyLossS.push_back(
m->energyLoss());
388 energyLossM.push_back(
m->energyLoss());
405 double chi2increment = tm->estimate();
406 LogTrace(
"TestTrackHits") <<
"tm->estimate()=" << tm->estimate();
408 title <<
"Chi2Increment_" << subdetId <<
"-" << layerId;
411 title <<
"Chi2IncrementVsEta_" << subdetId <<
"-" << layerId;
417 bool mergedhit =
false;
418 if (dynamic_cast<const SiPixelRecHit*>(rhit->hit())) {
419 clustersize = ((
const SiPixelRecHit*)(rhit->hit()))->cluster()->size();
423 if (simhitvecsize > 1)
425 }
else if (dynamic_cast<const SiStripRecHit2D*>(rhit->hit())) {
426 clustersize = ((
const SiStripRecHit2D*)(rhit->hit()))->cluster()->amplitudes().size();
430 if (simhitvecsize > 1)
432 }
else if (dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())) {
435 if (clsize1 > clsize2)
436 clustersize = clsize1;
438 clustersize = clsize2;
442 if (simhitvecsize > 2)
444 }
else if (dynamic_cast<const ProjectedSiStripRecHit2D*>(rhit->hit())) {
449 if (simhitvecsize > 1)
461 std::vector<SimHitIdpr> simTrackIds =
hitAssociator.associateHitId(*(rhit)->
hit());
462 bool goodhit =
false;
463 for (
size_t j = 0;
j < simTrackIds.size();
j++) {
464 LogTrace(
"TestTrackHits") <<
"hit id=" << simTrackIds[
j].first;
465 for (
size_t jj = 0;
jj < tpids.size();
jj++) {
466 if (simTrackIds[
j].
first == tpids[
jj])
472 bool ioniOnly =
true;
473 const SiPixelRecHit* pix = dynamic_cast<const SiPixelRecHit*>(rhit->hit());
475 if (energyLossM.size() > 1 && energyLossS.size() <= 1) {
476 sort(energyLossM.begin(), energyLossM.end(), greater<double>());
478 }
else if (energyLossS.size() > 1 && energyLossM.size() <= 1) {
479 sort(energyLossS.begin(), energyLossS.end(), greater<double>());
481 }
else if (energyLossS.size() > 1 && energyLossM.size() > 1) {
482 sort(energyLossM.begin(), energyLossM.end(), greater<double>());
483 sort(energyLossS.begin(), energyLossS.end(), greater<double>());
484 energyLossM[1] / energyLossM[0] > energyLossS[1] / energyLossS[0]
490 LogVerbatim(
"TestTrackHits") <<
"MERGED HIT" << std::endl;
491 unsigned int idc = 0;
492 for (
size_t jj = 0;
jj < tpids.size();
jj++) {
493 idc +=
std::count(trackIds.begin(), trackIds.end(), tpids[
jj]);
495 if (idc == trackIds.size()) {
498 for (std::vector<PSimHit>::const_iterator
m = assSimHits.begin() + 1;
m < assSimHits.end();
m++) {
499 if ((
m->processType() != 7 &&
m->processType() != 8 &&
m->processType() != 9) &&
500 abs(
m->particleType()) != 11) {
505 if (ioniOnly && !shared) {
507 title <<
"Chi2DeltaHit_" << subdetId <<
"-" << layerId;
514 }
else if (!ioniOnly && !shared) {
516 title <<
"Chi2NSharedHit_" << subdetId <<
"-" << layerId;
525 title <<
"Chi2SharedHit_" << subdetId <<
"-" << layerId;
534 for (std::vector<PSimHit>::const_iterator
m = assSimHits.begin();
m < assSimHits.end();
m++) {
535 unsigned int tId =
m->trackId();
536 LogVerbatim(
"TestTrackHits") <<
"component with id=" << tId <<
" eLoss=" <<
m->energyLoss()
537 <<
" pType=" <<
m->processType();
538 if (
find(tpids.begin(), tpids.end(), tId) == tpids.end())
540 if (
m->processType() == 2) {
544 LogVerbatim(
"TestTrackHits") << gpr <<
" " << gps <<
" " << ger;
545 ROOT::Math::SVector<double, 3>
delta;
552 LogVerbatim(
"TestTrackHits") <<
"hit pull=" << mpull;
560 title <<
"Chi2GoodHit_" << subdetId <<
"-" << layerId;
571 title <<
"Chi2BadHit_" << subdetId <<
"-" << layerId;
584 if (dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())) {
585 if (simhitvecsize > 2 && goodhit) {
586 if (ioniOnly && !shared)
588 else if (!ioniOnly && !shared)
594 double rechitmatchedx = rhit->localPosition().x();
595 double rechitmatchedy = rhit->localPosition().y();
596 double mindist = 999999;
598 std::pair<LocalPoint, LocalVector> closestPair;
600 const BoundPlane& plane = (rhit)->det()->surface();
601 for (std::vector<PSimHit>::const_iterator
m = assSimHits.begin();
m < assSimHits.end();
m++) {
603 std::pair<LocalPoint, LocalVector> hitPair =
projectHit((*
m), stripDet, plane);
604 distx = fabs(rechitmatchedx - hitPair.first.x());
605 disty = fabs(rechitmatchedy - hitPair.first.y());
606 double dist = distx * distx + disty * disty;
607 if (
sqrt(dist) < mindist) {
609 closestPair = hitPair;
612 shitLPos = closestPair.first;
613 shitLMom = closestPair.second;
615 if (simhitvecsize > 1 && goodhit) {
616 if (ioniOnly && !shared)
618 else if (!ioniOnly && !shared)
636 GlobalError rhitGPEr = (rhit)->globalPositionError();
638 LogVerbatim(
"TestTrackHits") <<
"assSimHits.size()=" << assSimHits.size();
639 LogVerbatim(
"TestTrackHits") <<
"tsos globalPos =" << tsosGPos;
640 LogVerbatim(
"TestTrackHits") <<
"sim hit globalPos=" << shitGPos;
641 LogVerbatim(
"TestTrackHits") <<
"rec hit globalPos=" << rhitGPos;
642 LogVerbatim(
"TestTrackHits") <<
"geographicalId =" << rhit->det()->geographicalId().rawId();
646 if (rhit->detUnit())
LogTrace(
"TestTrackHits") <<
"rhit->detUnit()->geographicalId()="
647 << rhit->detUnit()->geographicalId().rawId() ;
648 LogTrace(
"TestTrackHits") <<
"rhit->det()->surface().position()="
649 << rhit->det()->surface().position() ;
650 if (rhit->detUnit())
LogTrace(
"TestTrackHits") <<
"rhit->detUnit()->surface().position()="
651 << rhit->detUnit()->surface().position() ;
652 LogTrace(
"TestTrackHits") <<
"rhit->det()->position()=" << rhit->det()->position() ;
653 if (rhit->detUnit())
LogTrace(
"TestTrackHits") <<
"rhit->detUnit()->position()="
654 << rhit->detUnit()->position() ;
655 LogTrace(
"TestTrackHits") <<
"rhit->det()->surface().bounds().length()="
656 << rhit->det()->surface().bounds().length() ;
657 if (rhit->detUnit())
LogTrace(
"TestTrackHits") <<
"rhit->detUnit()->surface().bounds().length()="
658 << rhit->detUnit()->surface().bounds().length() ;
659 LogTrace(
"TestTrackHits") <<
"rhit->det()->surface().bounds().width()="
660 << rhit->det()->surface().bounds().width() ;
661 if (rhit->detUnit())
LogTrace(
"TestTrackHits") <<
"rhit->detUnit()->surface().bounds().width()="
662 << rhit->detUnit()->surface().bounds().width() ;
663 LogTrace(
"TestTrackHits") <<
"rhit->det()->surface().bounds().thickness()="
664 << rhit->det()->surface().bounds().thickness() ;
665 if (rhit->detUnit())
LogTrace(
"TestTrackHits") <<
"rhit->detUnit()->surface().bounds().thickness()="
666 << rhit->detUnit()->surface().bounds().thickness() ;
669 double pullGPX_rs = (rhitGPos.
x() - shitGPos.
x()) /
sqrt(rhitGPEr.
cxx());
670 double pullGPY_rs = (rhitGPos.
y() - shitGPos.
y()) /
sqrt(rhitGPEr.
cyy());
671 double pullGPZ_rs = (rhitGPos.
z() - shitGPos.
z()) /
sqrt(rhitGPEr.
czz());
679 title <<
"PullGP_X_" << subdetId <<
"-" << layerId <<
"_rs";
682 title <<
"PullGP_Y_" << subdetId <<
"-" << layerId <<
"_rs";
685 title <<
"PullGP_Z_" << subdetId <<
"-" << layerId <<
"_rs";
688 double pullGPX_tr = (tsosGPos.
x() - rhitGPos.
x()) /
sqrt(tsosGPEr.
cxx() + rhitGPEr.
cxx());
689 double pullGPY_tr = (tsosGPos.
y() - rhitGPos.
y()) /
sqrt(tsosGPEr.
cyy() + rhitGPEr.
cyy());
690 double pullGPZ_tr = (tsosGPos.
z() - rhitGPos.
z()) /
sqrt(tsosGPEr.
czz() + rhitGPEr.
czz());
698 title <<
"PullGP_X_" << subdetId <<
"-" << layerId <<
"_tr";
701 title <<
"PullGP_Y_" << subdetId <<
"-" << layerId <<
"_tr";
704 title <<
"PullGP_Z_" << subdetId <<
"-" << layerId <<
"_tr";
707 double pullGPX_ts = (tsosGPos.
x() - shitGPos.
x()) /
sqrt(tsosGPEr.
cxx());
708 double pullGPY_ts = (tsosGPos.
y() - shitGPos.
y()) /
sqrt(tsosGPEr.
cyy());
709 double pullGPZ_ts = (tsosGPos.
z() - shitGPos.
z()) /
sqrt(tsosGPEr.
czz());
717 title <<
"PullGP_X_" << subdetId <<
"-" << layerId <<
"_ts";
720 title <<
"PullGP_Y_" << subdetId <<
"-" << layerId <<
"_ts";
723 title <<
"PullGP_Z_" << subdetId <<
"-" << layerId <<
"_ts";
726 double pullGMX_ts = (tsosGMom.
x() - shitGMom.
x()) /
sqrt(tsosGMEr.
cxx());
727 double pullGMY_ts = (tsosGMom.
y() - shitGMom.
y()) /
sqrt(tsosGMEr.
cyy());
728 double pullGMZ_ts = (tsosGMom.
z() - shitGMom.
z()) /
sqrt(tsosGMEr.
czz());
736 title <<
"PullGM_X_" << subdetId <<
"-" << layerId <<
"_ts";
739 title <<
"PullGM_Y_" << subdetId <<
"-" << layerId <<
"_ts";
742 title <<
"PullGM_Z_" << subdetId <<
"-" << layerId <<
"_ts";
744 if (dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())) {
747 LogTrace(
"TestTrackHits") <<
"MONO HIT";
748 auto m = dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())->monoHit();
750 if (tMonoHit ==
nullptr)
752 vector<PSimHit> assMonoSimHits =
hitAssociator.associateHit(*tMonoHit->hit());
753 if (assMonoSimHits.empty())
755 const PSimHit sMonoHit = *(assMonoSimHits.begin());
756 const Surface* monoSurf = &(tMonoHit->det()->surface());
757 if (monoSurf ==
nullptr)
776 GlobalPoint monoRhitGPos = tMonoHit->globalPosition();
777 GlobalError monoRhitGPEr = tMonoHit->globalPositionError();
779 double pullGPX_rs_mono = (monoRhitGPos.
x() - monoShitGPos.
x()) /
sqrt(monoRhitGPEr.
cxx());
780 double pullGPY_rs_mono = (monoRhitGPos.
y() - monoShitGPos.
y()) /
sqrt(monoRhitGPEr.
cyy());
781 double pullGPZ_rs_mono = (monoRhitGPos.
z() - monoShitGPos.
z()) /
sqrt(monoRhitGPEr.
czz());
790 title <<
"Chi2Increment_mono_" << subdetId <<
"-" << layerId;
794 title <<
"PullGP_X_" << subdetId <<
"-" << layerId <<
"_rs_mono";
797 title <<
"PullGP_Y_" << subdetId <<
"-" << layerId <<
"_rs_mono";
800 title <<
"PullGP_Z_" << subdetId <<
"-" << layerId <<
"_rs_mono";
803 double pullGPX_tr_mono = (monoTsosGPos.
x() - monoRhitGPos.
x()) /
sqrt(monoTsosGPEr.
cxx() + monoRhitGPEr.
cxx());
804 double pullGPY_tr_mono = (monoTsosGPos.
y() - monoRhitGPos.
y()) /
sqrt(monoTsosGPEr.
cyy() + monoRhitGPEr.
cyy());
805 double pullGPZ_tr_mono = (monoTsosGPos.
z() - monoRhitGPos.
z()) /
sqrt(monoTsosGPEr.
czz() + monoRhitGPEr.
czz());
811 title <<
"PullGP_X_" << subdetId <<
"-" << layerId <<
"_tr_mono";
814 title <<
"PullGP_Y_" << subdetId <<
"-" << layerId <<
"_tr_mono";
817 title <<
"PullGP_Z_" << subdetId <<
"-" << layerId <<
"_tr_mono";
820 double pullGPX_ts_mono = (monoTsosGPos.
x() - monoShitGPos.
x()) /
sqrt(monoTsosGPEr.
cxx());
821 double pullGPY_ts_mono = (monoTsosGPos.
y() - monoShitGPos.
y()) /
sqrt(monoTsosGPEr.
cyy());
822 double pullGPZ_ts_mono = (monoTsosGPos.
z() - monoShitGPos.
z()) /
sqrt(monoTsosGPEr.
czz());
828 title <<
"PullGP_X_" << subdetId <<
"-" << layerId <<
"_ts_mono";
831 title <<
"PullGP_Y_" << subdetId <<
"-" << layerId <<
"_ts_mono";
834 title <<
"PullGP_Z_" << subdetId <<
"-" << layerId <<
"_ts_mono";
837 double pullGMX_ts_mono = (monoTsosGMom.
x() - monoShitGMom.
x()) /
sqrt(monoTsosGMEr.
cxx());
838 double pullGMY_ts_mono = (monoTsosGMom.
y() - monoShitGMom.
y()) /
sqrt(monoTsosGMEr.
cyy());
839 double pullGMZ_ts_mono = (monoTsosGMom.
z() - monoShitGMom.
z()) /
sqrt(monoTsosGMEr.
czz());
845 title <<
"PullGM_X_" << subdetId <<
"-" << layerId <<
"_ts_mono";
848 title <<
"PullGM_Y_" << subdetId <<
"-" << layerId <<
"_ts_mono";
851 title <<
"PullGM_Z_" << subdetId <<
"-" << layerId <<
"_ts_mono";
855 LogTrace(
"TestTrackHits") <<
"STEREO HIT";
856 auto s = dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())->stereoHit();
858 if (tStereoHit ==
nullptr)
860 vector<PSimHit> assStereoSimHits =
hitAssociator.associateHit(*tStereoHit->hit());
861 if (assStereoSimHits.empty())
863 const PSimHit sStereoHit = *(assStereoSimHits.begin());
864 const Surface* stereoSurf = &(tStereoHit->det()->surface());
865 if (stereoSurf ==
nullptr)
868 if (stereoState.
isValid() == 0)
884 GlobalPoint stereoRhitGPos = tStereoHit->globalPosition();
885 GlobalError stereoRhitGPEr = tStereoHit->globalPositionError();
891 title <<
"Chi2Increment_stereo_" << subdetId <<
"-" << layerId;
894 double pullGPX_rs_stereo = (stereoRhitGPos.
x() - stereoShitGPos.
x()) /
sqrt(stereoRhitGPEr.
cxx());
895 double pullGPY_rs_stereo = (stereoRhitGPos.
y() - stereoShitGPos.
y()) /
sqrt(stereoRhitGPEr.
cyy());
896 double pullGPZ_rs_stereo = (stereoRhitGPos.
z() - stereoShitGPos.
z()) /
sqrt(stereoRhitGPEr.
czz());
902 title <<
"PullGP_X_" << subdetId <<
"-" << layerId <<
"_rs_stereo";
905 title <<
"PullGP_Y_" << subdetId <<
"-" << layerId <<
"_rs_stereo";
908 title <<
"PullGP_Z_" << subdetId <<
"-" << layerId <<
"_rs_stereo";
911 double pullGPX_tr_stereo =
912 (stereoTsosGPos.
x() - stereoRhitGPos.
x()) /
sqrt(stereoTsosGPEr.
cxx() + stereoRhitGPEr.
cxx());
913 double pullGPY_tr_stereo =
914 (stereoTsosGPos.
y() - stereoRhitGPos.
y()) /
sqrt(stereoTsosGPEr.
cyy() + stereoRhitGPEr.
cyy());
915 double pullGPZ_tr_stereo =
916 (stereoTsosGPos.
z() - stereoRhitGPos.
z()) /
sqrt(stereoTsosGPEr.
czz() + stereoRhitGPEr.
czz());
922 title <<
"PullGP_X_" << subdetId <<
"-" << layerId <<
"_tr_stereo";
925 title <<
"PullGP_Y_" << subdetId <<
"-" << layerId <<
"_tr_stereo";
928 title <<
"PullGP_Z_" << subdetId <<
"-" << layerId <<
"_tr_stereo";
931 double pullGPX_ts_stereo = (stereoTsosGPos.
x() - stereoShitGPos.
x()) /
sqrt(stereoTsosGPEr.
cxx());
932 double pullGPY_ts_stereo = (stereoTsosGPos.
y() - stereoShitGPos.
y()) /
sqrt(stereoTsosGPEr.
cyy());
933 double pullGPZ_ts_stereo = (stereoTsosGPos.
z() - stereoShitGPos.
z()) /
sqrt(stereoTsosGPEr.
czz());
939 title <<
"PullGP_X_" << subdetId <<
"-" << layerId <<
"_ts_stereo";
942 title <<
"PullGP_Y_" << subdetId <<
"-" << layerId <<
"_ts_stereo";
945 title <<
"PullGP_Z_" << subdetId <<
"-" << layerId <<
"_ts_stereo";
948 double pullGMX_ts_stereo = (stereoTsosGMom.
x() - stereoShitGMom.
x()) /
sqrt(stereoTsosGMEr.
cxx());
949 double pullGMY_ts_stereo = (stereoTsosGMom.
y() - stereoShitGMom.
y()) /
sqrt(stereoTsosGMEr.
cyy());
950 double pullGMZ_ts_stereo = (stereoTsosGMom.
z() - stereoShitGMom.
z()) /
sqrt(stereoTsosGMEr.
czz());
956 title <<
"PullGM_X_" << subdetId <<
"-" << layerId <<
"_ts_stereo";
959 title <<
"PullGM_Y_" << subdetId <<
"-" << layerId <<
"_ts_stereo";
962 title <<
"PullGM_Z_" << subdetId <<
"-" << layerId <<
"_ts_stereo";
966 LogTrace(
"TestTrackHits") <<
"traj chi2=" << tchi2;
970 LogTrace(
"TestTrackHits") <<
"end of event: processd hits=" << evtHits;
975 TDirectory* chi2d =
file->mkdir(
"Chi2Increment");
977 TDirectory* gp_ts =
file->mkdir(
"GP_TSOS-SimHit");
978 TDirectory* gm_ts =
file->mkdir(
"GM_TSOS-SimHit");
979 TDirectory* gp_tr =
file->mkdir(
"GP_TSOS-RecHit");
980 TDirectory* gp_rs =
file->mkdir(
"GP_RecHit-SimHit");
982 TDirectory* gp_tsx = gp_ts->mkdir(
"X");
983 TDirectory* gp_tsy = gp_ts->mkdir(
"Y");
984 TDirectory* gp_tsz = gp_ts->mkdir(
"Z");
985 TDirectory* gm_tsx = gm_ts->mkdir(
"X");
986 TDirectory* gm_tsy = gm_ts->mkdir(
"Y");
987 TDirectory* gm_tsz = gm_ts->mkdir(
"Z");
988 TDirectory* gp_trx = gp_tr->mkdir(
"X");
989 TDirectory* gp_try = gp_tr->mkdir(
"Y");
990 TDirectory* gp_trz = gp_tr->mkdir(
"Z");
991 TDirectory* gp_rsx = gp_rs->mkdir(
"X");
992 TDirectory* gp_rsy = gp_rs->mkdir(
"Y");
993 TDirectory* gp_rsz = gp_rs->mkdir(
"Z");
995 TDirectory* gp_tsx_mono = gp_ts->mkdir(
"MONOX");
996 TDirectory* gp_tsy_mono = gp_ts->mkdir(
"MONOY");
997 TDirectory* gp_tsz_mono = gp_ts->mkdir(
"MONOZ");
998 TDirectory* gm_tsx_mono = gm_ts->mkdir(
"MONOX");
999 TDirectory* gm_tsy_mono = gm_ts->mkdir(
"MONOY");
1000 TDirectory* gm_tsz_mono = gm_ts->mkdir(
"MONOZ");
1001 TDirectory* gp_trx_mono = gp_tr->mkdir(
"MONOX");
1002 TDirectory* gp_try_mono = gp_tr->mkdir(
"MONOY");
1003 TDirectory* gp_trz_mono = gp_tr->mkdir(
"MONOZ");
1004 TDirectory* gp_rsx_mono = gp_rs->mkdir(
"MONOX");
1005 TDirectory* gp_rsy_mono = gp_rs->mkdir(
"MONOY");
1006 TDirectory* gp_rsz_mono = gp_rs->mkdir(
"MONOZ");
1008 TDirectory* gp_tsx_stereo = gp_ts->mkdir(
"STEREOX");
1009 TDirectory* gp_tsy_stereo = gp_ts->mkdir(
"STEREOY");
1010 TDirectory* gp_tsz_stereo = gp_ts->mkdir(
"STEREOZ");
1011 TDirectory* gm_tsx_stereo = gm_ts->mkdir(
"STEREOX");
1012 TDirectory* gm_tsy_stereo = gm_ts->mkdir(
"STEREOY");
1013 TDirectory* gm_tsz_stereo = gm_ts->mkdir(
"STEREOZ");
1014 TDirectory* gp_trx_stereo = gp_tr->mkdir(
"STEREOX");
1015 TDirectory* gp_try_stereo = gp_tr->mkdir(
"STEREOY");
1016 TDirectory* gp_trz_stereo = gp_tr->mkdir(
"STEREOZ");
1017 TDirectory* gp_rsx_stereo = gp_rs->mkdir(
"STEREOX");
1018 TDirectory* gp_rsy_stereo = gp_rs->mkdir(
"STEREOY");
1019 TDirectory* gp_rsz_stereo = gp_rs->mkdir(
"STEREOZ");
1058 for (
int i = 0;
i != 6;
i++)
1059 for (
int j = 0;
j != 9;
j++) {
1060 if (
i == 0 &&
j > 2)
1062 if (
i == 1 &&
j > 1)
1064 if (
i == 2 &&
j > 3)
1066 if (
i == 3 &&
j > 2)
1068 if (
i == 4 &&
j > 5)
1070 if (
i == 5 &&
j > 8)
1074 title <<
"Chi2Increment_" <<
i + 1 <<
"-" <<
j + 1;
1077 title <<
"Chi2IncrementVsEta_" <<
i + 1 <<
"-" <<
j + 1;
1080 title <<
"Chi2GoodHit_" <<
i + 1 <<
"-" <<
j + 1;
1083 title <<
"Chi2BadHit_" <<
i + 1 <<
"-" <<
j + 1;
1086 title <<
"Chi2DeltaHit_" <<
i + 1 <<
"-" <<
j + 1;
1089 title <<
"Chi2NSharedHit_" <<
i + 1 <<
"-" <<
j + 1;
1092 title <<
"Chi2SharedHit_" <<
i + 1 <<
"-" <<
j + 1;
1098 title <<
"PullGP_X_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts";
1102 title <<
"PullGP_Y_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts";
1106 title <<
"PullGP_Z_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts";
1112 title <<
"PullGM_X_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts";
1116 title <<
"PullGM_Y_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts";
1120 title <<
"PullGM_Z_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts";
1126 title <<
"PullGP_X_" <<
i + 1 <<
"-" <<
j + 1 <<
"_tr";
1130 title <<
"PullGP_Y_" <<
i + 1 <<
"-" <<
j + 1 <<
"_tr";
1134 title <<
"PullGP_Z_" <<
i + 1 <<
"-" <<
j + 1 <<
"_tr";
1140 title <<
"PullGP_X_" <<
i + 1 <<
"-" <<
j + 1 <<
"_rs";
1144 title <<
"PullGP_Y_" <<
i + 1 <<
"-" <<
j + 1 <<
"_rs";
1148 title <<
"PullGP_Z_" <<
i + 1 <<
"-" <<
j + 1 <<
"_rs";
1151 if (((
i == 2 ||
i == 4) && (
j == 0 ||
j == 1)) || (
i == 3 ||
i == 5)) {
1154 title <<
"Chi2Increment_mono_" <<
i + 1 <<
"-" <<
j + 1;
1157 title <<
"Chi2Increment_stereo_" <<
i + 1 <<
"-" <<
j + 1;
1163 title <<
"PullGP_X_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts_mono";
1167 title <<
"PullGP_Y_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts_mono";
1171 title <<
"PullGP_Z_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts_mono";
1177 title <<
"PullGM_X_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts_mono";
1181 title <<
"PullGM_Y_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts_mono";
1185 title <<
"PullGM_Z_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts_mono";
1191 title <<
"PullGP_X_" <<
i + 1 <<
"-" <<
j + 1 <<
"_tr_mono";
1195 title <<
"PullGP_Y_" <<
i + 1 <<
"-" <<
j + 1 <<
"_tr_mono";
1199 title <<
"PullGP_Z_" <<
i + 1 <<
"-" <<
j + 1 <<
"_tr_mono";
1205 title <<
"PullGP_X_" <<
i + 1 <<
"-" <<
j + 1 <<
"_rs_mono";
1209 title <<
"PullGP_Y_" <<
i + 1 <<
"-" <<
j + 1 <<
"_rs_mono";
1213 title <<
"PullGP_Z_" <<
i + 1 <<
"-" <<
j + 1 <<
"_rs_mono";
1218 gp_tsx_stereo->cd();
1220 title <<
"PullGP_X_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts_stereo";
1222 gp_tsy_stereo->cd();
1224 title <<
"PullGP_Y_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts_stereo";
1226 gp_tsz_stereo->cd();
1228 title <<
"PullGP_Z_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts_stereo";
1232 gm_tsx_stereo->cd();
1234 title <<
"PullGM_X_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts_stereo";
1236 gm_tsy_stereo->cd();
1238 title <<
"PullGM_Y_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts_stereo";
1240 gm_tsz_stereo->cd();
1242 title <<
"PullGM_Z_" <<
i + 1 <<
"-" <<
j + 1 <<
"_ts_stereo";
1246 gp_trx_stereo->cd();
1248 title <<
"PullGP_X_" <<
i + 1 <<
"-" <<
j + 1 <<
"_tr_stereo";
1250 gp_try_stereo->cd();
1252 title <<
"PullGP_Y_" <<
i + 1 <<
"-" <<
j + 1 <<
"_tr_stereo";
1254 gp_trz_stereo->cd();
1256 title <<
"PullGP_Z_" <<
i + 1 <<
"-" <<
j + 1 <<
"_tr_stereo";
1260 gp_rsx_stereo->cd();
1262 title <<
"PullGP_X_" <<
i + 1 <<
"-" <<
j + 1 <<
"_rs_stereo";
1264 gp_rsy_stereo->cd();
1266 title <<
"PullGP_Y_" <<
i + 1 <<
"-" <<
j + 1 <<
"_rs_stereo";
1268 gp_rsz_stereo->cd();
1270 title <<
"PullGP_Z_" <<
i + 1 <<
"-" <<
j + 1 <<
"_rs_stereo";
1285 LocalPoint localHit = plane.toLocal(globalpos);
1302 return std::pair<LocalPoint, LocalVector>(projectedPos, localStripDir);
1305 template <
unsigned int D>
1309 VecD
r = asSVector<D>(rhit->parameters()) -
me.measuredParameters<
D>(*rhit);
1311 SMatDD
R = asSMatrix<D>(rhit->parametersError()) +
me.measuredError<
D>(*rhit);
1313 return ROOT::Math::Similarity(
r,
R);