Method called before the event loop.
233 using namespace reco;
235 edm::LogInfo(
"MuonTrackValidator") <<
"\n====================================================" <<
"\n"
236 <<
"Analyzing new event" <<
"\n"
237 <<
"====================================================\n" <<
"\n";
240 std::unique_ptr<ParametersDefinerForTP> Cosmic_parametersDefinerTP;
250 Cosmic_parametersDefinerTP = _Cosmic_parametersDefinerTP->clone();
255 Cosmic_parametersDefinerTP->initEvent(simHitsTPAssoc);
274 std::vector<const reco::TrackToTrackingParticleAssociator*> associator;
279 associator.push_back( theAssociator.
product() );
284 for (
unsigned int ww=0;ww<
associators.size();ww++){
285 for (
unsigned int www=0;www<
label.size();www++){
293 unsigned int trackCollectionSize = 0;
304 trackCollectionSize = trackCollection->size();
309 <<
label[www].process()<<
":"
310 <<
label[www].label()<<
":"
311 <<
label[www].instance()<<
" with "
314 LogTrace(
"MuonTrackValidator") <<
"Calling associateRecoToSim method" <<
"\n";
315 recSimColl=associator[ww]->associateRecoToSim(trackCollection,
317 LogTrace(
"MuonTrackValidator") <<
"Calling associateSimToReco method" <<
"\n";
318 simRecColl=associator[ww]->associateSimToReco(trackCollection,
323 <<
label[www].process()<<
":"
324 <<
label[www].label()<<
":"
325 <<
label[www].instance()<<
" with "
332 simRecColl= *(simtorecoCollectionH.
product());
336 recSimColl= *(recotosimCollectionH.
product());
345 edm::LogVerbatim(
"MuonTrackValidator") <<
"\n# of TrackingParticles: " << tPCeff.size() <<
"\n";
349 bool TP_is_matched =
false;
351 bool Quality05 =
false;
352 bool Quality075 =
false;
366 LogTrace(
"MuonTrackValidator") <<
"TrackingParticle "<<
i;
372 dxySim = (-vertex.x()*
sin(momentum.phi())+vertex.y()*
cos(momentum.phi()));
373 dzSim = vertex.z() - (vertex.x()*momentum.x()+vertex.y()*momentum.y()) /
374 sqrt(momentum.perp2()) * momentum.z()/
sqrt(momentum.perp2());
381 momentumTP = Cosmic_parametersDefinerTP->momentum(event,setup,tpr);
382 vertexTP = Cosmic_parametersDefinerTP->vertex(event,setup,tpr);
383 dxySim = (-vertexTP.x()*
sin(momentumTP.phi())+vertexTP.y()*
cos(momentumTP.phi()));
384 dzSim = vertexTP.z() - (vertexTP.x()*momentumTP.x()+vertexTP.y()*momentumTP.y()) /
385 sqrt(momentumTP.perp2()) * momentumTP.z()/
sqrt(momentumTP.perp2());
387 edm::LogVerbatim(
"MuonTrackValidator") <<
"--------------------Selected TrackingParticle #"<<tpr.key();
388 edm::LogVerbatim(
"MuonTrackValidator") <<
"momentumTP: pt = "<<
sqrt(momentumTP.perp2())<<
", pz = "<<momentumTP.z()
389 <<
", \t vertexTP: radius = "<<
sqrt(vertexTP.perp2())<<
", z = "<<vertexTP.z() <<
"\n";
396 std::vector<std::pair<RefToBase<Track>,
double> > rt;
397 if(simRecColl.
find(tpr) != simRecColl.
end()){
398 rt = (std::vector<std::pair<RefToBase<Track>,
double> >) simRecColl[tpr];
401 edm::LogVerbatim(
"MuonTrackValidator")<<
"-----------------------------associated Track #"<<assoc_recoTrack.
key();
402 TP_is_matched =
true;
404 quality = rt.begin()->second;
406 <<
" with pt=" <<
sqrt(momentumTP.perp2())
407 <<
" associated with quality:" << quality <<
"\n";
409 if (quality > 0.75) {
413 else if (quality > 0.5) {
420 <<
"TrackingParticle #" << tpr.key()
421 <<
" with pt,eta,phi: "
422 <<
sqrt(momentumTP.perp2()) <<
" , "
423 << momentumTP.eta() <<
" , "
424 << momentumTP.phi() <<
" , "
425 <<
" NOT associated to any reco::Track" <<
"\n";
440 else if (Quality05) {
460 else if (Quality05) {
480 else if (Quality05) {
530 nSimHits= tpr.get()->numberOfHits();
533 nSimHits= tpr.get()->numberOfHits() - tpr.get()->numberOfTrackerHits();
536 nSimHits=tpr.get()->numberOfTrackerHits();
540 edm::LogVerbatim(
"MuonTrackValidator") <<
"\t N simhits = "<< nSimHits<<
"\n";
557 <<
label[www].process()<<
":"
558 <<
label[www].label()<<
":"
559 <<
label[www].instance()
560 <<
": " << trackCollectionSize <<
"\n";
564 bool Track_is_matched =
false;
568 std::vector<std::pair<TrackingParticleRef, double> > tp;
573 edm::LogVerbatim(
"MuonTrackValidator")<<
"----------------------------------------Track #"<< track.key();
575 if(recSimColl.
find(track) != recSimColl.
end()) {
576 tp = recSimColl[track];
577 if (tp.size() != 0) {
578 tpr = tp.begin()->first;
580 if(simRecColl.
find(tpr) != simRecColl.
end()) {
581 std::vector<std::pair<RefToBase<Track>,
double> > track_checkback = simRecColl[tpr];
583 assoc_track_checkback = track_checkback.begin()->first;
585 if ( assoc_track_checkback.
key() == track.key() ) {
586 edm::LogVerbatim(
"MuonTrackValidator")<<
"------------------associated TrackingParticle #"<<tpr.
key();
587 Track_is_matched =
true;
589 double Purity = tp.begin()->second;
590 double Quality = track_checkback.begin()->second;
591 edm::LogVerbatim(
"MuonTrackValidator") <<
"reco::Track #" << track.key() <<
" with pt=" << track->pt()
592 <<
" associated with quality:" << Purity <<
"\n";
599 if (!Track_is_matched)
601 <<
"reco::Track #" << track.key() <<
" with pt=" << track->pt() <<
" NOT associated to any TrackingParticle" <<
"\n";
605 if(recSimColl.
find(track) != recSimColl.
end()){
606 tp = recSimColl[track];
608 Track_is_matched =
true;
609 tpr = tp.begin()->first;
611 edm::LogVerbatim(
"MuonTrackValidator") <<
"reco::Track #" << track.key() <<
" with pt=" << track->pt()
612 <<
" associated with quality:" << tp.begin()->second <<
"\n";
615 edm::LogVerbatim(
"MuonTrackValidator") <<
"reco::Track #" << track.key() <<
" with pt=" << track->pt()
616 <<
" NOT associated to any TrackingParticle" <<
"\n";
625 if (Track_is_matched) {
635 if (Track_is_matched) {
645 if (Track_is_matched) {
655 if (Track_is_matched) {
665 if (Track_is_matched) {
675 edm::LogVerbatim(
"MuonTrackValidator") <<
"\t N valid rechits = "<< (int)track->found() <<
"\n";
685 if (! (Track_is_matched &&
tpSelector(*tpp)) )
continue;
687 momentumTP = Lhc_parametersDefinerTP->
momentum(event,setup,tpr) ;
688 vertexTP = Lhc_parametersDefinerTP->vertex(event,setup,tpr);
693 if (! (Track_is_matched &&
cosmictpSelector(tpr,&bs,event,setup)) )
continue;
695 momentumTP = Cosmic_parametersDefinerTP->momentum(event,setup,tpr) ;
696 vertexTP = Cosmic_parametersDefinerTP->vertex(event,setup,tpr);
702 double assocChi2 = -tp.begin()->second;
706 else if (
associators[ww]==
"trackAssociatorByHits"){
707 double fraction = tp.begin()->second;
713 h_nchi2[
w]->Fill(track->normalizedChi2());
714 h_nchi2_prob[
w]->Fill(TMath::Prob(track->chi2(),(int)track->ndof()));
715 h_hits[
w]->Fill(track->numberOfValidHits());
717 chi2_vs_nhits[
w]->Fill(track->numberOfValidHits(),track->normalizedChi2());
720 double ptSim =
sqrt(momentumTP.perp2());
721 double qoverpSim = tpr->charge()/
sqrt(momentumTP.x()*momentumTP.x()+momentumTP.y()*momentumTP.y()+momentumTP.z()*momentumTP.z());
722 double thetaSim = momentumTP.theta();
723 double lambdaSim =
M_PI/2-momentumTP.theta();
724 double phiSim = momentumTP.phi();
725 double dxySim = (-vertexTP.x()*
sin(momentumTP.phi())+vertexTP.y()*
cos(momentumTP.phi()));
726 double dzSim = vertexTP.z() - (vertexTP.x()*momentumTP.x()+vertexTP.y()*momentumTP.y()) /
727 sqrt(momentumTP.perp2()) * momentumTP.z()/
sqrt(momentumTP.perp2());
733 double qoverpErrorRec(0);
735 double ptErrorRec(0);
737 double lambdaErrorRec(0);
739 double phiErrorRec(0);
744 gsfTrack =
dynamic_cast<const GsfTrack*
>(&(*track));
745 if (gsfTrack==0)
edm::LogInfo(
"MuonTrackValidator") <<
"Trying to access mode for a non-GsfTrack";
750 getRecoMomentum(*gsfTrack, ptRec, ptErrorRec, qoverpRec, qoverpErrorRec,
751 lambdaRec,lambdaErrorRec, phiRec, phiErrorRec);
757 lambdaRec,lambdaErrorRec, phiRec, phiErrorRec);
760 double thetaRec = track->theta();
761 double ptError = ptErrorRec;
762 double ptres = ptRec - ptSim;
763 double etares = track->eta()-momentumTP.Eta();
764 double dxyRec = track->dxy(bs.position());
765 double dzRec = track->dz(bs.position());
767 double qoverpPull=(qoverpRec-qoverpSim)/qoverpErrorRec;
768 double thetaPull=(lambdaRec-lambdaSim)/lambdaErrorRec;
769 double phiDiff = phiRec - phiSim;
771 if (phiDiff >0.) phiDiff = phiDiff - 2.*
M_PI;
772 else phiDiff = phiDiff + 2.*
M_PI;
774 double phiPull=phiDiff/phiErrorRec;
775 double dxyPull=(dxyRec-dxySim)/track->dxyError();
776 double dzPull=(dzRec-dzSim)/track->dzError();
778 double contrib_Qoverp = ((qoverpRec-qoverpSim)/qoverpErrorRec)*((qoverpRec-qoverpSim)/qoverpErrorRec)/5;
779 double contrib_dxy = ((dxyRec-dxySim)/track->dxyError())*((dxyRec-dxySim)/track->dxyError())/5;
780 double contrib_dz = ((dzRec-dzSim)/track->dzError())*((dzRec-dzSim)/track->dzError())/5;
781 double contrib_theta = ((lambdaRec-lambdaSim)/lambdaErrorRec)*((lambdaRec-lambdaSim)/lambdaErrorRec)/5;
782 double contrib_phi = (phiDiff/phiErrorRec)*(phiDiff/phiErrorRec)/5;
784 edm::LogVerbatim(
"MuonTrackValidator") <<
"assocChi2=" << tp.begin()->second <<
"\n"
786 <<
"ptREC=" << ptRec <<
"\n"
787 <<
"etaREC=" << track->eta() <<
"\n"
788 <<
"qoverpREC=" << qoverpRec <<
"\n"
789 <<
"dxyREC=" << dxyRec <<
"\n"
790 <<
"dzREC=" << dzRec <<
"\n"
791 <<
"thetaREC=" << track->theta() <<
"\n"
792 <<
"phiREC=" << phiRec <<
"\n"
794 <<
"qoverpError()=" << qoverpErrorRec <<
"\n"
795 <<
"dxyError()=" << track->dxyError() <<
"\n"
796 <<
"dzError()=" << track->dzError() <<
"\n"
797 <<
"thetaError()=" << lambdaErrorRec <<
"\n"
798 <<
"phiError()=" << phiErrorRec <<
"\n"
800 <<
"ptSIM=" << ptSim <<
"\n"
801 <<
"etaSIM=" << momentumTP.Eta() <<
"\n"
802 <<
"qoverpSIM=" << qoverpSim <<
"\n"
803 <<
"dxySIM=" << dxySim <<
"\n"
804 <<
"dzSIM=" << dzSim <<
"\n"
805 <<
"thetaSIM=" <<
M_PI/2-lambdaSim <<
"\n"
806 <<
"phiSIM=" << phiSim <<
"\n"
808 <<
"contrib_Qoverp=" << contrib_Qoverp <<
"\n"
809 <<
"contrib_dxy=" << contrib_dxy <<
"\n"
810 <<
"contrib_dz=" << contrib_dz <<
"\n"
811 <<
"contrib_theta=" << contrib_theta <<
"\n"
812 <<
"contrib_phi=" << contrib_phi <<
"\n"
814 <<
"chi2PULL="<<contrib_Qoverp+contrib_dxy+contrib_dz+contrib_theta+contrib_phi<<
"\n";
822 h_pt[
w]->Fill(ptres/ptError);
870 nSimHits= tpr.
get()->numberOfHits();
873 nSimHits= tpr.
get()->numberOfHits() - tpr.
get()->numberOfTrackerHits();
876 nSimHits=tpr.
get()->numberOfTrackerHits();
886 <<
"Total Associated (simToReco): " << ats <<
"\n"
887 <<
"Total Reconstructed: " << rT <<
"\n"
888 <<
"Total Associated (recoToSim): " << at <<
"\n"
889 <<
"Total Fakes: " << rT-at <<
"\n";
std::vector< MonitorElement * > nrecHit_vs_nsimHit_sim2rec
std::vector< std::vector< int > > totASS_zpos
std::vector< MonitorElement * > h_PurityVsQuality
std::vector< std::vector< int > > totASSeta_Quality05
std::vector< std::vector< int > > totREC_phi
std::vector< MonitorElement * > etares_vs_eta
std::vector< MonitorElement * > nrecHit_vs_nsimHit_rec2sim
std::vector< MonitorElement * > ptres_vs_phi
std::vector< std::vector< int > > totASS2_hit
std::vector< edm::EDGetTokenT< edm::View< reco::Track > > > track_Collection_Token
std::vector< std::vector< int > > totRECpT
std::vector< TrackingParticle > TrackingParticleCollection
std::vector< std::vector< int > > totREC_dz
Vector momentum() const
spatial momentum vector
std::vector< MonitorElement * > nrec_vs_nsim
const_iterator end() const
last iterator over the map (read only)
std::vector< MonitorElement * > h_pullDz
std::vector< MonitorElement * > h_etaSIM
bool ignoremissingtkcollection_
bool getByToken(EDGetToken token, Handle< PROD > &result) const
std::vector< MonitorElement * > h_fakes
std::vector< MonitorElement * > ptres_vs_eta
std::vector< MonitorElement * > cotThetares_vs_eta
virtual double getEta(double eta)
std::vector< MonitorElement * > h_eta
std::vector< MonitorElement * > phimean_vs_eta_phi
Sin< T >::type sin(const T &t)
std::vector< std::vector< int > > totASS_dz
const_iterator find(const key_type &k) const
find element with specified reference key
std::vector< std::vector< double > > etaintervals
std::vector< MonitorElement * > h_pt
std::vector< MonitorElement * > phires_vs_pt
std::vector< MonitorElement * > h_pullQoverp
std::string parametersDefiner
std::vector< MonitorElement * > h_assochi2_prob
std::vector< std::vector< int > > totSIM_zpos
std::vector< std::vector< int > > totASS2_dz
std::vector< MonitorElement * > cotThetares_vs_pt
std::vector< MonitorElement * > dxyres_vs_eta
std::vector< std::vector< int > > totASS_phi_Quality05
std::vector< MonitorElement * > dzres_vs_pt
key_type key() const
Accessor for product key.
std::vector< MonitorElement * > h_tracks
std::vector< MonitorElement * > dzres_vs_eta
void getRecoMomentum(const reco::Track &track, double &pt, double &ptError, double &qoverp, double &qoverpError, double &lambda, double &lambdaError, double &phi, double &phiError) const
retrieval of reconstructed momentum components from reco::Track (== mean values for GSF) ...
std::vector< std::vector< double > > dxyintervals
std::vector< MonitorElement * > nDThits_vs_eta
std::vector< std::vector< int > > totASS_phi
std::vector< MonitorElement * > thetapull_vs_phi
std::vector< std::vector< int > > totASS_dxy
edm::EDGetTokenT< SimHitTPAssociationProducer::SimHitTPAssociationList > _simHitTpMapTag
std::vector< MonitorElement * > h_assochi2
std::vector< MonitorElement * > chi2_vs_phi
std::vector< std::vector< int > > totSIM_dxy
std::vector< MonitorElement * > chi2_vs_nhits
math::XYZPointD Point
point in the space
std::vector< std::vector< int > > totASSeta_Quality075
std::vector< MonitorElement * > nhits_vs_eta
std::vector< MonitorElement * > phires_vs_phi
virtual double getPt(double pt)
CosmicTrackingParticleSelector cosmictpSelector
std::vector< MonitorElement * > nhits_vs_phi
std::vector< std::vector< int > > totASS2pT
std::vector< MonitorElement * > h_vertposSIM
edm::EDGetTokenT< TrackingParticleCollection > tp_fake_Token
std::vector< MonitorElement * > dxyres_vs_pt
tuple quality
[pTError/pT]*max(1,normChi2) <= ptErrorCut
std::vector< edm::InputTag > label
std::vector< MonitorElement * > h_assocSharedHit
std::vector< MonitorElement * > h_nchi2
std::vector< std::vector< double > > pTintervals
void post_insert()
post insert action
std::vector< std::vector< int > > totASS2eta
std::vector< std::vector< int > > totSIMpT
std::vector< std::vector< int > > totREC_dxy
std::vector< MonitorElement * > h_pullPhi
std::vector< std::vector< int > > totSIM_hit
std::vector< MonitorElement * > h_nchi2_prob
std::vector< MonitorElement * > ptpull_vs_eta
Cos< T >::type cos(const T &t)
std::vector< std::vector< int > > totASSpT
std::vector< MonitorElement * > dxypull_vs_eta
Abs< T >::type abs(const T &t)
std::vector< std::vector< int > > totASS_phi_Quality075
void initEvent(edm::Handle< SimHitTPAssociationProducer::SimHitTPAssociationList > simHitsTPAssocToSet) const
std::vector< MonitorElement * > nME0hits_vs_eta
unsigned short numberOfValidHits() const
number of valid hits found
std::vector< MonitorElement * > h_pullTheta
std::vector< std::vector< int > > totSIM_vertpos
std::vector< std::vector< int > > totASSpT_Quality05
std::vector< std::vector< int > > totREC_hit
T const * get() const
Returns C++ pointer to the item.
std::vector< MonitorElement * > phipull_vs_eta
std::vector< std::vector< int > > totSIMeta
std::vector< std::vector< int > > totASS2_dxy
std::vector< MonitorElement * > thetapull_vs_eta
std::vector< MonitorElement * > h_charge
std::vector< MonitorElement * > ptmean_vs_eta_phi
std::vector< std::vector< int > > totASS2_phi
std::vector< std::vector< double > > phiintervals
edm::EDGetTokenT< reco::RecoToSimCollection > recoToSimCollection_Token
edm::EDGetTokenT< reco::SimToRecoCollection > simToRecoCollection_Token
std::vector< std::vector< int > > totASSpT_Quality075
std::vector< MonitorElement * > h_pullDxy
T const * product() const
std::vector< MonitorElement * > h_tracksSIM
std::vector< MonitorElement * > chi2_vs_eta
TrackingParticleSelector tpSelector
std::vector< std::vector< int > > totSIM_phi
Point vertex() const
Parent vertex position.
std::vector< MonitorElement * > h_hits
std::vector< MonitorElement * > h_losthits
std::vector< MonitorElement * > ptpull_vs_phi
std::vector< std::vector< double > > tmp
edm::InputTag associatormap
std::vector< std::vector< int > > totASSeta
std::vector< MonitorElement * > ptres_vs_pt
std::vector< std::vector< int > > totRECeta
std::vector< std::vector< double > > vertposintervals
std::vector< std::vector< double > > dzintervals
std::vector< MonitorElement * > phires_vs_eta
std::vector< MonitorElement * > h_ptSIM
Monte Carlo truth information used for tracking validation.
std::vector< std::vector< double > > zposintervals
std::vector< MonitorElement * > h_assocFraction
math::XYZVectorD Vector
point in the space
std::vector< std::vector< int > > totSIM_dz
std::vector< std::string > associators
bool BiDirectional_RecoToSim_association
std::vector< MonitorElement * > nlosthits_vs_eta
std::vector< MonitorElement * > dzpull_vs_eta
edm::EDGetTokenT< TrackingParticleCollection > tp_effic_Token
std::vector< MonitorElement * > nRPChits_vs_eta
std::vector< std::vector< int > > totASS_vertpos
std::vector< std::vector< int > > totASS_hit
std::vector< MonitorElement * > nGEMhits_vs_eta
std::vector< MonitorElement * > nCSChits_vs_eta
std::vector< MonitorElement * > phipull_vs_phi
edm::EDGetTokenT< reco::BeamSpot > bsSrc_Token