Go to the documentation of this file.00001 #include "RecoTracker/FinalTrackSelectors/interface/DuplicateTrackMerger.h"
00002
00003 #include "FWCore/Framework/interface/Event.h"
00004 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
00005 #include "TrackingTools/PatternTools/interface/TSCPBuilderNoMaterial.h"
00006 #include "TrackingTools/PatternTools/interface/TwoTrackMinimumDistance.h"
00007 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00008 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
00009 #include "CondFormats/DataRecord/interface/GBRWrapperRcd.h"
00010 #include "FWCore/Framework/interface/EventSetup.h"
00011 #include "FWCore/Framework/interface/ESHandle.h"
00012 #include "TFile.h"
00013
00014 using reco::modules::DuplicateTrackMerger;
00015
00016 DuplicateTrackMerger::DuplicateTrackMerger(const edm::ParameterSet& iPara) : merger_(iPara)
00017 {
00018 forestLabel_ = "MVADuplicate";
00019 useForestFromDB_ = true;
00020
00021 minDeltaR3d_ = -4.0;
00022 minBDTG_ = -0.96;
00023 minpT_ = 0.2;
00024 minP_ = 0.4;
00025 maxDCA_ = 30.0;
00026 maxDPhi_ = 0.30;
00027 maxDLambda_ = 0.30;
00028 maxDdsz_ = 10.0;
00029 maxDdxy_ = 10.0;
00030 maxDQoP_ = 0.25;
00031 if(iPara.exists("minpT"))minpT_ = iPara.getParameter<double>("minpT");
00032 if(iPara.exists("minP"))minP_ = iPara.getParameter<double>("minP");
00033 if(iPara.exists("maxDCA"))maxDCA_ = iPara.getParameter<double>("maxDCA");
00034 if(iPara.exists("maxDPhi"))maxDPhi_ = iPara.getParameter<double>("maxDPhi");
00035 if(iPara.exists("maxDLambda"))maxDLambda_ = iPara.getParameter<double>("maxDLambda");
00036 if(iPara.exists("maxDdsz"))maxDdsz_ = iPara.getParameter<double>("maxDdsz");
00037 if(iPara.exists("maxDdxy"))maxDdxy_ = iPara.getParameter<double>("maxDdxy");
00038 if(iPara.exists("maxDQoP"))maxDQoP_ = iPara.getParameter<double>("maxDQoP");
00039 if(iPara.exists("source"))trackSource_ = iPara.getParameter<edm::InputTag>("source");
00040 if(iPara.exists("minDeltaR3d"))minDeltaR3d_ = iPara.getParameter<double>("minDeltaR3d");
00041 if(iPara.exists("minBDTG"))minBDTG_ = iPara.getParameter<double>("minBDTG");
00042
00043 produces<std::vector<TrackCandidate> >("candidates");
00044 produces<CandidateToDuplicate>("candidateMap");
00045
00046 forest_ = 0;
00047 gbrVals_ = new float[9];
00048 dbFileName_ = "";
00049 if(iPara.exists("GBRForestLabel"))forestLabel_ = iPara.getParameter<std::string>("GBRForestLabel");
00050 if(iPara.exists("GBRForestFileName")){
00051 dbFileName_ = iPara.getParameter<std::string>("GBRForestFileName");
00052 useForestFromDB_ = false;
00053 }
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069 }
00070
00071 DuplicateTrackMerger::~DuplicateTrackMerger()
00072 {
00073
00074 if(gbrVals_)delete [] gbrVals_;
00075 if(!useForestFromDB_ && forest_) delete forest_;
00076
00077
00078 }
00079
00080 void DuplicateTrackMerger::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
00081 {
00082
00083 merger_.init(iSetup);
00084
00085 if(!forest_){
00086 if(useForestFromDB_){
00087 edm::ESHandle<GBRForest> forestHandle;
00088 iSetup.get<GBRWrapperRcd>().get(forestLabel_,forestHandle);
00089 forest_ = (GBRForest*)forestHandle.product();
00090 }else{
00091 TFile gbrfile(dbFileName_.c_str());
00092 forest_ = (GBRForest*)gbrfile.Get(forestLabel_.c_str());
00093 }
00094 }
00095
00096
00097 edm::Handle<reco::TrackCollection >handle;
00098 iEvent.getByLabel(trackSource_,handle);
00099 reco::TrackRefProd refTrks(handle);
00100
00101 iSetup.get<IdealMagneticFieldRecord>().get(magfield_);
00102 TwoTrackMinimumDistance ttmd;
00103 TSCPBuilderNoMaterial tscpBuilder;
00104 std::auto_ptr<std::vector<TrackCandidate> > out_duplicateCandidates(new std::vector<TrackCandidate>());
00105
00106 std::auto_ptr<CandidateToDuplicate> out_candidateMap(new CandidateToDuplicate());
00107
00108 for(int i = 0; i < (int)handle->size(); i++){
00109 const reco::Track *rt1 = &(handle->at(i));
00110 if(rt1->innerMomentum().Rho() < minpT_)continue;
00111 if(rt1->innerMomentum().R() < minP_)continue;
00112 for(int j = i+1; j < (int)handle->size();j++){
00113 const reco::Track *rt2 = &(handle->at(j));
00114 if(rt2->innerMomentum().Rho() < minpT_)continue;
00115 if(rt2->innerMomentum().R() < minP_)continue;
00116 if(rt1->charge() != rt2->charge())continue;
00117 const reco::Track* t1,*t2;
00118 if(rt1->outerPosition().Rho() < rt2->outerPosition().Rho()){
00119 t1 = rt1;
00120 t2 = rt2;
00121 }else{
00122 t1 = rt2;
00123 t2 = rt1;
00124 }
00125 double deltaR3d = sqrt(pow(t1->outerPosition().x()-t2->innerPosition().x(),2) + pow(t1->outerPosition().y()-t2->innerPosition().y(),2) + pow(t1->outerPosition().z()-t2->innerPosition().z(),2));
00126
00127 if(t1->outerPosition().Rho() > t2->innerPosition().Rho())deltaR3d *= -1.0;
00128 if(deltaR3d < minDeltaR3d_)continue;
00129
00130 FreeTrajectoryState fts1 = trajectoryStateTransform::outerFreeState(*t1, &*magfield_);
00131 FreeTrajectoryState fts2 = trajectoryStateTransform::innerFreeState(*t2, &*magfield_);
00132 GlobalPoint avgPoint((t1->outerPosition().x()+t2->innerPosition().x())*0.5,(t1->outerPosition().y()+t2->innerPosition().y())*0.5,(t1->outerPosition().z()+t2->innerPosition().z())*0.5);
00133 TrajectoryStateClosestToPoint TSCP1 = tscpBuilder(fts1, avgPoint);
00134 TrajectoryStateClosestToPoint TSCP2 = tscpBuilder(fts2, avgPoint);
00135 if(!TSCP1.isValid())continue;
00136 if(!TSCP2.isValid())continue;
00137
00138 const FreeTrajectoryState ftsn1 = TSCP1.theState();
00139 const FreeTrajectoryState ftsn2 = TSCP2.theState();
00140
00141 if ( (ftsn2.position()-ftsn1.position()).mag() > maxDCA_ ) continue;
00142
00143 double qoverp1 = ftsn1.signedInverseMomentum();
00144 double qoverp2 = ftsn2.signedInverseMomentum();
00145 tmva_dqoverp_ = qoverp1-qoverp2;
00146 if ( fabs(tmva_dqoverp_) > maxDQoP_ ) continue;
00147
00148 double lambda1 = M_PI/2 - ftsn1.momentum().theta();
00149 double lambda2 = M_PI/2 - ftsn2.momentum().theta();
00150 tmva_dlambda_ = lambda1-lambda2;
00151 if ( fabs(tmva_dlambda_) > maxDLambda_ ) continue;
00152
00153 double phi1 = ftsn1.momentum().phi();
00154 double phi2 = ftsn2.momentum().phi();
00155 tmva_dphi_ = phi1-phi2;
00156 if(fabs(tmva_dphi_) > M_PI) tmva_dphi_ = 2*M_PI - fabs(tmva_dphi_);
00157 if ( fabs(tmva_dphi_) > maxDPhi_ ) continue;
00158
00159 double dxy1 = (-ftsn1.position().x() * ftsn1.momentum().y() + ftsn1.position().y() * ftsn1.momentum().x())/TSCP1.pt();
00160 double dxy2 = (-ftsn2.position().x() * ftsn2.momentum().y() + ftsn2.position().y() * ftsn2.momentum().x())/TSCP2.pt();
00161 tmva_ddxy_ = dxy1-dxy2;
00162 if ( fabs(tmva_ddxy_) > maxDdxy_ ) continue;
00163
00164 double dsz1 = ftsn1.position().z() * TSCP1.pt() / TSCP1.momentum().mag() - (ftsn1.position().x() * ftsn1.momentum().y() + ftsn1.position().y() * ftsn1.momentum().x())/TSCP1.pt() * ftsn1.momentum().z()/ftsn1.momentum().mag();
00165 double dsz2 = ftsn2.position().z() * TSCP2.pt() / TSCP2.momentum().mag() - (ftsn2.position().x() * ftsn2.momentum().y() + ftsn2.position().y() * ftsn2.momentum().x())/TSCP2.pt() * ftsn2.momentum().z()/ftsn2.momentum().mag();
00166 tmva_ddsz_ = dsz1-dsz2;
00167 if ( fabs(tmva_ddsz_) > maxDdsz_ ) continue;
00168
00169 tmva_d3dr_ = avgPoint.perp();
00170 tmva_d3dz_ = avgPoint.z();
00171 tmva_outer_nMissingInner_ = t2->trackerExpectedHitsInner().numberOfLostHits();
00172 tmva_inner_nMissingOuter_ = t1->trackerExpectedHitsOuter().numberOfLostHits();
00173
00174
00175 gbrVals_[0] = tmva_ddsz_;
00176 gbrVals_[1] = tmva_ddxy_;
00177 gbrVals_[2] = tmva_dphi_;
00178 gbrVals_[3] = tmva_dlambda_;
00179 gbrVals_[4] = tmva_dqoverp_;
00180 gbrVals_[5] = tmva_d3dr_;
00181 gbrVals_[6] = tmva_d3dz_;
00182 gbrVals_[7] = tmva_outer_nMissingInner_;
00183 gbrVals_[8] = tmva_inner_nMissingOuter_;
00184
00185
00186 double mvaBDTG = forest_->GetClassifier(gbrVals_);
00187 if(mvaBDTG < minBDTG_)continue;
00188
00189
00190 TrackCandidate mergedTrack = merger_.merge(*t1,*t2);
00191 out_duplicateCandidates->push_back(mergedTrack);
00192 std::pair<TrackRef,TrackRef> trackPair(TrackRef(refTrks,i),TrackRef(refTrks,j));
00193 std::pair<TrackCandidate, std::pair<TrackRef,TrackRef> > cp(mergedTrack,trackPair);
00194 out_candidateMap->push_back(cp);
00195 }
00196 }
00197 iEvent.put(out_duplicateCandidates,"candidates");
00198 iEvent.put(out_candidateMap,"candidateMap");
00199
00200 }