ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitProd/TreeFiller/src/FillerDCASig.cc
Revision: 1.2
Committed: Wed May 23 03:24:07 2012 UTC (12 years, 11 months ago) by paus
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_032, Mit_031, Mit_025c_branch2, Mit_025c_branch1, Mit_030, Mit_029c, Mit_029b, Mit_030_pre1, Mit_029a, Mit_029, Mit_029_pre1, Mit_028a, Mit_025c_branch0, HEAD
Branch point for: Mit_025c_branch
Changes since 1.1: +352 -276 lines
Log Message: 
Backport.

File Contents

# User Rev Content
1 pharris 1.1 #include "MitProd/TreeFiller/interface/FillerDCASig.h"
2     #include "DataFormats/Common/interface/RefToPtr.h"
3     #include "MitEdm/DataFormats/interface/RefToBaseToPtr.h"
4     #include "MitAna/DataTree/interface/DCASigCol.h"
5     #include "MitAna/DataTree/interface/Names.h"
6     #include "MitProd/ObjectService/interface/ObjectService.h"
7    
8     #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
9     #include "TrackingTools/Records/interface/TransientTrackRecord.h"
10     #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
11     #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
12     #include "TrackingTools/PatternTools/interface/Trajectory.h"
13     #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
14     #include "TrackingTools/PatternTools/interface/ClosestApproachInRPhi.h"
15     #include "TrackingTools/PatternTools/interface/TwoTrackMinimumDistance.h"
16     #include "RecoVertex/KinematicFitPrimitives/interface/KinematicParticleFactoryFromTransientTrack.h"
17    
18     using namespace std;
19     using namespace edm;
20     using namespace mithep;
21    
22     //--------------------------------------------------------------------------------------------------
23     FillerDCASig::FillerDCASig(const ParameterSet &cfg, const char *name, bool active) :
24 paus 1.2 BaseFiller (cfg,name,active),
25     edmElectronName_ (Conf().getUntrackedParameter<string>("edmElectronName","gsfElectrons")),
26     edmMuonName_ (Conf().getUntrackedParameter<string>("edmMuonName", "muons")),
27     edmTauName_ (Conf().getUntrackedParameter<string>("edmTauName", "hpsPFTauProducer")),
28     mitName_ (Conf().getUntrackedParameter<string>("mitName", "DCASig")),
29     electronMapName_ (Conf().getUntrackedParameter<string>("electronMapName","")),
30     muonMapName_ (Conf().getUntrackedParameter<string>("muonMapName", "")),
31     tauMapName_ (Conf().getUntrackedParameter<string>("tauMapName", "")),
32     electronPtMin_ (Conf().getUntrackedParameter<double>("electronPtMin", 5.)),
33     muonPtMin_ (Conf().getUntrackedParameter<double>("muonPtMin", 5.)),
34     tauPtMin_ (Conf().getUntrackedParameter<double>("tauPtMin", 10.)),
35     checkDCARef_ (Conf().getUntrackedParameter<bool>("checkDCARef", "False")),
36     DCASigs_ (new mithep::DCASigArr),
37 pharris 1.1 transientTrackBuilder_(0)
38     {
39     }
40    
41     //--------------------------------------------------------------------------------------------------
42     FillerDCASig::~FillerDCASig()
43     {
44     // Destructor.
45     delete DCASigs_;
46     }
47    
48     //--------------------------------------------------------------------------------------------------
49     void FillerDCASig::BookDataBlock(TreeWriter &tws)
50     {
51     // Add DCA to the tree
52     tws.AddBranch(mitName_,&DCASigs_);
53     OS()->add<mithep::DCASigArr>(DCASigs_,mitName_);
54 paus 1.2
55     //Load lepton maps
56     if (!electronMapName_.empty()) {
57     electronMap_ = OS()->get<ElectronMap>(electronMapName_);
58     if (electronMap_)
59     AddBranchDep(mitName_,electronMap_->GetBrName());
60     else
61     printf("\n FillerDCASig::BookDataBlock Electron map (Name: %s) not found!!\n\n",electronMapName_.data());
62     }
63 pharris 1.1 if (!muonMapName_.empty()) {
64     muonMap_ = OS()->get<MuonMap>(muonMapName_);
65     if (muonMap_)
66     AddBranchDep(mitName_,muonMap_->GetBrName());
67 paus 1.2 else
68     printf("\n FillerDCASig::BookDataBlock Muon map (Name: %s) not found!!\n\n",muonMapName_.data());
69 pharris 1.1 }
70     if (!tauMapName_.empty()) {
71     tauMap_ = OS()->get<PFTauMap>(tauMapName_);
72     if (tauMap_)
73     AddBranchDep(mitName_,tauMap_->GetBrName());
74 paus 1.2 else
75     printf("\n FillerDCASig::BookDataBlock Tau map (Name: %s) not found!!\n\n",tauMapName_.data());
76 pharris 1.1 }
77     }
78     //--------------------------------------------------------------------------------------------------
79 paus 1.2 void FillerDCASig::FillDataBlock(const edm::Event &event,
80     const edm::EventSetup &setup)
81 pharris 1.1 {
82     // Fill conversions data structure and maps.
83    
84     DCASigs_->Delete();
85    
86 paus 1.2 //Transient Track Buildder
87 pharris 1.1 edm::ESHandle<TransientTrackBuilder> builder;
88     setup.get<TransientTrackRecord>().get("TransientTrackBuilder",builder);
89     transientTrackBuilder_ = builder.product();
90 paus 1.2
91 pharris 1.1 //handle for the electron collection
92     Handle<reco::GsfElectronCollection> hElectronProduct;
93     GetProduct(edmElectronName_, hElectronProduct, event);
94 paus 1.2 const reco::GsfElectronCollection inElectrons = *(hElectronProduct.product());
95 pharris 1.1
96     //handle for the muon collection
97     Handle<reco::MuonCollection> hMuonProduct;
98 paus 1.2 GetProduct(edmMuonName_, hMuonProduct, event);
99     const reco::MuonCollection inMuons = *(hMuonProduct.product());
100 pharris 1.1
101     // handle for the tau collection
102     Handle<reco::PFTauCollection> hTauProduct;
103     GetProduct(edmTauName_, hTauProduct, event);
104 paus 1.2 const reco::PFTauCollection inTaus = *(hTauProduct.product());
105    
106 pharris 1.1 //Declare some variables
107     Double_t lDCA3D = 0; //3D x-y DCA
108     Double_t lDCA3DE = 0; //3D x-y DCA Err
109     Double_t lDCA2D = 0; //2D x-y DCA
110     Double_t lDCA2DE = 0; //2D x-y DCA Err
111     Double_t lDCARPhi3D = 0; //3D r-phi DCA
112     Double_t lDCARPhi3DE = 0; //3D r-phi DCA Err
113     Double_t lDCARPhi2D = 0; //2D r-phi DCA
114     Double_t lDCARPhi2DE = 0; //2D r-phi DCA Err
115 paus 1.2
116     //Loop through the electrons
117     for (reco::GsfElectronCollection::const_iterator iElectron = inElectrons.begin();
118     iElectron != inElectrons.end(); ++iElectron) {
119     if (iElectron->pt() < electronPtMin_)
120     continue;
121    
122     //E + Mu Objects
123     for (reco::MuonCollection::const_iterator iMuon = inMuons.begin(); iMuon != inMuons.end(); ++iMuon) {
124     if (iMuon->pt() < muonPtMin_)
125     continue;
126    
127 pharris 1.1 //Declare the Object
128 paus 1.2 const reco::Track * trackElec = (const reco::Track*)(iElectron->gsfTrack().get()); // electron Track
129     const reco::Track * trackMuon = (const reco::Track*)(iMuon->innerTrack() .get()); // muon inner track
130     if (trackMuon == NULL || trackElec == NULL || trackElec == trackMuon) continue;
131 pharris 1.1 mithep::DCASig *outDCASig = DCASigs_->AddNew();
132     calculateDCA(lDCA3D ,lDCA3DE ,lDCA2D ,lDCA2DE,
133 paus 1.2 lDCARPhi3D,lDCARPhi3DE,lDCARPhi2D,lDCARPhi2DE,
134     trackElec,trackMuon,DCASig::eEMu);
135    
136     outDCASig->SetDCA2D (lDCA2D); outDCASig->SetDCA2DErr (lDCA2DE);
137     outDCASig->SetDCA3D (lDCA3D); outDCASig->SetDCA3DErr (lDCA3DE);
138     outDCASig->SetDCA2DRPhi(lDCARPhi2D); outDCASig->SetDCA2DRPhiErr(lDCARPhi2DE);
139     outDCASig->SetDCA3DRPhi(lDCARPhi3D); outDCASig->SetDCA3DRPhiErr(lDCARPhi3DE);
140 pharris 1.1 outDCASig->SetType(DCASig::eEMu);
141    
142     //Now the references
143     edm::Ptr<reco::GsfElectron> pElectronPtr(hElectronProduct, iElectron - inElectrons.begin());
144 paus 1.2 try {
145     outDCASig->SetElectron(electronMap_->GetMit(pElectronPtr));
146     }
147     catch(...) {
148     if (checkDCARef_) throw edm::Exception(edm::errors::Configuration, "FillerDCASig:FillDataBlock()\n")
149     << "Error! Electron unmapped collection " << edmElectronName_ << endl;
150 pharris 1.1 }
151     edm::Ptr<reco::Muon> pMuonPtr(hMuonProduct, iMuon - inMuons.begin());
152 paus 1.2 try {
153     outDCASig->SetMuon (muonMap_->GetMit(pMuonPtr));
154     }
155     catch(...) {
156     if (checkDCARef_) throw edm::Exception(edm::errors::Configuration, "FillerDCASig:FillDataBlock()\n")
157     << "Error! Muon unmapped collection " << edmMuonName_ << endl;
158 pharris 1.1 }
159     }
160    
161 paus 1.2 //E + Tau objects
162     for (reco::PFTauCollection::const_iterator iTau = inTaus.begin(); iTau != inTaus.end(); ++iTau) {
163     if (iTau->pt() < tauPtMin_)
164     continue;
165    
166     //declare the Object
167     const reco::Track * trackElec = (const reco::Track*)(iElectron->gsfTrack().get()); // electron Track
168     if (iTau->leadPFChargedHadrCand().isNull())
169     continue;
170    
171     const reco::Track * trackTau = (const reco::Track*)(iTau->leadPFChargedHadrCand()->trackRef().get()); // Tau lead track
172     if (trackTau == NULL || trackElec == NULL || trackTau == trackElec)
173     continue;
174    
175     mithep::DCASig *outDCASig = DCASigs_->AddNew();
176 pharris 1.1 calculateDCA(lDCA3D ,lDCA3DE ,lDCA2D ,lDCA2DE,
177 paus 1.2 lDCARPhi3D,lDCARPhi3DE,lDCARPhi2D,lDCARPhi2DE,
178     trackElec,trackTau,DCASig::eETau);
179 pharris 1.1
180 paus 1.2 outDCASig->SetDCA2D (lDCA2D); outDCASig->SetDCA2DErr (lDCA2DE);
181     outDCASig->SetDCA3D (lDCA3D); outDCASig->SetDCA3DErr (lDCA3DE);
182     outDCASig->SetDCA2DRPhi(lDCARPhi2D); outDCASig->SetDCA2DRPhiErr(lDCARPhi2DE);
183     outDCASig->SetDCA3DRPhi(lDCARPhi3D); outDCASig->SetDCA3DRPhiErr(lDCARPhi3DE);
184 pharris 1.1 outDCASig->SetType(DCASig::eETau);
185    
186 paus 1.2 //references
187     edm::Ptr<reco::GsfElectron> pElectronPtr(hElectronProduct, iElectron-inElectrons.begin());
188     try {
189     outDCASig->SetElectron(electronMap_->GetMit(pElectronPtr));
190     }
191     catch(...) {
192     if (checkDCARef_) throw edm::Exception(edm::errors::Configuration, "FillerDCASig:FillDataBlock()\n")
193     << "Error! Electron unmapped collection " << edmElectronName_ << endl;
194     }
195     edm::Ptr<reco::PFTau> pTauPtr(hTauProduct,iTau - inTaus.begin());
196     try {
197     outDCASig->SetTau(tauMap_->GetMit(pTauPtr));
198     }
199     catch(...) {
200     if (checkDCARef_)
201     throw edm::Exception(edm::errors::Configuration, "FillerDCASig:FillDataBlock()\n")
202     << "Error! Tau unmapped collection " << edmTauName_ << endl;
203 pharris 1.1 }
204     //std::cout << " ===> Tau Ref : " << &pTauPtr << " -- " << ((tauMap_->GetMit(pTauPtr))) << endl;
205     }
206 paus 1.2
207 pharris 1.1 //E + E Objects
208     for (reco::GsfElectronCollection::const_iterator iElectron1 = iElectron+1;
209 paus 1.2 iElectron1 != inElectrons.end(); ++iElectron1) {
210     if (iElectron1->pt() < electronPtMin_)
211     continue;
212    
213     //declare the Object
214     const reco::Track * trackElec0 = (const reco::Track*)(iElectron ->gsfTrack().get()); // electron Track
215     const reco::Track * trackElec1 = (const reco::Track*)(iElectron1->gsfTrack().get()); // electron Track
216 pharris 1.1
217 paus 1.2 if (trackElec1 == NULL || trackElec0 == NULL || trackElec0 == trackElec1)
218     continue;
219 pharris 1.1
220 paus 1.2 mithep::DCASig *outDCASig = DCASigs_->AddNew();
221 pharris 1.1 calculateDCA(lDCA3D ,lDCA3DE ,lDCA2D ,lDCA2DE,
222 paus 1.2 lDCARPhi3D,lDCARPhi3DE,lDCARPhi2D,lDCARPhi2DE,
223     trackElec0,trackElec1,DCASig::eEE);
224 pharris 1.1
225 paus 1.2 outDCASig->SetDCA2D (lDCA2D); outDCASig->SetDCA2DErr (lDCA2DE);
226     outDCASig->SetDCA3D (lDCA3D); outDCASig->SetDCA3DErr (lDCA3DE);
227     outDCASig->SetDCA2DRPhi(lDCARPhi2D); outDCASig->SetDCA2DRPhiErr(lDCARPhi2DE);
228     outDCASig->SetDCA3DRPhi(lDCARPhi3D); outDCASig->SetDCA3DRPhiErr(lDCARPhi3DE);
229 pharris 1.1 outDCASig->SetType(DCASig::eEE);
230 paus 1.2
231     //references
232     edm::Ptr<reco::GsfElectron> pElectronPtr(hElectronProduct, iElectron - inElectrons.begin());
233     try {
234     outDCASig->SetElectron(electronMap_->GetMit(pElectronPtr));
235     }
236     catch(...) {
237     if (checkDCARef_)
238     throw edm::Exception(edm::errors::Configuration, "FillerDCASig:FillDataBlock()\n")
239     << "Error! Electron unmapped collection " << edmElectronName_ << endl;
240 pharris 1.1 }
241     edm::Ptr<reco::GsfElectron> pElectron1Ptr(hElectronProduct, iElectron1 - inElectrons.begin());
242 paus 1.2 try {
243     outDCASig->SetElectron(electronMap_->GetMit(pElectron1Ptr),true);
244     }
245     catch(...) {
246     if (checkDCARef_)
247     throw edm::Exception(edm::errors::Configuration, "FillerDCASig:FillDataBlock()\n")
248     << "Error! Electron unmapped collection " << edmElectronName_ << endl;
249 pharris 1.1 }
250     }
251     }
252    
253 paus 1.2 // Mu + Tau combination
254     for (reco::MuonCollection::const_iterator iMuon = inMuons.begin(); iMuon != inMuons.end();
255     ++iMuon) {
256     if (iMuon->pt() < muonPtMin_)
257     continue;
258    
259     for (reco::PFTauCollection::const_iterator iTau = inTaus.begin();
260     iTau != inTaus.end(); ++iTau) {
261     if (iTau->pt() < tauPtMin_)
262     continue;
263    
264     const reco::Track * trackMuon = (const reco::Track*)(iMuon->innerTrack().get()); // muon inner track
265     if (iTau->leadPFChargedHadrCand().isNull())
266     continue;
267     const reco::Track * trackTau = (const reco::Track*)(iTau->leadPFChargedHadrCand()->trackRef().get()); // Tau lead track
268    
269     if (trackTau == NULL || trackMuon == NULL || trackTau == trackMuon)
270     continue;
271    
272 pharris 1.1 calculateDCA(lDCA3D ,lDCA3DE ,lDCA2D ,lDCA2DE,
273 paus 1.2 lDCARPhi3D,lDCARPhi3DE,lDCARPhi2D,lDCARPhi2DE,
274     trackMuon,trackTau,DCASig::eMuTau);
275 pharris 1.1
276 paus 1.2 mithep::DCASig *outDCASig = DCASigs_->AddNew();
277     outDCASig->SetDCA2D (lDCA2D); outDCASig->SetDCA2DErr (lDCA2DE);
278     outDCASig->SetDCA3D (lDCA3D); outDCASig->SetDCA3DErr (lDCA3DE);
279     outDCASig->SetDCA2DRPhi(lDCARPhi2D); outDCASig->SetDCA2DRPhiErr(lDCARPhi2DE);
280     outDCASig->SetDCA3DRPhi(lDCARPhi3D); outDCASig->SetDCA3DRPhiErr(lDCARPhi3DE);
281 pharris 1.1 outDCASig->SetType(DCASig::eMuTau);
282    
283 paus 1.2 //references
284     edm::Ptr<reco::Muon> pMuonPtr(hMuonProduct, iMuon - inMuons.begin());
285     try {
286     outDCASig->SetMuon(muonMap_->GetMit(pMuonPtr));
287     }
288     catch(...) {
289     if (checkDCARef_)
290     throw edm::Exception(edm::errors::Configuration, "FillerDCASig:FillDataBlock()\n")
291     << "Error! Muon unmapped collection " << edmMuonName_ << endl;
292     }
293     edm::Ptr<reco::PFTau> pTauPtr(hTauProduct, iTau - inTaus.begin());
294     try {
295     outDCASig->SetTau(tauMap_->GetMit(pTauPtr));
296     }
297     catch(...) {
298     if (checkDCARef_)
299     throw edm::Exception(edm::errors::Configuration, "FillerDCASig:FillDataBlock()\n")
300     << "Error! Tau unmapped collection " << edmTauName_ << endl;
301 pharris 1.1 }
302     }
303 paus 1.2
304     //Mu + Mu
305 pharris 1.1 for (reco::MuonCollection::const_iterator iMuon1 = iMuon+1;
306 paus 1.2 iMuon1 != inMuons.end(); ++iMuon1) {
307     if (iMuon1->pt() < muonPtMin_)
308     continue;
309 pharris 1.1
310 paus 1.2 const reco::Track * trackMuon0 = (const reco::Track*)(iMuon ->innerTrack().get()); // muon inner track
311     const reco::Track * trackMuon1 = (const reco::Track*)(iMuon1->innerTrack().get()); // muon inner track
312    
313     if (trackMuon1 == NULL || trackMuon0 == NULL || trackMuon0 == trackMuon1)
314     continue;
315 pharris 1.1
316     calculateDCA(lDCA3D ,lDCA3DE ,lDCA2D ,lDCA2DE,
317 paus 1.2 lDCARPhi3D,lDCARPhi3DE,lDCARPhi2D,lDCARPhi2DE,
318     trackMuon0,trackMuon1,DCASig::eMuMu);
319    
320     mithep::DCASig *outDCASig = DCASigs_->AddNew();
321     outDCASig->SetDCA2D (lDCA2D); outDCASig->SetDCA2DErr (lDCA2DE);
322     outDCASig->SetDCA3D (lDCA3D); outDCASig->SetDCA3DErr (lDCA3DE);
323     outDCASig->SetDCA2DRPhi(lDCARPhi2D); outDCASig->SetDCA2DRPhiErr(lDCARPhi2DE);
324     outDCASig->SetDCA3DRPhi(lDCARPhi3D); outDCASig->SetDCA3DRPhiErr(lDCARPhi3DE);
325 pharris 1.1 outDCASig->SetType(DCASig::eMuMu);
326 paus 1.2
327 pharris 1.1 //Now the references
328 paus 1.2 edm::Ptr<reco::Muon> pMuon0Ptr(hMuonProduct, iMuon - inMuons.begin());
329     try {
330     outDCASig->SetMuon(muonMap_->GetMit(pMuon0Ptr));
331     }
332     catch(...) {
333     if (checkDCARef_)
334     throw edm::Exception(edm::errors::Configuration, "FillerDCASig:FillDataBlock()\n")
335     << "Error! Muon unmapped collection " << edmMuonName_ << endl;
336     }
337     edm::Ptr<reco::Muon> pMuon1Ptr(hMuonProduct, iMuon1 - inMuons.begin());
338     try {
339     outDCASig->SetMuon(muonMap_->GetMit(pMuon1Ptr),true);
340     }
341     catch(...) {
342     if (checkDCARef_)
343     throw edm::Exception(edm::errors::Configuration, "FillerDCASig:FillDataBlock()\n")
344     << "Error! Muon unmapped collection " << edmMuonName_ << endl;
345 pharris 1.1 }
346     }
347     }
348 paus 1.2
349     // Tau + Tau combination
350     for (reco::PFTauCollection::const_iterator iTau0 = inTaus.begin();
351     iTau0 != inTaus.end(); ++iTau0) {
352     if (iTau0->pt() < tauPtMin_)
353     continue;
354 pharris 1.1 for (reco::PFTauCollection::const_iterator iTau1 = iTau0+1;
355 paus 1.2 iTau1 != inTaus.end(); ++iTau1) {
356     if (iTau0 == iTau1)
357     continue;
358     if (iTau1->pt() < tauPtMin_)
359     continue;
360    
361     if (iTau0->leadPFChargedHadrCand().isNull())
362     continue;
363     if (iTau1->leadPFChargedHadrCand().isNull())
364     continue;
365     const reco::Track *trackTau0 = (const reco::Track*)(iTau0->leadPFChargedHadrCand()->trackRef().get()); // Tau lead track
366     const reco::Track *trackTau1 = (const reco::Track*)(iTau1->leadPFChargedHadrCand()->trackRef().get()); // Tau lead track
367    
368     if (trackTau1 == NULL || trackTau0 == NULL || trackTau0 == trackTau1)
369     continue;
370    
371 pharris 1.1 calculateDCA(lDCA3D ,lDCA3DE ,lDCA2D ,lDCA2DE,
372 paus 1.2 lDCARPhi3D,lDCARPhi3DE,lDCARPhi2D,lDCARPhi2DE,
373     trackTau0,trackTau1,DCASig::eTauTau);
374 pharris 1.1
375 paus 1.2 mithep::DCASig *outDCASig = DCASigs_->AddNew();
376 pharris 1.1
377 paus 1.2 outDCASig->SetDCA2D (lDCA2D); outDCASig->SetDCA2DErr (lDCA2DE);
378     outDCASig->SetDCA3D (lDCA3D); outDCASig->SetDCA3DErr (lDCA3DE);
379     outDCASig->SetDCA2DRPhi(lDCARPhi2D); outDCASig->SetDCA2DRPhiErr(lDCARPhi2DE);
380     outDCASig->SetDCA3DRPhi(lDCARPhi3D); outDCASig->SetDCA3DRPhiErr(lDCARPhi3DE);
381 pharris 1.1 outDCASig->SetType(DCASig::eTauTau);
382 paus 1.2
383     //references
384     edm::Ptr<reco::PFTau> pTau0Ptr(hTauProduct, iTau0 - inTaus.begin());
385     try {
386     outDCASig->SetTau(tauMap_->GetMit(pTau0Ptr));}
387     catch(...) {
388     if (checkDCARef_)
389     throw edm::Exception(edm::errors::Configuration, "FillerDCASig:FillDataBlock()\n")
390     << "Error! Tau unmapped collection " << edmTauName_ << endl;
391     }
392     edm::Ptr<reco::PFTau> pTau1Ptr(hTauProduct, iTau1 - inTaus.begin());
393     try {
394     outDCASig->SetTau(tauMap_->GetMit(pTau1Ptr),true);
395     }
396     catch(...) {
397     if (checkDCARef_)
398     throw edm::Exception(edm::errors::Configuration, "FillerDCASig:FillDataBlock()\n")
399     << "Error! Tau unmapped collection " << edmTauName_ << endl;
400 pharris 1.1 }
401     }
402     }
403     DCASigs_->Trim();
404     }
405 paus 1.2
406 pharris 1.1 //--------------------------------------------------------------------------------------------------
407     //Get DCA Significane
408 paus 1.2 //--------------------------------------------------------------------------------------------------
409     void FillerDCASig::calculateDCA(Double_t &iDCA3D ,Double_t &iDCA3DE ,Double_t &iDCA2D ,
410     Double_t &iDCA2DE,
411     Double_t &iDCARPhi3D,Double_t &iDCARPhi3DE,Double_t &iDCARPhi2D,
412     Double_t &iDCARPhi2DE,
413     const reco::Track *iTrack1,const reco::Track *iTrack2,
414     mithep::DCASig::EDCAType iType) {
415 pharris 1.1 iDCA3D = -1.0;
416     iDCA3DE = -1.0;
417     iDCA2D = -1.0;
418     iDCA2DE = -1.0;
419     iDCARPhi3D = -1.0;
420     iDCARPhi3DE = -1.0;
421     iDCARPhi2D = -1.0;
422     iDCARPhi2DE = -1.0;
423 paus 1.2
424     if (iTrack1==NULL)
425     return;
426     if (iTrack2==NULL)
427     return;
428 pharris 1.1
429     ParticleMass pion_mass = 0.139657;
430 paus 1.2 float pion_sigma = pion_mass*1.e-6;
431 pharris 1.1 ParticleMass muon_mass = 0.105658;
432 paus 1.2 float muon_sigma = muon_mass*1.e-6;
433 pharris 1.1 ParticleMass elec_mass = 0.000511;
434 paus 1.2 float elec_sigma = elec_mass*1.e-6;
435     float mass1 = elec_mass;
436     float mass2 = muon_mass;
437     float mass_sigma1 = elec_sigma;
438     float mass_sigma2 = muon_sigma;
439    
440     if (iType == DCASig::eETau ) mass1 = pion_sigma;
441     if (iType == DCASig::eEE ) mass2 = elec_sigma;
442     if (iType == DCASig::eMuTau ) { mass1 = muon_sigma; mass2 = pion_sigma; }
443     if (iType == DCASig::eMuMu ) { mass1 = muon_sigma; mass2 = muon_sigma; }
444     if (iType == DCASig::eTauTau) { mass1 = pion_sigma; mass2 = pion_sigma; }
445     if (iType == DCASig::eETau ) mass_sigma2 = pion_sigma;
446     if (iType == DCASig::eEE ) mass_sigma2 = elec_sigma;
447     if (iType == DCASig::eMuTau ) { mass_sigma1 = muon_sigma; mass_sigma2 = pion_sigma; }
448     if (iType == DCASig::eMuMu ) { mass_sigma1 = muon_sigma; mass_sigma2 = muon_sigma; }
449     if (iType == DCASig::eTauTau) { mass_sigma1 = pion_sigma; mass_sigma2 = pion_sigma; }
450 pharris 1.1
451     reco::TransientTrack transientTrack1;
452     reco::TransientTrack transientTrack2;
453 paus 1.2
454 pharris 1.1 transientTrack1 = transientTrackBuilder_->build(*iTrack1);
455     transientTrack2 = transientTrackBuilder_->build(*iTrack2);
456     reco::TransientTrack * trk1Ptr = & transientTrack1;
457     reco::TransientTrack * trk2Ptr = & transientTrack2;
458    
459     FreeTrajectoryState track1State = trk1Ptr->impactPointTSCP().theState();
460     FreeTrajectoryState track2State = trk2Ptr->impactPointTSCP().theState();
461 paus 1.2
462 pharris 1.1 if (trk1Ptr->impactPointTSCP().isValid() && trk2Ptr->impactPointTSCP().isValid()) {
463 paus 1.2
464 pharris 1.1 ClosestApproachInRPhi cApp;
465     TwoTrackMinimumDistance minDist;
466 paus 1.2
467 pharris 1.1 typedef ROOT::Math::SVector<double, 3> SVector3;
468 paus 1.2 typedef ROOT::Math::SMatrix<double, 3, 3, ROOT::Math::MatRepSym<double, 3> > SMatrixSym3D;
469 pharris 1.1 cApp.calculate(track1State,track2State);
470     minDist.calculate(track1State,track2State);
471     if (minDist.status()) {
472 paus 1.2
473 pharris 1.1 minDist.distance();
474     std::pair<GlobalPoint,GlobalPoint> pcaLeptons = minDist.points();
475     GlobalPoint track1PCA = pcaLeptons.first;
476     GlobalPoint track2PCA = pcaLeptons.second;
477 paus 1.2
478     //create a KinematicParticleFactory
479 pharris 1.1 KinematicParticleFactoryFromTransientTrack pFactory;
480 paus 1.2
481 pharris 1.1 //initial chi2 and ndf before kinematic fits.
482     float chi = 0.;
483     float ndf = 0.;
484 paus 1.2 RefCountedKinematicParticle track1Part = pFactory.particle(transientTrack1,mass1,chi,ndf,mass_sigma1);
485     RefCountedKinematicParticle track2Part = pFactory.particle(transientTrack2,mass2,chi,ndf,mass_sigma2);
486 pharris 1.1
487     SVector3 distanceVector(track1PCA.x()-track2PCA.x(),
488 paus 1.2 track1PCA.y()-track2PCA.y(),
489     track1PCA.z()-track2PCA.z());
490    
491 pharris 1.1 iDCA3D = ROOT::Math::Mag(distanceVector);
492 paus 1.2
493 pharris 1.1 std::vector<float> vvv(6);
494 paus 1.2
495 pharris 1.1 vvv[0] = track1Part->stateAtPoint(track1PCA).kinematicParametersError().matrix()(0,0);
496 paus 1.2 vvv[1] = track1Part->stateAtPoint(track1PCA).kinematicParametersError().matrix()(0,1);
497 pharris 1.1 vvv[2] = track1Part->stateAtPoint(track1PCA).kinematicParametersError().matrix()(1,1);
498     vvv[3] = track1Part->stateAtPoint(track1PCA).kinematicParametersError().matrix()(0,2);
499 paus 1.2 vvv[4] = track1Part->stateAtPoint(track1PCA).kinematicParametersError().matrix()(1,2);
500 pharris 1.1 vvv[5] = track1Part->stateAtPoint(track1PCA).kinematicParametersError().matrix()(2,2);
501 paus 1.2
502 pharris 1.1 SMatrixSym3D track1PCACov(vvv.begin(),vvv.end());
503 paus 1.2
504 pharris 1.1 vvv[0] = track2Part->stateAtPoint(track2PCA).kinematicParametersError().matrix()(0,0);
505 paus 1.2 vvv[1] = track2Part->stateAtPoint(track2PCA).kinematicParametersError().matrix()(0,1);
506 pharris 1.1 vvv[2] = track2Part->stateAtPoint(track2PCA).kinematicParametersError().matrix()(1,1);
507     vvv[3] = track2Part->stateAtPoint(track2PCA).kinematicParametersError().matrix()(0,2);
508 paus 1.2 vvv[4] = track2Part->stateAtPoint(track2PCA).kinematicParametersError().matrix()(1,2);
509 pharris 1.1 vvv[5] = track2Part->stateAtPoint(track2PCA).kinematicParametersError().matrix()(2,2);
510 paus 1.2
511 pharris 1.1 SMatrixSym3D track2PCACov(vvv.begin(),vvv.end());
512 paus 1.2
513 pharris 1.1 SMatrixSym3D totCov = track1PCACov + track2PCACov;
514 paus 1.2
515     if (iDCA3D != 0) iDCA3DE = sqrt(ROOT::Math::Similarity(totCov, distanceVector))/iDCA3D;
516     if (iDCA3D == 0) iDCA3DE = 0.;
517 pharris 1.1
518     distanceVector(2) = 0.0;
519     iDCA2D = ROOT::Math::Mag(distanceVector);
520 paus 1.2 if (iDCA2D != 0) iDCA2DE = sqrt(ROOT::Math::Similarity(totCov, distanceVector))/iDCA2D;
521     if (iDCA2D == 0) iDCA2DE = 0.;
522 pharris 1.1 }
523     if (cApp.status()) {
524 paus 1.2
525 pharris 1.1 cApp.distance();
526     std::pair<GlobalPoint,GlobalPoint> pcaLeptons = cApp.points();
527     GlobalPoint track1PCA = pcaLeptons.first;
528     GlobalPoint track2PCA = pcaLeptons.second;
529    
530     //Creating a KinematicParticleFactory
531     KinematicParticleFactoryFromTransientTrack pFactory;
532 paus 1.2
533 pharris 1.1 //initial chi2 and ndf before kinematic fits.
534     float chi = 0.;
535     float ndf = 0.;
536 paus 1.2 RefCountedKinematicParticle track1Part = pFactory.particle(transientTrack1,mass1,chi,ndf,mass_sigma1);
537     RefCountedKinematicParticle track2Part = pFactory.particle(transientTrack2,mass2,chi,ndf,mass_sigma2);
538    
539 pharris 1.1 SVector3 distanceVector(track1PCA.x()-track2PCA.x(),
540 paus 1.2 track1PCA.y()-track2PCA.y(),
541     track1PCA.z()-track2PCA.z());
542 pharris 1.1 iDCARPhi3D = ROOT::Math::Mag(distanceVector);
543 paus 1.2
544 pharris 1.1 std::vector<float> vvv(6);
545 paus 1.2
546 pharris 1.1 vvv[0] = track1Part->stateAtPoint(track1PCA).kinematicParametersError().matrix()(0,0);
547 paus 1.2 vvv[1] = track1Part->stateAtPoint(track1PCA).kinematicParametersError().matrix()(0,1);
548 pharris 1.1 vvv[2] = track1Part->stateAtPoint(track1PCA).kinematicParametersError().matrix()(1,1);
549     vvv[3] = track1Part->stateAtPoint(track1PCA).kinematicParametersError().matrix()(0,2);
550 paus 1.2 vvv[4] = track1Part->stateAtPoint(track1PCA).kinematicParametersError().matrix()(1,2);
551 pharris 1.1 vvv[5] = track1Part->stateAtPoint(track1PCA).kinematicParametersError().matrix()(2,2);
552 paus 1.2
553 pharris 1.1 SMatrixSym3D track1PCACov(vvv.begin(),vvv.end());
554    
555     vvv[0] = track2Part->stateAtPoint(track2PCA).kinematicParametersError().matrix()(0,0);
556 paus 1.2 vvv[1] = track2Part->stateAtPoint(track2PCA).kinematicParametersError().matrix()(0,1);
557 pharris 1.1 vvv[2] = track2Part->stateAtPoint(track2PCA).kinematicParametersError().matrix()(1,1);
558     vvv[3] = track2Part->stateAtPoint(track2PCA).kinematicParametersError().matrix()(0,2);
559 paus 1.2 vvv[4] = track2Part->stateAtPoint(track2PCA).kinematicParametersError().matrix()(1,2);
560 pharris 1.1 vvv[5] = track2Part->stateAtPoint(track2PCA).kinematicParametersError().matrix()(2,2);
561 paus 1.2
562 pharris 1.1 SMatrixSym3D track2PCACov(vvv.begin(),vvv.end());
563    
564     SMatrixSym3D totCov = track1PCACov + track2PCACov;
565 paus 1.2
566     if (iDCARPhi3D != 0)
567     iDCARPhi3DE = sqrt(ROOT::Math::Similarity(totCov, distanceVector))/iDCARPhi3D;
568     if (iDCARPhi3D == 0)
569     iDCARPhi3DE = 0.;
570 pharris 1.1
571     distanceVector(2) = 0.0;
572     iDCARPhi2D = ROOT::Math::Mag(distanceVector);
573 paus 1.2 if (iDCARPhi2D != 0)
574     iDCARPhi2DE = sqrt(ROOT::Math::Similarity(totCov, distanceVector))/iDCARPhi2D;
575     if (iDCARPhi2D == 0)
576     iDCARPhi2DE = 0.;
577 pharris 1.1 }
578     }
579     }
580