CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/RecoTracker/FinalTrackSelectors/src/DuplicateTrackMerger.cc

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   tmvaReader_ = new TMVA::Reader("!Color:Silent");
00057   tmvaReader_->AddVariable("ddsz",&tmva_ddsz_);
00058   tmvaReader_->AddVariable("ddxy",&tmva_ddxy_);
00059   tmvaReader_->AddVariable("dphi",&tmva_dphi_);
00060   tmvaReader_->AddVariable("dlambda",&tmva_dlambda_);
00061   tmvaReader_->AddVariable("dqoverp",&tmva_dqoverp_);
00062   tmvaReader_->AddVariable("delta3d_r",&tmva_d3dr_);
00063   tmvaReader_->AddVariable("delta3d_z",&tmva_d3dz_);
00064   tmvaReader_->AddVariable("outer_nMissingInner",&tmva_outer_nMissingInner_);
00065   tmvaReader_->AddVariable("inner_nMissingOuter",&tmva_inner_nMissingOuter_);
00066   tmvaReader_->BookMVA("BDTG",mvaFilePath);
00067   */
00068 
00069 }
00070 
00071 DuplicateTrackMerger::~DuplicateTrackMerger()
00072 {
00073 
00074   if(gbrVals_)delete [] gbrVals_;
00075   if(!useForestFromDB_ && forest_) delete forest_;
00076   /* no op */
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   //edm::Handle<edm::View<reco::Track> >handle;
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 }