CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
DuplicateTrackMerger.cc
Go to the documentation of this file.
2 
12 #include "TFile.h"
13 
15 
16 DuplicateTrackMerger::DuplicateTrackMerger(const edm::ParameterSet& iPara) : merger_(iPara)
17 {
18  forestLabel_ = "MVADuplicate";
19  useForestFromDB_ = true;
20 
21  minDeltaR3d_ = -4.0;
22  minBDTG_ = -0.96;
23  minpT_ = 0.2;
24  minP_ = 0.4;
25  maxDCA_ = 30.0;
26  maxDPhi_ = 0.30;
27  maxDLambda_ = 0.30;
28  maxDdsz_ = 10.0;
29  maxDdxy_ = 10.0;
30  maxDQoP_ = 0.25;
31  if(iPara.exists("minpT"))minpT_ = iPara.getParameter<double>("minpT");
32  if(iPara.exists("minP"))minP_ = iPara.getParameter<double>("minP");
33  if(iPara.exists("maxDCA"))maxDCA_ = iPara.getParameter<double>("maxDCA");
34  if(iPara.exists("maxDPhi"))maxDPhi_ = iPara.getParameter<double>("maxDPhi");
35  if(iPara.exists("maxDLambda"))maxDLambda_ = iPara.getParameter<double>("maxDLambda");
36  if(iPara.exists("maxDdsz"))maxDdsz_ = iPara.getParameter<double>("maxDdsz");
37  if(iPara.exists("maxDdxy"))maxDdxy_ = iPara.getParameter<double>("maxDdxy");
38  if(iPara.exists("maxDQoP"))maxDQoP_ = iPara.getParameter<double>("maxDQoP");
39  if(iPara.exists("source"))trackSource_ = consumes<reco::TrackCollection>(iPara.getParameter<edm::InputTag>("source"));
40  if(iPara.exists("minDeltaR3d"))minDeltaR3d_ = iPara.getParameter<double>("minDeltaR3d");
41  if(iPara.exists("minBDTG"))minBDTG_ = iPara.getParameter<double>("minBDTG");
42 
43  produces<std::vector<TrackCandidate> >("candidates");
44  produces<CandidateToDuplicate>("candidateMap");
45 
46  forest_ = 0;
47  gbrVals_ = new float[9];
48  dbFileName_ = "";
49  if(iPara.exists("GBRForestLabel"))forestLabel_ = iPara.getParameter<std::string>("GBRForestLabel");
50  if(iPara.exists("GBRForestFileName")){
51  dbFileName_ = iPara.getParameter<std::string>("GBRForestFileName");
52  useForestFromDB_ = false;
53  }
54 
55  /*
56  tmvaReader_ = new TMVA::Reader("!Color:Silent");
57  tmvaReader_->AddVariable("ddsz",&tmva_ddsz_);
58  tmvaReader_->AddVariable("ddxy",&tmva_ddxy_);
59  tmvaReader_->AddVariable("dphi",&tmva_dphi_);
60  tmvaReader_->AddVariable("dlambda",&tmva_dlambda_);
61  tmvaReader_->AddVariable("dqoverp",&tmva_dqoverp_);
62  tmvaReader_->AddVariable("delta3d_r",&tmva_d3dr_);
63  tmvaReader_->AddVariable("delta3d_z",&tmva_d3dz_);
64  tmvaReader_->AddVariable("outer_nMissingInner",&tmva_outer_nMissingInner_);
65  tmvaReader_->AddVariable("inner_nMissingOuter",&tmva_inner_nMissingOuter_);
66  tmvaReader_->BookMVA("BDTG",mvaFilePath);
67  */
68 
69 }
70 
72 {
73 
74  if(gbrVals_)delete [] gbrVals_;
75  if(!useForestFromDB_ && forest_) delete forest_;
76  /* no op */
77 
78 }
79 
81 {
82 
83  merger_.init(iSetup);
84 
85  if(!forest_){
86  if(useForestFromDB_){
87  edm::ESHandle<GBRForest> forestHandle;
88  iSetup.get<GBRWrapperRcd>().get(forestLabel_,forestHandle);
89  forest_ = (GBRForest*)forestHandle.product();
90  }else{
91  TFile gbrfile(dbFileName_.c_str());
92  forest_ = (GBRForest*)gbrfile.Get(forestLabel_.c_str());
93  }
94  }
95 
96  //edm::Handle<edm::View<reco::Track> >handle;
98  iEvent.getByToken(trackSource_,handle);
99  reco::TrackRefProd refTrks(handle);
100 
101  iSetup.get<IdealMagneticFieldRecord>().get(magfield_);
102  TwoTrackMinimumDistance ttmd;
103  TSCPBuilderNoMaterial tscpBuilder;
104  std::auto_ptr<std::vector<TrackCandidate> > out_duplicateCandidates(new std::vector<TrackCandidate>());
105 
106  std::auto_ptr<CandidateToDuplicate> out_candidateMap(new CandidateToDuplicate());
107 
108  for(int i = 0; i < (int)handle->size(); i++){
109  const reco::Track *rt1 = &(handle->at(i));
110  if(rt1->innerMomentum().Rho() < minpT_)continue;
111  if(rt1->innerMomentum().R() < minP_)continue;
112  for(int j = i+1; j < (int)handle->size();j++){
113  const reco::Track *rt2 = &(handle->at(j));
114  if(rt2->innerMomentum().Rho() < minpT_)continue;
115  if(rt2->innerMomentum().R() < minP_)continue;
116  if(rt1->charge() != rt2->charge())continue;
117  const reco::Track* t1,*t2;
118  if(rt1->outerPosition().Rho() < rt2->outerPosition().Rho()){
119  t1 = rt1;
120  t2 = rt2;
121  }else{
122  t1 = rt2;
123  t2 = rt1;
124  }
125  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));
126 
127  if(t1->outerPosition().Rho() > t2->innerPosition().Rho())deltaR3d *= -1.0;
128  if(deltaR3d < minDeltaR3d_)continue;
129 
132  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);
133  TrajectoryStateClosestToPoint TSCP1 = tscpBuilder(fts1, avgPoint);
134  TrajectoryStateClosestToPoint TSCP2 = tscpBuilder(fts2, avgPoint);
135  if(!TSCP1.isValid())continue;
136  if(!TSCP2.isValid())continue;
137 
138  const FreeTrajectoryState ftsn1 = TSCP1.theState();
139  const FreeTrajectoryState ftsn2 = TSCP2.theState();
140 
141  if ( (ftsn2.position()-ftsn1.position()).mag() > maxDCA_ ) continue;
142 
143  double qoverp1 = ftsn1.signedInverseMomentum();
144  double qoverp2 = ftsn2.signedInverseMomentum();
145  tmva_dqoverp_ = qoverp1-qoverp2;
146  if ( fabs(tmva_dqoverp_) > maxDQoP_ ) continue;
147 
148  double lambda1 = M_PI/2 - ftsn1.momentum().theta();
149  double lambda2 = M_PI/2 - ftsn2.momentum().theta();
150  tmva_dlambda_ = lambda1-lambda2;
151  if ( fabs(tmva_dlambda_) > maxDLambda_ ) continue;
152 
153  double phi1 = ftsn1.momentum().phi();
154  double phi2 = ftsn2.momentum().phi();
155  tmva_dphi_ = phi1-phi2;
156  if(fabs(tmva_dphi_) > M_PI) tmva_dphi_ = 2*M_PI - fabs(tmva_dphi_);
157  if ( fabs(tmva_dphi_) > maxDPhi_ ) continue;
158 
159  double dxy1 = (-ftsn1.position().x() * ftsn1.momentum().y() + ftsn1.position().y() * ftsn1.momentum().x())/TSCP1.pt();
160  double dxy2 = (-ftsn2.position().x() * ftsn2.momentum().y() + ftsn2.position().y() * ftsn2.momentum().x())/TSCP2.pt();
161  tmva_ddxy_ = dxy1-dxy2;
162  if ( fabs(tmva_ddxy_) > maxDdxy_ ) continue;
163 
164  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();
165  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();
166  tmva_ddsz_ = dsz1-dsz2;
167  if ( fabs(tmva_ddsz_) > maxDdsz_ ) continue;
168 
169  tmva_d3dr_ = avgPoint.perp();
170  tmva_d3dz_ = avgPoint.z();
173 
174 
175  gbrVals_[0] = tmva_ddsz_;
176  gbrVals_[1] = tmva_ddxy_;
177  gbrVals_[2] = tmva_dphi_;
178  gbrVals_[3] = tmva_dlambda_;
179  gbrVals_[4] = tmva_dqoverp_;
180  gbrVals_[5] = tmva_d3dr_;
181  gbrVals_[6] = tmva_d3dz_;
184 
185 
186  double mvaBDTG = forest_->GetClassifier(gbrVals_);
187  if(mvaBDTG < minBDTG_)continue;
188 
189 
190  TrackCandidate mergedTrack = merger_.merge(*t1,*t2);
191  out_duplicateCandidates->push_back(mergedTrack);
192  std::pair<TrackRef,TrackRef> trackPair(TrackRef(refTrks,i),TrackRef(refTrks,j));
193  std::pair<TrackCandidate, std::pair<TrackRef,TrackRef> > cp(mergedTrack,trackPair);
194  out_candidateMap->push_back(cp);
195  }
196  }
197  iEvent.put(out_duplicateCandidates,"candidates");
198  iEvent.put(out_candidateMap,"candidateMap");
199 
200 }
T getParameter(std::string const &) const
int i
Definition: DBlmapReader.cc:9
double maxDQoP_
max difference in q/p between two tracks
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:434
const FreeTrajectoryState & theState() const
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
FreeTrajectoryState innerFreeState(const reco::Track &tk, const MagneticField *field)
int numberOfLostHits() const
Definition: HitPattern.h:639
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
T y() const
Definition: PV3DBase.h:63
bool exists(std::string const &parameterName) const
checks if a parameter exists
void produce(edm::Event &, const edm::EventSetup &) override
produce one event
const math::XYZPoint & outerPosition() const
position of the outermost hit
Definition: Track.h:46
double maxDCA_
max distance between two tracks at closest approach
void init(const edm::EventSetup &iSetup)
Definition: TrackMerger.cc:29
TrackCandidate merge(const reco::Track &inner, const reco::Track &outer) const
Definition: TrackMerger.cc:37
Geom::Theta< T > theta() const
Definition: PV3DBase.h:75
const math::XYZPoint & innerPosition() const
position of the innermost hit
Definition: Track.h:41
int iEvent
Definition: GenABIO.cc:243
const HitPattern & trackerExpectedHitsOuter() const
Access the hit pattern counting (in the Tracker) the number of expected crossed layers after the last...
Definition: TrackBase.h:225
T mag() const
Definition: PV3DBase.h:67
std::vector< std::pair< TrackCandidate, std::pair< reco::TrackRef, reco::TrackRef > > > CandidateToDuplicate
auto const T2 &decltype(t1.eta()) t2
Definition: deltaR.h:18
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:116
T sqrt(T t)
Definition: SSEVec.h:48
std::string dbFileName_
MVA weights file.
edm::ESHandle< MagneticField > magfield_
T z() const
Definition: PV3DBase.h:64
GBRForest * forest_
MVA discriminator.
tuple handle
Definition: patZpeak.py:22
int j
Definition: DBlmapReader.cc:9
const HitPattern & trackerExpectedHitsInner() const
Access the hit pattern counting (in the Tracker) the number of expected crossed layers before the fir...
Definition: TrackBase.h:223
float tmva_ddsz_
MVA input variables.
GlobalVector momentum() const
GlobalPoint position() const
#define M_PI
Definition: BFit3D.cc:3
const T & get() const
Definition: EventSetup.h:55
edm::Ref< TrackCollection > TrackRef
persistent reference to a Track
Definition: TrackFwd.h:14
T const * product() const
Definition: ESHandle.h:62
double maxDLambda_
max difference in Lambda between two tracks
double maxDdxy_
max difference in transverse impact parameter between two tracks
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< float >, ROOT::Math::GlobalCoordinateSystemTag > GlobalPoint
point in global coordinate system
Definition: Point3D.h:17
double minDeltaR3d_
minDeltaR3d cut value
const math::XYZVector & innerMomentum() const
momentum vector at the innermost hit position
Definition: Track.h:44
int charge() const
track electric charge
Definition: TrackBase.h:111
FreeTrajectoryState outerFreeState(const reco::Track &tk, const MagneticField *field)
double maxDdsz_
max difference in longitudinal impact parameter between two tracks
double maxDPhi_
max difference in phi between two tracks
double GetClassifier(const float *vector) const
Definition: GBRForest.h:60
T x() const
Definition: PV3DBase.h:62
double signedInverseMomentum() const
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
edm::EDGetTokenT< reco::TrackCollection > trackSource_
track input collection