8 #include <TDirectory.h>
49 file =
new TFile(
out.c_str(),
"recreate");
50 for (
int i=0;
i!=6;
i++)
51 for (
int j=0;
j!=9;
j++){
52 if (
i==0 &&
j>2)
break;
53 if (
i==1 &&
j>1)
break;
54 if (
i==2 &&
j>3)
break;
55 if (
i==3 &&
j>2)
break;
56 if (
i==4 &&
j>5)
break;
57 if (
i==5 &&
j>8)
break;
60 title <<
"Chi2Increment_" <<
i+1 <<
"-" <<
j+1 ;
63 title <<
"Chi2IncrementVsEta_" <<
i+1 <<
"-" << j+1 ;
66 title <<
"Chi2GoodHit_" <<
i+1 <<
"-" << j+1 ;
69 title <<
"Chi2BadHit_" <<
i+1 <<
"-" << j+1 ;
72 title <<
"Chi2DeltaHit_" <<
i+1 <<
"-" << j+1 ;
75 title <<
"Chi2NSharedHit_" <<
i+1 <<
"-" << j+1 ;
78 title <<
"Chi2SharedHit_" <<
i+1 <<
"-" << j+1 ;
82 title <<
"PullGP_X_" <<
i+1 <<
"-" << j+1 <<
"_ts";
85 title <<
"PullGP_Y_" <<
i+1 <<
"-" << j+1 <<
"_ts";
88 title <<
"PullGP_Z_" <<
i+1 <<
"-" << j+1 <<
"_ts";
92 title <<
"PullGM_X_" <<
i+1 <<
"-" << j+1 <<
"_ts";
95 title <<
"PullGM_Y_" <<
i+1 <<
"-" << j+1 <<
"_ts";
98 title <<
"PullGM_Z_" <<
i+1 <<
"-" << j+1 <<
"_ts";
102 title <<
"PullGP_X_" <<
i+1 <<
"-" << j+1 <<
"_tr";
105 title <<
"PullGP_Y_" <<
i+1 <<
"-" << j+1 <<
"_tr";
108 title <<
"PullGP_Z_" <<
i+1 <<
"-" << j+1 <<
"_tr";
112 title <<
"PullGP_X_" <<
i+1 <<
"-" << j+1 <<
"_rs";
115 title <<
"PullGP_Y_" <<
i+1 <<
"-" << j+1 <<
"_rs";
118 title <<
"PullGP_Z_" <<
i+1 <<
"-" << j+1 <<
"_rs";
121 if ( ((
i==2||
i==4)&&(j==0||j==1)) || (
i==3||
i==5) ){
124 title <<
"Chi2Increment_mono_" <<
i+1 <<
"-" << j+1 ;
128 title <<
"PullGP_X_" <<
i+1 <<
"-" << j+1 <<
"_ts_mono";
131 title <<
"PullGP_Y_" <<
i+1 <<
"-" << j+1 <<
"_ts_mono";
134 title <<
"PullGP_Z_" <<
i+1 <<
"-" << j+1 <<
"_ts_mono";
138 title <<
"PullGM_X_" <<
i+1 <<
"-" << j+1 <<
"_ts_mono";
141 title <<
"PullGM_Y_" <<
i+1 <<
"-" << j+1 <<
"_ts_mono";
144 title <<
"PullGM_Z_" <<
i+1 <<
"-" << j+1 <<
"_ts_mono";
148 title <<
"PullGP_X_" <<
i+1 <<
"-" << j+1 <<
"_tr_mono";
151 title <<
"PullGP_Y_" <<
i+1 <<
"-" << j+1 <<
"_tr_mono";
154 title <<
"PullGP_Z_" <<
i+1 <<
"-" << j+1 <<
"_tr_mono";
158 title <<
"PullGP_X_" <<
i+1 <<
"-" << j+1 <<
"_rs_mono";
161 title <<
"PullGP_Y_" <<
i+1 <<
"-" << j+1 <<
"_rs_mono";
164 title <<
"PullGP_Z_" <<
i+1 <<
"-" << j+1 <<
"_rs_mono";
169 title <<
"Chi2Increment_stereo_" <<
i+1 <<
"-" << j+1 ;
173 title <<
"PullGP_X_" <<
i+1 <<
"-" << j+1 <<
"_ts_stereo";
176 title <<
"PullGP_Y_" <<
i+1 <<
"-" << j+1 <<
"_ts_stereo";
179 title <<
"PullGP_Z_" <<
i+1 <<
"-" << j+1 <<
"_ts_stereo";
183 title <<
"PullGM_X_" <<
i+1 <<
"-" << j+1 <<
"_ts_stereo";
186 title <<
"PullGM_Y_" <<
i+1 <<
"-" << j+1 <<
"_ts_stereo";
189 title <<
"PullGM_Z_" <<
i+1 <<
"-" << j+1 <<
"_ts_stereo";
193 title <<
"PullGP_X_" <<
i+1 <<
"-" << j+1 <<
"_tr_stereo";
196 title <<
"PullGP_Y_" <<
i+1 <<
"-" << j+1 <<
"_tr_stereo";
199 title <<
"PullGP_Z_" <<
i+1 <<
"-" << j+1 <<
"_tr_stereo";
203 title <<
"PullGP_X_" <<
i+1 <<
"-" << j+1 <<
"_rs_stereo";
206 title <<
"PullGP_Y_" <<
i+1 <<
"-" << j+1 <<
"_rs_stereo";
209 title <<
"PullGP_Z_" <<
i+1 <<
"-" << j+1 <<
"_rs_stereo";
214 hTotChi2GoodHit =
new TH1F(
"TotChi2GoodHit",
"TotChi2GoodHit",1000,0,100);
215 hTotChi2BadHit =
new TH1F(
"TotChi2BadHit",
"TotChi2BadHit",1000,0,100);
216 hTotChi2DeltaHit =
new TH1F(
"TotChi2DeltaHit",
"TotChi2DeltaHit",1000,0,100);
219 hProcess_vs_Chi2 =
new TH2F(
"Process_vs_Chi2",
"Process_vs_Chi2",1000,0,100,17,-0.5,16.5);
220 hClsize_vs_Chi2 =
new TH2F(
"Clsize_vs_Chi2",
"Clsize_vs_Chi2",1000,0,100,17,-0.5,16.5);
221 hPixClsize_vs_Chi2=
new TH2F(
"PixClsize_vs_Chi2",
"PixClsize_vs_Chi2",1000,0,100,17,-0.5,16.5);
222 hPrjClsize_vs_Chi2=
new TH2F(
"PrjClsize_vs_Chi2",
"PrjClsize_vs_Chi2",1000,0,100,17,-0.5,16.5);
223 hSt1Clsize_vs_Chi2=
new TH2F(
"St1Clsize_vs_Chi2",
"St1Clsize_vs_Chi2",1000,0,100,17,-0.5,16.5);
224 hSt2Clsize_vs_Chi2=
new TH2F(
"St2Clsize_vs_Chi2",
"St2Clsize_vs_Chi2",1000,0,100,17,-0.5,16.5);
225 hGoodHit_vs_Chi2 =
new TH2F(
"GoodHit_vs_Chi2",
"GoodHit_vs_Chi2",10000,0,1000,2,-0.5,1.5);
226 hClusterSize =
new TH1F(
"ClusterSize",
"ClusterSize",40,-0.5,39.5);
227 hPixClusterSize =
new TH1F(
"PixClusterSize",
"PixClusterSize",40,-0.5,39.5);
228 hPrjClusterSize =
new TH1F(
"PrjClusterSize",
"PrjClusterSize",40,-0.5,39.5);
229 hSt1ClusterSize =
new TH1F(
"St1ClusterSize",
"St1ClusterSize",40,-0.5,39.5);
230 hSt2ClusterSize =
new TH1F(
"St2ClusterSize",
"St2ClusterSize",40,-0.5,39.5);
231 hSimHitVecSize =
new TH1F(
"hSimHitVecSize",
"hSimHitVecSize",40,-0.5,39.5);
232 hPixSimHitVecSize =
new TH1F(
"PixSimHitVecSize",
"PixSimHitVecSize",40,-0.5,39.5);
233 hPrjSimHitVecSize =
new TH1F(
"PrjSimHitVecSize",
"PrjSimHitVecSize",40,-0.5,39.5);
234 hSt1SimHitVecSize =
new TH1F(
"St1SimHitVecSize",
"St1SimHitVecSize",40,-0.5,39.5);
235 hSt2SimHitVecSize =
new TH1F(
"St2SimHitVecSize",
"St2SimHitVecSize",40,-0.5,39.5);
236 goodbadmerged =
new TH1F(
"goodbadmerged",
"goodbadmerged",5,0.5,5.5);
237 energyLossRatio =
new TH1F(
"energyLossRatio",
"energyLossRatio",100,0,1);
238 mergedPull =
new TH1F(
"mergedPull",
"mergedPull",200,0,20);
239 probXgood =
new TH1F(
"probXgood",
"probXgood",110,0,1.1);
240 probXbad =
new TH1F(
"probXbad",
"probXbad",110,0,1.1);
241 probXdelta =
new TH1F(
"probXdelta",
"probXdelta",110,0,1.1);
242 probXshared =
new TH1F(
"probXshared",
"probXshared",110,0,1.1);
243 probXnoshare =
new TH1F(
"probXnoshare",
"probXnoshare",110,0,1.1);
244 probYgood =
new TH1F(
"probYgood",
"probYgood",110,0,1.1);
245 probYbad =
new TH1F(
"probYbad",
"probYbad",110,0,1.1);
246 probYdelta =
new TH1F(
"probYdelta",
"probYdelta",110,0,1.1);
247 probYshared =
new TH1F(
"probYshared",
"probYshared",110,0,1.1);
248 probYnoshare =
new TH1F(
"probYnoshare",
"probYnoshare",110,0,1.1);
258 LogDebug(
"TestTrackHits") <<
"new event" ;
264 iEvent.
getByLabel(
"offlineBeamSpot",beamSpot);
283 LogTrace(
"TestTrackHits") <<
"\n*****************new trajectory********************" ;
286 std::vector<TrajectoryMeasurement> tmColl = it->measurements();
297 std::vector<std::pair<TrackingParticleRef, double> > tP;
298 if(recSimColl.
find(track) != recSimColl.
end()){
299 tP = recSimColl[track];
301 edm::LogVerbatim(
"TestTrackHits") <<
"reco::Track #" << ++yyy <<
" with pt=" << track->
pt()
302 <<
" associated with quality:" << tP.begin()->second <<
" good track #" << ++yy <<
" has hits:" << track->
numberOfValidHits() <<
"\n";
305 edm::LogVerbatim(
"TestTrackHits") <<
"reco::Track #" << ++yyy <<
" with pt=" << track->
pt()
306 <<
" NOT associated to any TrackingParticle" <<
"\n";
321 LogTrace(
"TestTrackHits") <<
"a tp is associated with fraction=" << tP.begin()->second;
323 std::vector<unsigned int> tpids;
325 LogTrace(
"TestTrackHits") <<
"tp id=" << g4T->trackId();
326 tpids.push_back(g4T->trackId());
331 for (std::vector<TrajectoryMeasurement>::iterator tm=tmColl.begin();tm!=tmColl.end();++tm){
333 tchi2+=tm->estimate();
335 LogTrace(
"TestTrackHits") <<
"+++++++++++++++++new hit+++++++++++++++++" ;
336 CTTRHp rhit = tm->recHit();
339 TSOS state =
combiner(tm->backwardPredictedState(), tm->forwardPredictedState());
341 if (rhit->isValid()==0 && rhit->det()!=0)
continue;
345 int subdetId = rhit->det()->geographicalId().subdetId();
346 DetId id = rhit->det()->geographicalId();
347 int layerId = tTopo->layer(
id);
348 LogTrace(
"TestTrackHits") <<
"subdetId=" << subdetId <<
" layerId=" << layerId ;
350 const Surface * surf = rhit->surface();
351 if (surf==0)
continue;
353 double energyLoss_ = 0.;
354 unsigned int monoId = 0;
355 std::vector<double> energyLossM;
356 std::vector<double> energyLossS;
358 unsigned int simhitvecsize = assSimHits.size();
359 if (simhitvecsize==0)
continue;
361 std::vector<unsigned int> trackIds;
364 for(std::vector<PSimHit>::const_iterator
m=assSimHits.begin();
m<assSimHits.end();
m++){
365 unsigned int tId =
m->trackId();
366 if (
find(trackIds.begin(),trackIds.end(),tId)==trackIds.end()) trackIds.push_back(tId);
367 if (
m->energyLoss()>energyLoss_) {
369 energyLoss_ =
m->energyLoss();
371 if (dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())) {
372 if (monoId==0) monoId =
m->detUnitId();
373 if (monoId ==
m->detUnitId()){
374 energyLossM.push_back(
m->energyLoss());
377 energyLossS.push_back(
m->energyLoss());
381 energyLossM.push_back(
m->energyLoss());
398 double chi2increment = tm->estimate();
399 LogTrace(
"TestTrackHits") <<
"tm->estimate()=" << tm->estimate();
401 title <<
"Chi2Increment_" << subdetId <<
"-" << layerId;
404 title <<
"Chi2IncrementVsEta_" << subdetId <<
"-" << layerId;
410 bool mergedhit =
false;
411 if (dynamic_cast<const SiPixelRecHit*>(rhit->hit())){
412 clustersize = ((
const SiPixelRecHit*)(rhit->hit()))->cluster()->size() ;
416 if (simhitvecsize>1) mergedhit =
true;
418 else if (dynamic_cast<const SiStripRecHit2D*>(rhit->hit())){
419 clustersize = ((
const SiStripRecHit2D*)(rhit->hit()))->cluster()->amplitudes().size() ;
423 if (simhitvecsize>1) mergedhit =
true;
425 else if (dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())){
426 int clsize1 = ((
const SiStripMatchedRecHit2D*)(rhit->hit()))->monoCluster().amplitudes().size();
427 int clsize2 = ((
const SiStripMatchedRecHit2D*)(rhit->hit()))->stereoCluster().amplitudes().size();
428 if (clsize1>clsize2) clustersize = clsize1;
429 else clustersize = clsize2;
433 if (simhitvecsize>2) mergedhit =
true;
435 else if (dynamic_cast<const ProjectedSiStripRecHit2D*>(rhit->hit())){
440 if (simhitvecsize>1) mergedhit =
true;
452 bool goodhit =
false;
453 for(
size_t j=0;
j<simTrackIds.size();
j++){
454 LogTrace(
"TestTrackHits") <<
"hit id=" << simTrackIds[
j].first;
455 for (
size_t jj=0;
jj<tpids.size();
jj++){
456 if (simTrackIds[
j].
first == tpids[
jj]) goodhit =
true;
461 bool ioniOnly =
true;
464 if (energyLossM.size()>1&&energyLossS.size()<=1) {
465 sort(energyLossM.begin(),energyLossM.end(),greater<double>());
468 else if (energyLossS.size()>1&&energyLossM.size()<=1) {
469 sort(energyLossS.begin(),energyLossS.end(),greater<double>());
472 else if (energyLossS.size()>1&&energyLossM.size()>1) {
473 sort(energyLossM.begin(),energyLossM.end(),greater<double>());
474 sort(energyLossS.begin(),energyLossS.end(),greater<double>());
475 energyLossM[1]/energyLossM[0] > energyLossS[1]/energyLossS[0]
481 LogVerbatim(
"TestTrackHits") <<
"MERGED HIT" << std::endl;
482 unsigned int idc = 0;
483 for (
size_t jj=0;
jj<tpids.size();
jj++){
484 idc +=
std::count(trackIds.begin(),trackIds.end(),tpids[
jj]);
486 if (idc==trackIds.size()) {
489 for(std::vector<PSimHit>::const_iterator
m=assSimHits.begin()+1;
m<assSimHits.end();
m++){
490 if ((
m->processType()!=7&&
m->processType()!=8&&
m->processType()!=9)&&
abs(
m->particleType())!=11){
495 if (ioniOnly&&!shared) {
497 title <<
"Chi2DeltaHit_" << subdetId <<
"-" << layerId;
504 }
else if(!ioniOnly&&!shared) {
506 title <<
"Chi2NSharedHit_" << subdetId <<
"-" << layerId;
515 title <<
"Chi2SharedHit_" << subdetId <<
"-" << layerId;
524 for(std::vector<PSimHit>::const_iterator
m=assSimHits.begin();
m<assSimHits.end();
m++) {
525 unsigned int tId =
m->trackId();
526 LogVerbatim(
"TestTrackHits") <<
"component with id=" << tId <<
" eLoss=" <<
m->energyLoss() <<
" pType=" <<
m->processType();
527 if (
find(tpids.begin(),tpids.end(),tId)==tpids.end())
continue;
528 if (
m->processType()==2) {
532 LogVerbatim(
"TestTrackHits") << gpr <<
" " << gps <<
" " << ger;
533 ROOT::Math::SVector<double,3>
delta;
534 delta[0]=gpr.
x()-gps.
x();
535 delta[1]=gpr.
y()-gps.
y();
536 delta[2]=gpr.
z()-gps.
z();
537 LogVerbatim(
"TestTrackHits") << delta <<
" " << ger ;
538 double mpull =
sqrt(delta[0]*delta[0]/ger[0][0]+delta[1]*delta[1]/ger[1][1]+delta[2]*delta[2]/ger[2][2]);
539 LogVerbatim(
"TestTrackHits") <<
"hit pull=" << mpull;
547 title <<
"Chi2GoodHit_" << subdetId <<
"-" << layerId;
558 title <<
"Chi2BadHit_" << subdetId <<
"-" << layerId;
563 probXbad->Fill(pix->probabilityX());
564 probYbad->Fill(pix->probabilityY());
571 if (dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())) {
572 if (simhitvecsize>2 && goodhit) {
578 double rechitmatchedx = rhit->localPosition().x();
579 double rechitmatchedy = rhit->localPosition().y();
580 double mindist = 999999;
582 std::pair<LocalPoint,LocalVector> closestPair;
584 const BoundPlane& plane = (rhit)->det()->surface();
585 for(std::vector<PSimHit>::const_iterator
m=assSimHits.begin();
m<assSimHits.end();
m++){
587 std::pair<LocalPoint,LocalVector> hitPair =
projectHit((*
m),stripDet,plane);
588 distx = fabs(rechitmatchedx - hitPair.first.x());
589 disty = fabs(rechitmatchedy - hitPair.first.y());
590 double dist = distx*distx+disty*disty;
591 if(
sqrt(dist)<mindist){
593 closestPair = hitPair;
596 shitLPos = closestPair.first;
597 shitLMom = closestPair.second;
599 if (simhitvecsize>1 && goodhit) {
617 GlobalError rhitGPEr = (rhit)->globalPositionError();
619 LogVerbatim(
"TestTrackHits") <<
"assSimHits.size()=" << assSimHits.size() ;
620 LogVerbatim(
"TestTrackHits") <<
"tsos globalPos =" << tsosGPos ;
621 LogVerbatim(
"TestTrackHits") <<
"sim hit globalPos=" << shitGPos ;
622 LogVerbatim(
"TestTrackHits") <<
"rec hit globalPos=" << rhitGPos ;
623 LogVerbatim(
"TestTrackHits") <<
"geographicalId =" << rhit->det()->geographicalId().rawId() ;
627 if (rhit->detUnit())
LogTrace(
"TestTrackHits") <<
"rhit->detUnit()->geographicalId()="
628 << rhit->detUnit()->geographicalId().rawId() ;
629 LogTrace(
"TestTrackHits") <<
"rhit->det()->surface().position()="
630 << rhit->det()->surface().position() ;
631 if (rhit->detUnit())
LogTrace(
"TestTrackHits") <<
"rhit->detUnit()->surface().position()="
632 << rhit->detUnit()->surface().position() ;
633 LogTrace(
"TestTrackHits") <<
"rhit->det()->position()=" << rhit->det()->position() ;
634 if (rhit->detUnit())
LogTrace(
"TestTrackHits") <<
"rhit->detUnit()->position()="
635 << rhit->detUnit()->position() ;
636 LogTrace(
"TestTrackHits") <<
"rhit->det()->surface().bounds().length()="
637 << rhit->det()->surface().bounds().length() ;
638 if (rhit->detUnit())
LogTrace(
"TestTrackHits") <<
"rhit->detUnit()->surface().bounds().length()="
639 << rhit->detUnit()->surface().bounds().length() ;
640 LogTrace(
"TestTrackHits") <<
"rhit->det()->surface().bounds().width()="
641 << rhit->det()->surface().bounds().width() ;
642 if (rhit->detUnit())
LogTrace(
"TestTrackHits") <<
"rhit->detUnit()->surface().bounds().width()="
643 << rhit->detUnit()->surface().bounds().width() ;
644 LogTrace(
"TestTrackHits") <<
"rhit->det()->surface().bounds().thickness()="
645 << rhit->det()->surface().bounds().thickness() ;
646 if (rhit->detUnit())
LogTrace(
"TestTrackHits") <<
"rhit->detUnit()->surface().bounds().thickness()="
647 << rhit->detUnit()->surface().bounds().thickness() ;
650 double pullGPX_rs = (rhitGPos.
x()-shitGPos.
x())/
sqrt(rhitGPEr.
cxx());
651 double pullGPY_rs = (rhitGPos.
y()-shitGPos.
y())/
sqrt(rhitGPEr.
cyy());
652 double pullGPZ_rs = (rhitGPos.
z()-shitGPos.
z())/
sqrt(rhitGPEr.
czz());
660 title <<
"PullGP_X_" << subdetId <<
"-" << layerId <<
"_rs";
663 title <<
"PullGP_Y_" << subdetId <<
"-" << layerId <<
"_rs";
666 title <<
"PullGP_Z_" << subdetId <<
"-" << layerId <<
"_rs";
669 double pullGPX_tr = (tsosGPos.
x()-rhitGPos.
x())/
sqrt(tsosGPEr.
cxx()+rhitGPEr.
cxx());
670 double pullGPY_tr = (tsosGPos.
y()-rhitGPos.
y())/
sqrt(tsosGPEr.
cyy()+rhitGPEr.
cyy());
671 double pullGPZ_tr = (tsosGPos.
z()-rhitGPos.
z())/
sqrt(tsosGPEr.
czz()+rhitGPEr.
czz());
679 title <<
"PullGP_X_" << subdetId <<
"-" << layerId <<
"_tr";
682 title <<
"PullGP_Y_" << subdetId <<
"-" << layerId <<
"_tr";
685 title <<
"PullGP_Z_" << subdetId <<
"-" << layerId <<
"_tr";
688 double pullGPX_ts = (tsosGPos.
x()-shitGPos.
x())/
sqrt(tsosGPEr.
cxx());
689 double pullGPY_ts = (tsosGPos.
y()-shitGPos.
y())/
sqrt(tsosGPEr.
cyy());
690 double pullGPZ_ts = (tsosGPos.
z()-shitGPos.
z())/
sqrt(tsosGPEr.
czz());
695 LogTrace(
"TestTrackHits") <<
"ts1" ;
698 title <<
"PullGP_X_" << subdetId <<
"-" << layerId <<
"_ts";
701 title <<
"PullGP_Y_" << subdetId <<
"-" << layerId <<
"_ts";
704 title <<
"PullGP_Z_" << subdetId <<
"-" << layerId <<
"_ts";
707 double pullGMX_ts = (tsosGMom.
x()-shitGMom.
x())/
sqrt(tsosGMEr.
cxx());
708 double pullGMY_ts = (tsosGMom.
y()-shitGMom.
y())/
sqrt(tsosGMEr.
cyy());
709 double pullGMZ_ts = (tsosGMom.
z()-shitGMom.
z())/
sqrt(tsosGMEr.
czz());
714 LogTrace(
"TestTrackHits") <<
"ts2" ;
717 title <<
"PullGM_X_" << subdetId <<
"-" << layerId <<
"_ts";
720 title <<
"PullGM_Y_" << subdetId <<
"-" << layerId <<
"_ts";
723 title <<
"PullGM_Z_" << subdetId <<
"-" << layerId <<
"_ts";
725 if (dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())) {
728 LogTrace(
"TestTrackHits") <<
"MONO HIT" ;
729 auto m =
dynamic_cast<const SiStripMatchedRecHit2D*
>(rhit->hit())->monoHit();
731 if (tMonoHit==0)
continue;
733 if (assMonoSimHits.size()==0)
continue;
734 const PSimHit sMonoHit = *(assMonoSimHits.begin());
735 const Surface * monoSurf = &( tMonoHit->det()->surface() );
736 if (monoSurf==0)
continue;
737 TSOS monoState = thePropagatorAnyDir->
propagate(state,*monoSurf);
738 if (monoState.
isValid()==0)
continue;
753 GlobalPoint monoRhitGPos = tMonoHit->globalPosition();
754 GlobalError monoRhitGPEr = tMonoHit->globalPositionError();
756 double pullGPX_rs_mono = (monoRhitGPos.
x()-monoShitGPos.
x())/
sqrt(monoRhitGPEr.
cxx());
757 double pullGPY_rs_mono = (monoRhitGPos.
y()-monoShitGPos.
y())/
sqrt(monoRhitGPEr.
cyy());
758 double pullGPZ_rs_mono = (monoRhitGPos.
z()-monoShitGPos.
z())/
sqrt(monoRhitGPEr.
czz());
767 title <<
"Chi2Increment_mono_" << subdetId <<
"-" << layerId ;
771 title <<
"PullGP_X_" << subdetId <<
"-" << layerId <<
"_rs_mono";
774 title <<
"PullGP_Y_" << subdetId <<
"-" << layerId <<
"_rs_mono";
777 title <<
"PullGP_Z_" << subdetId <<
"-" << layerId <<
"_rs_mono";
780 double pullGPX_tr_mono = (monoTsosGPos.
x()-monoRhitGPos.
x())/
sqrt(monoTsosGPEr.
cxx()+monoRhitGPEr.
cxx());
781 double pullGPY_tr_mono = (monoTsosGPos.
y()-monoRhitGPos.
y())/
sqrt(monoTsosGPEr.
cyy()+monoRhitGPEr.
cyy());
782 double pullGPZ_tr_mono = (monoTsosGPos.
z()-monoRhitGPos.
z())/
sqrt(monoTsosGPEr.
czz()+monoRhitGPEr.
czz());
788 title <<
"PullGP_X_" << subdetId <<
"-" << layerId <<
"_tr_mono";
791 title <<
"PullGP_Y_" << subdetId <<
"-" << layerId <<
"_tr_mono";
794 title <<
"PullGP_Z_" << subdetId <<
"-" << layerId <<
"_tr_mono";
797 double pullGPX_ts_mono = (monoTsosGPos.
x()-monoShitGPos.
x())/
sqrt(monoTsosGPEr.
cxx());
798 double pullGPY_ts_mono = (monoTsosGPos.
y()-monoShitGPos.
y())/
sqrt(monoTsosGPEr.
cyy());
799 double pullGPZ_ts_mono = (monoTsosGPos.
z()-monoShitGPos.
z())/
sqrt(monoTsosGPEr.
czz());
805 title <<
"PullGP_X_" << subdetId <<
"-" << layerId <<
"_ts_mono";
808 title <<
"PullGP_Y_" << subdetId <<
"-" << layerId <<
"_ts_mono";
811 title <<
"PullGP_Z_" << subdetId <<
"-" << layerId <<
"_ts_mono";
814 double pullGMX_ts_mono = (monoTsosGMom.
x()-monoShitGMom.
x())/
sqrt(monoTsosGMEr.
cxx());
815 double pullGMY_ts_mono = (monoTsosGMom.
y()-monoShitGMom.
y())/
sqrt(monoTsosGMEr.
cyy());
816 double pullGMZ_ts_mono = (monoTsosGMom.
z()-monoShitGMom.
z())/
sqrt(monoTsosGMEr.
czz());
822 title <<
"PullGM_X_" << subdetId <<
"-" << layerId <<
"_ts_mono";
825 title <<
"PullGM_Y_" << subdetId <<
"-" << layerId <<
"_ts_mono";
828 title <<
"PullGM_Z_" << subdetId <<
"-" << layerId <<
"_ts_mono";
832 LogTrace(
"TestTrackHits") <<
"STEREO HIT" ;
833 auto s =
dynamic_cast<const SiStripMatchedRecHit2D*
>(rhit->hit())->stereoHit();
835 if (tStereoHit==0)
continue;
837 if (assStereoSimHits.size()==0)
continue;
838 const PSimHit sStereoHit = *(assStereoSimHits.begin());
839 const Surface * stereoSurf = &( tStereoHit->det()->surface() );
840 if (stereoSurf==0)
continue;
841 TSOS stereoState = thePropagatorAnyDir->
propagate(state,*stereoSurf);
842 if (stereoState.
isValid()==0)
continue;
857 GlobalPoint stereoRhitGPos = tStereoHit->globalPosition();
858 GlobalError stereoRhitGPEr = tStereoHit->globalPositionError();
864 title <<
"Chi2Increment_stereo_" << subdetId <<
"-" << layerId ;
867 double pullGPX_rs_stereo = (stereoRhitGPos.
x()-stereoShitGPos.
x())/
sqrt(stereoRhitGPEr.
cxx());
868 double pullGPY_rs_stereo = (stereoRhitGPos.
y()-stereoShitGPos.
y())/
sqrt(stereoRhitGPEr.
cyy());
869 double pullGPZ_rs_stereo = (stereoRhitGPos.
z()-stereoShitGPos.
z())/
sqrt(stereoRhitGPEr.
czz());
875 title <<
"PullGP_X_" << subdetId <<
"-" << layerId <<
"_rs_stereo";
878 title <<
"PullGP_Y_" << subdetId <<
"-" << layerId <<
"_rs_stereo";
881 title <<
"PullGP_Z_" << subdetId <<
"-" << layerId <<
"_rs_stereo";
884 double pullGPX_tr_stereo = (stereoTsosGPos.
x()-stereoRhitGPos.
x())/
sqrt(stereoTsosGPEr.
cxx()+stereoRhitGPEr.
cxx());
885 double pullGPY_tr_stereo = (stereoTsosGPos.
y()-stereoRhitGPos.
y())/
sqrt(stereoTsosGPEr.
cyy()+stereoRhitGPEr.
cyy());
886 double pullGPZ_tr_stereo = (stereoTsosGPos.
z()-stereoRhitGPos.
z())/
sqrt(stereoTsosGPEr.
czz()+stereoRhitGPEr.
czz());
892 title <<
"PullGP_X_" << subdetId <<
"-" << layerId <<
"_tr_stereo";
895 title <<
"PullGP_Y_" << subdetId <<
"-" << layerId <<
"_tr_stereo";
898 title <<
"PullGP_Z_" << subdetId <<
"-" << layerId <<
"_tr_stereo";
901 double pullGPX_ts_stereo = (stereoTsosGPos.
x()-stereoShitGPos.
x())/
sqrt(stereoTsosGPEr.
cxx());
902 double pullGPY_ts_stereo = (stereoTsosGPos.
y()-stereoShitGPos.
y())/
sqrt(stereoTsosGPEr.
cyy());
903 double pullGPZ_ts_stereo = (stereoTsosGPos.
z()-stereoShitGPos.
z())/
sqrt(stereoTsosGPEr.
czz());
909 title <<
"PullGP_X_" << subdetId <<
"-" << layerId <<
"_ts_stereo";
912 title <<
"PullGP_Y_" << subdetId <<
"-" << layerId <<
"_ts_stereo";
915 title <<
"PullGP_Z_" << subdetId <<
"-" << layerId <<
"_ts_stereo";
918 double pullGMX_ts_stereo = (stereoTsosGMom.
x()-stereoShitGMom.
x())/
sqrt(stereoTsosGMEr.
cxx());
919 double pullGMY_ts_stereo = (stereoTsosGMom.
y()-stereoShitGMom.
y())/
sqrt(stereoTsosGMEr.
cyy());
920 double pullGMZ_ts_stereo = (stereoTsosGMom.
z()-stereoShitGMom.
z())/
sqrt(stereoTsosGMEr.
czz());
926 title <<
"PullGM_X_" << subdetId <<
"-" << layerId <<
"_ts_stereo";
929 title <<
"PullGM_Y_" << subdetId <<
"-" << layerId <<
"_ts_stereo";
932 title <<
"PullGM_Z_" << subdetId <<
"-" << layerId <<
"_ts_stereo";
936 LogTrace(
"TestTrackHits") <<
"traj chi2=" << tchi2 ;
937 LogTrace(
"TestTrackHits") <<
"track chi2=" << track->
chi2() ;
940 LogTrace(
"TestTrackHits") <<
"end of event: processd hits=" << evtHits ;
946 TDirectory * chi2d =
file->mkdir(
"Chi2Increment");
948 TDirectory * gp_ts =
file->mkdir(
"GP_TSOS-SimHit");
949 TDirectory * gm_ts =
file->mkdir(
"GM_TSOS-SimHit");
950 TDirectory * gp_tr =
file->mkdir(
"GP_TSOS-RecHit");
951 TDirectory * gp_rs =
file->mkdir(
"GP_RecHit-SimHit");
953 TDirectory * gp_tsx = gp_ts->mkdir(
"X");
954 TDirectory * gp_tsy = gp_ts->mkdir(
"Y");
955 TDirectory * gp_tsz = gp_ts->mkdir(
"Z");
956 TDirectory * gm_tsx = gm_ts->mkdir(
"X");
957 TDirectory * gm_tsy = gm_ts->mkdir(
"Y");
958 TDirectory * gm_tsz = gm_ts->mkdir(
"Z");
959 TDirectory * gp_trx = gp_tr->mkdir(
"X");
960 TDirectory * gp_try = gp_tr->mkdir(
"Y");
961 TDirectory * gp_trz = gp_tr->mkdir(
"Z");
962 TDirectory * gp_rsx = gp_rs->mkdir(
"X");
963 TDirectory * gp_rsy = gp_rs->mkdir(
"Y");
964 TDirectory * gp_rsz = gp_rs->mkdir(
"Z");
966 TDirectory * gp_tsx_mono = gp_ts->mkdir(
"MONOX");
967 TDirectory * gp_tsy_mono = gp_ts->mkdir(
"MONOY");
968 TDirectory * gp_tsz_mono = gp_ts->mkdir(
"MONOZ");
969 TDirectory * gm_tsx_mono = gm_ts->mkdir(
"MONOX");
970 TDirectory * gm_tsy_mono = gm_ts->mkdir(
"MONOY");
971 TDirectory * gm_tsz_mono = gm_ts->mkdir(
"MONOZ");
972 TDirectory * gp_trx_mono = gp_tr->mkdir(
"MONOX");
973 TDirectory * gp_try_mono = gp_tr->mkdir(
"MONOY");
974 TDirectory * gp_trz_mono = gp_tr->mkdir(
"MONOZ");
975 TDirectory * gp_rsx_mono = gp_rs->mkdir(
"MONOX");
976 TDirectory * gp_rsy_mono = gp_rs->mkdir(
"MONOY");
977 TDirectory * gp_rsz_mono = gp_rs->mkdir(
"MONOZ");
979 TDirectory * gp_tsx_stereo = gp_ts->mkdir(
"STEREOX");
980 TDirectory * gp_tsy_stereo = gp_ts->mkdir(
"STEREOY");
981 TDirectory * gp_tsz_stereo = gp_ts->mkdir(
"STEREOZ");
982 TDirectory * gm_tsx_stereo = gm_ts->mkdir(
"STEREOX");
983 TDirectory * gm_tsy_stereo = gm_ts->mkdir(
"STEREOY");
984 TDirectory * gm_tsz_stereo = gm_ts->mkdir(
"STEREOZ");
985 TDirectory * gp_trx_stereo = gp_tr->mkdir(
"STEREOX");
986 TDirectory * gp_try_stereo = gp_tr->mkdir(
"STEREOY");
987 TDirectory * gp_trz_stereo = gp_tr->mkdir(
"STEREOZ");
988 TDirectory * gp_rsx_stereo = gp_rs->mkdir(
"STEREOX");
989 TDirectory * gp_rsy_stereo = gp_rs->mkdir(
"STEREOY");
990 TDirectory * gp_rsz_stereo = gp_rs->mkdir(
"STEREOZ");
1029 for (
int i=0;
i!=6;
i++)
1030 for (
int j=0;
j!=9;
j++){
1031 if (
i==0 &&
j>2)
break;
1032 if (
i==1 &&
j>1)
break;
1033 if (
i==2 &&
j>3)
break;
1034 if (
i==3 &&
j>2)
break;
1035 if (
i==4 &&
j>5)
break;
1036 if (
i==5 &&
j>8)
break;
1039 title <<
"Chi2Increment_" <<
i+1 <<
"-" <<
j+1 ;
1042 title <<
"Chi2IncrementVsEta_" <<
i+1 <<
"-" << j+1 ;
1045 title <<
"Chi2GoodHit_" <<
i+1 <<
"-" << j+1 ;
1048 title <<
"Chi2BadHit_" <<
i+1 <<
"-" << j+1 ;
1051 title <<
"Chi2DeltaHit_" <<
i+1 <<
"-" << j+1 ;
1054 title <<
"Chi2NSharedHit_" <<
i+1 <<
"-" << j+1 ;
1057 title <<
"Chi2SharedHit_" <<
i+1 <<
"-" << j+1 ;
1063 title <<
"PullGP_X_" <<
i+1 <<
"-" << j+1 <<
"_ts";
1067 title <<
"PullGP_Y_" <<
i+1 <<
"-" << j+1 <<
"_ts";
1071 title <<
"PullGP_Z_" <<
i+1 <<
"-" << j+1 <<
"_ts";
1077 title <<
"PullGM_X_" <<
i+1 <<
"-" << j+1 <<
"_ts";
1081 title <<
"PullGM_Y_" <<
i+1 <<
"-" << j+1 <<
"_ts";
1085 title <<
"PullGM_Z_" <<
i+1 <<
"-" << j+1 <<
"_ts";
1091 title <<
"PullGP_X_" <<
i+1 <<
"-" << j+1 <<
"_tr";
1095 title <<
"PullGP_Y_" <<
i+1 <<
"-" << j+1 <<
"_tr";
1099 title <<
"PullGP_Z_" <<
i+1 <<
"-" << j+1 <<
"_tr";
1105 title <<
"PullGP_X_" <<
i+1 <<
"-" << j+1 <<
"_rs";
1109 title <<
"PullGP_Y_" <<
i+1 <<
"-" << j+1 <<
"_rs";
1113 title <<
"PullGP_Z_" <<
i+1 <<
"-" << j+1 <<
"_rs";
1116 if ( ((
i==2||
i==4)&&(j==0||j==1)) || (
i==3||
i==5) ){
1119 title <<
"Chi2Increment_mono_" <<
i+1 <<
"-" << j+1 ;
1122 title <<
"Chi2Increment_stereo_" <<
i+1 <<
"-" << j+1 ;
1128 title <<
"PullGP_X_" <<
i+1 <<
"-" << j+1 <<
"_ts_mono";
1132 title <<
"PullGP_Y_" <<
i+1 <<
"-" << j+1 <<
"_ts_mono";
1136 title <<
"PullGP_Z_" <<
i+1 <<
"-" << j+1 <<
"_ts_mono";
1142 title <<
"PullGM_X_" <<
i+1 <<
"-" << j+1 <<
"_ts_mono";
1146 title <<
"PullGM_Y_" <<
i+1 <<
"-" << j+1 <<
"_ts_mono";
1150 title <<
"PullGM_Z_" <<
i+1 <<
"-" << j+1 <<
"_ts_mono";
1156 title <<
"PullGP_X_" <<
i+1 <<
"-" << j+1 <<
"_tr_mono";
1160 title <<
"PullGP_Y_" <<
i+1 <<
"-" << j+1 <<
"_tr_mono";
1164 title <<
"PullGP_Z_" <<
i+1 <<
"-" << j+1 <<
"_tr_mono";
1170 title <<
"PullGP_X_" <<
i+1 <<
"-" << j+1 <<
"_rs_mono";
1174 title <<
"PullGP_Y_" <<
i+1 <<
"-" << j+1 <<
"_rs_mono";
1178 title <<
"PullGP_Z_" <<
i+1 <<
"-" << j+1 <<
"_rs_mono";
1183 gp_tsx_stereo->cd();
1185 title <<
"PullGP_X_" <<
i+1 <<
"-" << j+1 <<
"_ts_stereo";
1187 gp_tsy_stereo->cd();
1189 title <<
"PullGP_Y_" <<
i+1 <<
"-" << j+1 <<
"_ts_stereo";
1191 gp_tsz_stereo->cd();
1193 title <<
"PullGP_Z_" <<
i+1 <<
"-" << j+1 <<
"_ts_stereo";
1197 gm_tsx_stereo->cd();
1199 title <<
"PullGM_X_" <<
i+1 <<
"-" << j+1 <<
"_ts_stereo";
1201 gm_tsy_stereo->cd();
1203 title <<
"PullGM_Y_" <<
i+1 <<
"-" << j+1 <<
"_ts_stereo";
1205 gm_tsz_stereo->cd();
1207 title <<
"PullGM_Z_" <<
i+1 <<
"-" << j+1 <<
"_ts_stereo";
1211 gp_trx_stereo->cd();
1213 title <<
"PullGP_X_" <<
i+1 <<
"-" << j+1 <<
"_tr_stereo";
1215 gp_try_stereo->cd();
1217 title <<
"PullGP_Y_" <<
i+1 <<
"-" << j+1 <<
"_tr_stereo";
1219 gp_trz_stereo->cd();
1221 title <<
"PullGP_Z_" <<
i+1 <<
"-" << j+1 <<
"_tr_stereo";
1225 gp_rsx_stereo->cd();
1227 title <<
"PullGP_X_" <<
i+1 <<
"-" << j+1 <<
"_rs_stereo";
1229 gp_rsy_stereo->cd();
1231 title <<
"PullGP_Y_" <<
i+1 <<
"-" << j+1 <<
"_rs_stereo";
1233 gp_rsz_stereo->cd();
1235 title <<
"PullGP_Z_" <<
i+1 <<
"-" << j+1 <<
"_rs_stereo";
1246 std::pair<LocalPoint,LocalVector>
1251 LocalPoint localHit = plane.toLocal(globalpos);
1258 float scale = -localHit.
z() / dir.
z();
1268 return std::pair<LocalPoint,LocalVector>( projectedPos, localStripDir);
1271 template<
unsigned int D>
1278 SMatDD
R = asSMatrix<D>(rhit->parametersError()) + me.
measuredError<
D>(*rhit);
1280 return ROOT::Math::Similarity(r,R) ;
GlobalPoint toGlobal(const Point2DBase< Scalar, LocalTag > lp) const
T getParameter(std::string const &) const
std::map< std::string, TH1F * > hPullGP_X_ts_mono
TH2F * hPrjClsize_vs_Chi2
std::map< std::string, TH1F * > hPullGM_X_ts_mono
std::map< std::string, TH1F * > hPullGM_Y_ts_stereo
ROOT::Math::SMatrix< double, D1, D1, ROOT::Math::MatRepSym< double, D1 > > SymMatrix
std::map< std::string, TH1F * > hPullGP_Y_ts_mono
virtual float stripAngle(float strip) const =0
const_iterator end() const
last iterator over the map (read only)
LocalVector momentumAtEntry() const
The momentum of the track that produced the hit, at entry point.
std::map< std::string, TH1F * > hPullGP_Z_rs_mono
std::map< std::string, TH1F * > hPullGP_Z_rs_stereo
#define DEFINE_FWK_MODULE(type)
TH2F * hSt2Clsize_vs_Chi2
TrackerHitAssociator * hitAssociator
std::map< std::string, TH1F * > hPullGP_Y_rs_stereo
Sin< T >::type sin(const T &t)
std::map< std::string, TH1F * > hPullGP_Z_ts_mono
const_iterator find(const key_type &k) const
find element with specified reference key
std::map< std::string, TH1F * > hPullGP_Y_rs
std::map< std::string, TH1F * > hPullGP_Z_tr_mono
edm::Handle< std::vector< Trajectory > > trajCollectionHandle
const CartesianTrajectoryError cartesianError() const
const edm::ParameterSet conf_
std::map< std::string, TH1F * > hChi2Increment
std::map< std::string, TH1F * > hChi2BadHit
TH2F * hSt1Clsize_vs_Chi2
GlobalPoint globalPosition() const
std::map< std::string, TH1F * > hChi2SharedHit
std::map< std::string, TH1F * > hChi2Increment_mono
std::map< std::string, TH1F * > hPullGP_Y_ts
std::map< std::string, TH1F * > hPullGP_X_rs
std::map< std::string, TH1F * > hPullGP_X_tr_stereo
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
virtual const StripTopology & specificTopology() const
Returns a reference to the strip proxy topology.
std::vector< SimHitIdpr > associateHitId(const TrackingRecHit &thit)
const Plane & surface() const
The nominal surface of the GeomDet.
std::map< std::string, TH1F * > hPullGP_Y_tr
std::map< std::string, TH1F * > hChi2DeltaHit
TH2F * hPixClsize_vs_Chi2
virtual float strip(const LocalPoint &) const =0
edm::ESHandle< Propagator > thePropagator
edm::Handle< edm::View< reco::Track > > trackCollectionHandle
std::map< std::string, TH1F * > hPullGM_Z_ts_stereo
double eta() const
pseudorapidity of momentum vector
std::map< std::string, TH1F * > hPullGP_Z_ts
ROOT::Math::SMatrix< double, 3, 3, ROOT::Math::MatRepSym< double, 3 > > AlgebraicSymMatrix33
edm::ESHandle< MagneticField > theMF
Local3DPoint localPosition() const
double chi2() const
chi-squared of the fit
std::map< std::string, TH1F * > hPullGP_X_ts
std::map< std::string, TH1F * > hPullGP_Z_rs
TestTrackHits(const edm::ParameterSet &)
edm::ESHandle< TrajectoryStateUpdator > theUpdator
std::map< std::string, TH1F * > hChi2Increment_stereo
double pt() const
track transverse momentum
std::map< std::string, TH1F * > hChi2GoodHit
std::map< std::string, TH1F * > hPullGM_Z_ts
Cos< T >::type cos(const T &t)
std::map< std::string, TH1F * > hPullGP_Z_tr
Abs< T >::type abs(const T &t)
TransientTrackingRecHit::ConstRecHitPointer CTTRHp
std::map< std::string, TH1F * > hPullGP_X_ts_stereo
unsigned short numberOfValidHits() const
number of valid hits found
virtual void beginRun(edm::Run &run, const edm::EventSetup &)
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
const AlgebraicSymMatrix66 & matrix() const
virtual TrajectoryStateOnSurface propagate(const FreeTrajectoryState &, const Surface &) const
edm::ESHandle< TransientTrackingRecHitBuilder > theBuilder
LocalVector localDirection() const
Obsolete. Same as momentumAtEntry().unit(), for backward compatibility.
std::vector< SimTrack >::const_iterator g4t_iterator
edm::Handle< TrackingParticleCollection > trackingParticleCollectionHandle
std::map< std::string, TH1F * > hPullGP_Y_rs_mono
TrajectoryStateOnSurface TSOS
ROOT::Math::SVector< double, D1 > Vector
TH1F * hTotChi2NSharedHit
std::map< std::string, TH1F * > hPullGM_Y_ts_mono
const GlobalError position() const
Position error submatrix.
std::map< std::string, TH1F * > hPullGP_X_rs_mono
virtual void analyze(const edm::Event &, const edm::EventSetup &)
std::map< std::string, TH1F * > hPullGP_Y_tr_stereo
std::map< std::string, TH1F * > hPullGP_Z_tr_stereo
edm::ESHandle< TrackAssociatorBase > trackAssociator
T const * product() const
std::map< std::string, TH1F * > hPullGP_X_tr
std::map< std::string, TH1F * > hPullGP_X_tr_mono
unsigned short processType() const
edm::ESHandle< TrackerGeometry > theG
std::map< std::string, TH1F * > hChi2NSharedHit
GlobalVector globalMomentum() const
std::pair< LocalPoint, LocalVector > projectHit(const PSimHit &, const StripGeomDetUnit *, const BoundPlane &)
std::map< std::string, TH1F * > hPullGM_X_ts
std::vector< PSimHit > associateHit(const TrackingRecHit &thit)
std::map< std::string, TH1F * > hPullGP_Z_ts_stereo
std::map< std::string, TH1F * > hPullGP_X_rs_stereo
std::map< std::string, TH1F * > hPullGP_Y_ts_stereo
edm::Handle< TrajTrackAssociationCollection > trajTrackAssociationCollectionHandle
std::map< std::string, TH1F * > hPullGM_Z_ts_mono
double computeChi2Increment(MeasurementExtractor, TransientTrackingRecHit::ConstRecHitPointer)
std::map< std::string, TH2F * > hChi2IncrementVsEta
const PositionType & position() const
std::string propagatorName
DecomposeProduct< arg, typename Div::arg > D
std::map< std::string, TH1F * > hPullGM_X_ts_stereo
std::map< std::string, TH1F * > hPullGM_Y_ts
std::map< std::string, TH1F * > hPullGP_Y_tr_mono