ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitEdm/Producers/src/ProducerConversions.cc
Revision: 1.21
Committed: Fri Dec 11 17:46:24 2009 UTC (15 years, 5 months ago) by bendavid
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_012f, Mit_012e
Changes since 1.20: +4 -13 lines
Log Message:
Remove track quality cuts, which do not really belong here

File Contents

# User Rev Content
1 bendavid 1.21 // $Id: ProducerConversions.cc,v 1.20 2009/12/08 17:40:04 bendavid Exp $
2 bendavid 1.1
3 loizides 1.5 #include "MitEdm/Producers/interface/ProducerConversions.h"
4 bendavid 1.1 #include "DataFormats/Common/interface/Handle.h"
5     #include "DataFormats/TrackReco/interface/Track.h"
6     #include "DataFormats/TrackReco/interface/TrackFwd.h"
7 bendavid 1.6 #include "DataFormats/VertexReco/interface/Vertex.h"
8     #include "DataFormats/VertexReco/interface/VertexFwd.h"
9 bendavid 1.18 #include "MagneticField/Engine/interface/MagneticField.h"
10     #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
11 bendavid 1.19 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
12     #include "TrackingTools/PatternTools/interface/ClosestApproachInRPhi.h"
13 bendavid 1.7 #include "MitEdm/Producers/interface/HitDropperRecord.h"
14     #include "MitEdm/Producers/interface/HitDropper.h"
15 bendavid 1.1 #include "MitEdm/DataFormats/interface/Types.h"
16 loizides 1.5 #include "MitEdm/DataFormats/interface/Collections.h"
17 bendavid 1.1 #include "MitEdm/DataFormats/interface/DecayPart.h"
18     #include "MitEdm/DataFormats/interface/StablePart.h"
19 bendavid 1.6 #include "MitEdm/DataFormats/interface/StableData.h"
20 bendavid 1.1 #include "MitEdm/VertexFitInterface/interface/MvfInterface.h"
21 bendavid 1.20 #include <TMath.h>
22 bendavid 1.1
23     using namespace std;
24     using namespace edm;
25     using namespace reco;
26     using namespace mitedm;
27     using namespace mithep;
28    
29     //--------------------------------------------------------------------------------------------------
30     ProducerConversions::ProducerConversions(const ParameterSet& cfg) :
31 paus 1.11 BaseCandProducer (cfg),
32     iStables1_ (cfg.getUntrackedParameter<string>("iStables1","")),
33     iStables2_ (cfg.getUntrackedParameter<string>("iStables2","")),
34     iPVertexes_ (cfg.getUntrackedParameter<string>("iPVertexes","offlinePrimaryVerticesWithBS")),
35     usePVertex_ (cfg.getUntrackedParameter<bool> ("usePVertex",true)),
36     convConstraint_ (cfg.getUntrackedParameter<bool> ("convConstraint",false)),
37     convConstraint3D_(cfg.getUntrackedParameter<bool> ("convConstraint3D",true)),
38 bendavid 1.15 rhoMin_ (cfg.getUntrackedParameter<double>("rhoMin",0.0)),
39 bendavid 1.20 useHitDropper_ (cfg.getUntrackedParameter<bool> ("useHitDropper",true)),
40 bendavid 1.21 applyChargeConstraint_(cfg.getUntrackedParameter<bool> ("applyChargeConstraint",false))
41 bendavid 1.1 {
42 loizides 1.5 // Constructor.
43    
44 bendavid 1.1 produces<DecayPartCol>();
45     }
46    
47     //--------------------------------------------------------------------------------------------------
48     ProducerConversions::~ProducerConversions()
49     {
50 loizides 1.5 // Destructor.
51 bendavid 1.1 }
52    
53     //--------------------------------------------------------------------------------------------------
54     void ProducerConversions::produce(Event &evt, const EventSetup &setup)
55     {
56 loizides 1.5 // Produce our DecayPartCol.
57    
58 bendavid 1.1 // First input collection
59     Handle<StablePartCol> hStables1;
60 bendavid 1.9 if (!GetProduct(iStables1_, hStables1, evt)) {
61     printf("Stable collection 1 not found in ProducerConversions\n");
62 bendavid 1.1 return;
63 bendavid 1.9 }
64 bendavid 1.1 const StablePartCol *pS1 = hStables1.product();
65     // Second input collection
66     Handle<StablePartCol> hStables2;
67 bendavid 1.9 if (!GetProduct(iStables2_, hStables2, evt)) {
68     printf("Stable collection 2 not found in ProducerConversions\n");
69 bendavid 1.1 return;
70 bendavid 1.9 }
71 bendavid 1.1 const StablePartCol *pS2 = hStables2.product();
72    
73 bendavid 1.6 const reco::Vertex *vertex = 0;
74     mitedm::VertexPtr vPtr;
75     if (usePVertex_) {
76     // Primary vertex collection
77     Handle<reco::VertexCollection> hVertexes;
78     if (!GetProduct(iPVertexes_, hVertexes, evt))
79     return;
80     const reco::VertexCollection *pV = hVertexes.product();
81    
82     //Choose the primary vertex with the largest number of tracks
83     UInt_t maxTracks=0;
84     for (UInt_t i=0; i<pV->size(); ++i) {
85     const reco::Vertex &v = pV->at(i);
86     UInt_t nTracks = v.tracksSize();
87     if (nTracks >= maxTracks) {
88     maxTracks = nTracks;
89     vertex = &v;
90     vPtr = mitedm::VertexPtr(hVertexes,i);
91     }
92     }
93     }
94 bendavid 1.7
95 paus 1.11 // Get hit dropper
96 bendavid 1.7 ESHandle<HitDropper> hDropper;
97     setup.get<HitDropperRecord>().get("HitDropper",hDropper);
98     const HitDropper *dropper = hDropper.product();
99 bendavid 1.6
100 bendavid 1.18 //Get Magnetic Field from event setup, taking value at (0,0,0)
101     edm::ESHandle<MagneticField> magneticField;
102     setup.get<IdealMagneticFieldRecord>().get(magneticField);
103     const double bfield = magneticField->inTesla(GlobalPoint(0.,0.,0.)).z();
104    
105 bendavid 1.1 // Create the output collection
106     auto_ptr<DecayPartCol> pD(new DecayPartCol());
107 bendavid 1.6
108 bendavid 1.19 ClosestApproachInRPhi helixIntersector;
109    
110 bendavid 1.20 int nFits = 0;
111    
112     //printf("S1 size = %i\n", pS1->size());
113    
114 bendavid 1.1 // Simple double loop
115     for (UInt_t i = 0; i<pS1->size(); ++i) {
116     const StablePart &s1 = pS1->at(i);
117    
118 bendavid 1.20 const reco::Track * t1 = s1.track();
119 bendavid 1.21
120 bendavid 1.1 UInt_t j;
121     if (iStables1_ == iStables2_)
122     j = i+1;
123     else
124     j = 0;
125    
126 bendavid 1.19 TrajectoryStateTransform tsTransform;
127    
128     FreeTrajectoryState initialState1 = tsTransform.initialFreeState(*s1.track(),&*magneticField);
129    
130 bendavid 1.1 for (; j<pS2->size(); ++j) {
131     const StablePart &s2 = pS2->at(j);
132    
133 bendavid 1.21 //Do fast helix fit to check if there's any hope
134 bendavid 1.13 const reco::Track * t2 = s2.track();
135    
136     int trackCharge = t1->charge() + t2->charge();
137 bendavid 1.19
138 bendavid 1.13 double dR0 = 0.0;
139    
140 bendavid 1.20 if (!applyChargeConstraint_ || trackCharge==0) {
141 bendavid 1.19 FreeTrajectoryState initialState2 = tsTransform.initialFreeState(*s2.track(),&*magneticField);
142     helixIntersector.calculate(initialState1, initialState2);
143     if (helixIntersector.status())
144     dR0 = helixIntersector.crossingPoint().perp();
145 bendavid 1.3 }
146 bendavid 1.13
147     int fitStatus = 0;
148     MultiVertexFitterD fit;
149 bendavid 1.20 if ( (!applyChargeConstraint_ || trackCharge==0) && dR0 > rhoMin_) {
150 bendavid 1.13
151     // Vertex fit now, possibly with conversion constraint
152 bendavid 1.20 nFits++;
153 bendavid 1.13
154 bendavid 1.18 fit.init(bfield); // Reset to the magnetic field from the event setup
155 bendavid 1.13 MvfInterface fitInt(&fit);
156     fitInt.addTrack(s1.track(),1,s1.mass(),MultiVertexFitterD::VERTEX_1);
157     fitInt.addTrack(s2.track(),2,s2.mass(),MultiVertexFitterD::VERTEX_1);
158     if (convConstraint3D_) {
159     fit.conversion_3d(MultiVertexFitterD::VERTEX_1);
160     //printf("applying 3d conversion constraint\n");
161     }
162     else if (convConstraint_) {
163     fit.conversion_2d(MultiVertexFitterD::VERTEX_1);
164     //printf("applying 2d conversion constraint\n");
165     }
166     //initialize primary vertex parameters in the fitter
167     if (usePVertex_) {
168     float vErr[3][3];
169     for (UInt_t vi=0; vi<3; ++vi)
170     for (UInt_t vj=0; vj<3; ++vj)
171     vErr[vi][vj] = vertex->covariance(vi,vj);
172    
173     fit.setPrimaryVertex(vertex->x(),vertex->y(),vertex->z());
174     fit.setPrimaryVertexError(vErr);
175     }
176    
177     fitStatus = fit.fit();
178 bendavid 1.6 }
179    
180     if (fitStatus) {
181 bendavid 1.3 DecayPart *d = new DecayPart(oPid_,DecayPart::Fast);
182 bendavid 1.1
183 bendavid 1.3 // Update temporarily some of the quantities (prob, chi2, nDoF, mass, lxy, pt, fourMomentum)
184     d->setProb(fit.prob());
185     d->setChi2(fit.chisq());
186     d->setNdof(fit.ndof());
187    
188     FourVector p4Fitted(0.,0.,0.,0.);
189     p4Fitted += fit.getTrackP4(1);
190     p4Fitted += fit.getTrackP4(2);
191     d->setFourMomentum(p4Fitted);
192 paus 1.12 d->setPosition(fit.getVertex (MultiVertexFitterD::VERTEX_1));
193     d->setError (fit.getErrorMatrix(MultiVertexFitterD::VERTEX_1));
194 bendavid 1.3 float mass, massErr;
195     const int trksIds[2] = { 1, 2 };
196     mass = fit.getMass(2,trksIds,massErr);
197    
198 bendavid 1.6 ThreeVector p3Fitted(p4Fitted.px(), p4Fitted.py(), p4Fitted.pz());
199    
200 paus 1.11 // Get decay length in xy plane
201 bendavid 1.6 float dl, dlErr;
202 paus 1.12 dl = fit.getDecayLength(MultiVertexFitterD::PRIMARY_VERTEX, MultiVertexFitterD::VERTEX_1,
203 bendavid 1.6 p3Fitted, dlErr);
204    
205 paus 1.11 // Get Z decay length
206 bendavid 1.6 float dlz, dlzErr;
207 paus 1.12 dlz = fit.getZDecayLength(MultiVertexFitterD::PRIMARY_VERTEX, MultiVertexFitterD::VERTEX_1,
208 bendavid 1.6 p3Fitted, dlzErr);
209    
210 paus 1.11 // Get impact parameter
211 bendavid 1.6 float dxy, dxyErr;
212 paus 1.12 dxy = fit.getImpactPar(MultiVertexFitterD::PRIMARY_VERTEX, MultiVertexFitterD::VERTEX_1,
213 bendavid 1.6 p3Fitted, dxyErr);
214 bendavid 1.1
215 bendavid 1.8 BasePartPtr ptr1(hStables1,i);
216     BasePartPtr ptr2(hStables2,j);
217    
218     StableData c1(fit.getTrackP4(1).px(),fit.getTrackP4(1).py(), fit.getTrackP4(1).pz(), ptr1);
219     StableData c2(fit.getTrackP4(2).px(),fit.getTrackP4(2).py(), fit.getTrackP4(2).pz(), ptr2);
220    
221 paus 1.12 const ThreeVector vtxPos = fit.getVertex(MultiVertexFitterD::VERTEX_1);
222 loizides 1.14 const ThreeVector trkMom1(fit.getTrackP4(1).px(),
223     fit.getTrackP4(1).py(),
224     fit.getTrackP4(1).pz());
225     const ThreeVector trkMom2(fit.getTrackP4(2).px(),
226     fit.getTrackP4(2).py(),
227     fit.getTrackP4(2).pz());
228 bendavid 1.8
229 paus 1.11 // Build corrected HitPattern for StableData, removing hits before the fit vertex
230 bendavid 1.15 if (useHitDropper_) {
231 bendavid 1.17 reco::HitPattern hits1 = dropper->CorrectedHitsAOD(s1.track(), vtxPos, trkMom1, dlErr, dlzErr);
232     reco::HitPattern hits2 = dropper->CorrectedHitsAOD(s2.track(), vtxPos, trkMom2, dlErr, dlzErr);
233    
234 bendavid 1.15 c1.SetHits(hits1);
235     c2.SetHits(hits2);
236     c1.SetHitsFilled();
237     c2.SetHitsFilled();
238     }
239 bendavid 1.8
240     d->addStableChild(c1);
241     d->addStableChild(c2);
242    
243    
244 bendavid 1.6 d->setFittedMass (mass);
245     d->setFittedMassError(massErr);
246    
247     d->setLxy(dl);
248     d->setLxyError(dlErr);
249     d->setLxyToPv(dl);
250     d->setLxyToPvError(dlErr);
251    
252     d->setLz(dlz);
253     d->setLzError(dlzErr);
254     d->setLzToPv(dlz);
255     d->setLzToPvError(dlzErr);
256    
257     d->setDxy(dxy);
258     d->setDxyError(dxyErr);
259     d->setDxyToPv(dxy);
260     d->setDxyToPvError(dxyErr);
261    
262     if (usePVertex_)
263     d->setPrimaryVertex(vPtr);
264    
265     // Put the result into our collection
266 bendavid 1.13 pD->push_back(*d);
267 loizides 1.5
268 bendavid 1.1 delete d;
269     }
270     }
271     }
272 bendavid 1.20
273     //printf("nConversionFits = %i\n",nFits);
274 bendavid 1.1
275     // Write the collection even if it is empty
276 loizides 1.5 if (0) {
277     cout << " ProducerConversions::produce - " << pD->size() << " entries collection created -"
278     << " (Pid: " << oPid_ << ")\n";
279     }
280 bendavid 1.1 evt.put(pD);
281     }
282    
283     //define this as a plug-in
284     DEFINE_FWK_MODULE(ProducerConversions);