Producers/src/ProducerConversions.cc

// $Id: ProducerConversions.cc,v 1.13 2008/11/14 16:25:43 bendavid Exp $

#include "MitEdm/Producers/interface/ProducerConversions.h"
#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/VertexReco/interface/VertexFwd.h"
#include "MitEdm/Producers/interface/HitDropperRecord.h"
#include "MitEdm/Producers/interface/HitDropper.h"
#include "MitEdm/DataFormats/interface/Types.h"
#include "MitEdm/DataFormats/interface/Collections.h"
#include "MitEdm/DataFormats/interface/DecayPart.h"
#include "MitEdm/DataFormats/interface/StablePart.h"
#include "MitEdm/DataFormats/interface/StableData.h"
#include "MitEdm/VertexFitInterface/interface/MvfInterface.h"
#include "MitEdm/VertexFitInterface/interface/HisInterface.h"

using namespace std;
using namespace edm;
using namespace reco;
using namespace mitedm;
using namespace mithep;

//--------------------------------------------------------------------------------------------------
ProducerConversions::ProducerConversions(const ParameterSet& cfg) :
  BaseCandProducer (cfg),
  iStables1_       (cfg.getUntrackedParameter<string>("iStables1","")),
  iStables2_       (cfg.getUntrackedParameter<string>("iStables2","")),
  iPVertexes_      (cfg.getUntrackedParameter<string>("iPVertexes","offlinePrimaryVerticesWithBS")),
  usePVertex_      (cfg.getUntrackedParameter<bool>  ("usePVertex",true)),
  convConstraint_  (cfg.getUntrackedParameter<bool>  ("convConstraint",false)),
  convConstraint3D_(cfg.getUntrackedParameter<bool>  ("convConstraint3D",true)),
  rhoMin_          (cfg.getUntrackedParameter<double>("rhoMin",0.0))
{
  // Constructor.

  produces<DecayPartCol>();
}

//--------------------------------------------------------------------------------------------------
ProducerConversions::~ProducerConversions()
{
  // Destructor.
}

//--------------------------------------------------------------------------------------------------
void ProducerConversions::produce(Event &evt, const EventSetup &setup)
{
  // Produce our DecayPartCol.

  // First input collection
  Handle<StablePartCol> hStables1;
  if (!GetProduct(iStables1_, hStables1, evt)) {
    printf("Stable collection 1 not found in ProducerConversions\n");
    return;  
  }
  const StablePartCol *pS1 = hStables1.product();
  // Second input collection
  Handle<StablePartCol> hStables2;
  if (!GetProduct(iStables2_, hStables2, evt)) {
    printf("Stable collection 2 not found in ProducerConversions\n");
    return;
  }
  const StablePartCol *pS2 = hStables2.product();

  const reco::Vertex *vertex = 0;
  mitedm::VertexPtr vPtr;
  if (usePVertex_) {
    // Primary vertex collection
    Handle<reco::VertexCollection> hVertexes;
    if (!GetProduct(iPVertexes_, hVertexes, evt))
      return;
    const reco::VertexCollection *pV = hVertexes.product();
  
    //Choose the primary vertex with the largest number of tracks
    UInt_t maxTracks=0;
    for (UInt_t i=0; i<pV->size(); ++i) {
      const reco::Vertex &v = pV->at(i);
      UInt_t nTracks = v.tracksSize();
      if (nTracks >= maxTracks) {
        maxTracks = nTracks;
        vertex = &v;
        vPtr = mitedm::VertexPtr(hVertexes,i);
      }
    }
  }
  
  // Get hit dropper
  ESHandle<HitDropper> hDropper;
  setup.get<HitDropperRecord>().get("HitDropper",hDropper);
  const HitDropper *dropper = hDropper.product();
  
  // Create the output collection
  auto_ptr<DecayPartCol> pD(new DecayPartCol());
  
  // Simple double loop
  for (UInt_t i = 0; i<pS1->size(); ++i) {
    const StablePart &s1 =  pS1->at(i);
    
    UInt_t j;
    if (iStables1_ == iStables2_)
      j = i+1; 
    else
      j = 0;
    
    for (; j<pS2->size(); ++j) {
      const StablePart &s2 = pS2->at(j);

            //Do fast helix fit to check if there's any hope
      const reco::Track * t1 = s1.track();
      const reco::Track * t2 = s2.track();
      
      int trackCharge = t1->charge() + t2->charge();
      double dR0 = 0.0;
      
      if (trackCharge==0) { 
        HisInterface            hisInt(t1,t2);
        const mithep::HelixIntersector *his = hisInt.hISector();
        dR0 = sqrt(his->IntersectionI(0).Location().X()*his->IntersectionI(0).Location().X()+
                  his->IntersectionI(0).Location().Y()*his->IntersectionI(0).Location().Y());
      }

                 
      int fitStatus = 0;
      MultiVertexFitterD fit;      
      if (trackCharge==0 && dR0 > rhoMin_) {
      
        // Vertex fit now, possibly with conversion constraint
  
        fit.init(3.8); // Reset to the MC magnetic field of 3.8 Tesla
        MvfInterface fitInt(&fit);
        fitInt.addTrack(s1.track(),1,s1.mass(),MultiVertexFitterD::VERTEX_1);
        fitInt.addTrack(s2.track(),2,s2.mass(),MultiVertexFitterD::VERTEX_1);
        if (convConstraint3D_) {
          fit.conversion_3d(MultiVertexFitterD::VERTEX_1);
          //printf("applying 3d conversion constraint\n");
        }
        else if (convConstraint_) {
          fit.conversion_2d(MultiVertexFitterD::VERTEX_1);
          //printf("applying 2d conversion constraint\n");
        }
        //initialize primary vertex parameters in the fitter
        if (usePVertex_) {
          float vErr[3][3];
          for (UInt_t vi=0; vi<3; ++vi)
            for (UInt_t vj=0; vj<3; ++vj)
              vErr[vi][vj] = vertex->covariance(vi,vj);
              
          fit.setPrimaryVertex(vertex->x(),vertex->y(),vertex->z());
          fit.setPrimaryVertexError(vErr);
        }
        
        fitStatus = fit.fit();     
      }
        
      if (fitStatus) {
        DecayPart *d = new DecayPart(oPid_,DecayPart::Fast);
        
        // Update temporarily some of the quantities (prob, chi2, nDoF, mass, lxy, pt, fourMomentum)
        d->setProb(fit.prob());
        d->setChi2(fit.chisq());
        d->setNdof(fit.ndof());

        FourVector p4Fitted(0.,0.,0.,0.);
        p4Fitted += fit.getTrackP4(1);
        p4Fitted += fit.getTrackP4(2);
        d->setFourMomentum(p4Fitted);
        d->setPosition(fit.getVertex     (MultiVertexFitterD::VERTEX_1));
        d->setError   (fit.getErrorMatrix(MultiVertexFitterD::VERTEX_1));
        float mass, massErr;
        const int trksIds[2] = { 1, 2 };
        mass = fit.getMass(2,trksIds,massErr);
        
        ThreeVector p3Fitted(p4Fitted.px(), p4Fitted.py(), p4Fitted.pz());
        
        // Get decay length in xy plane
        float dl, dlErr;
        dl = fit.getDecayLength(MultiVertexFitterD::PRIMARY_VERTEX, MultiVertexFitterD::VERTEX_1,
               p3Fitted, dlErr);
               
        // Get Z decay length               
        float dlz, dlzErr;
        dlz = fit.getZDecayLength(MultiVertexFitterD::PRIMARY_VERTEX, MultiVertexFitterD::VERTEX_1,
               p3Fitted, dlzErr);
               
        // Get impact parameter               
        float dxy, dxyErr;
        dxy = fit.getImpactPar(MultiVertexFitterD::PRIMARY_VERTEX, MultiVertexFitterD::VERTEX_1,
               p3Fitted, dxyErr);

        BasePartPtr ptr1(hStables1,i);
        BasePartPtr ptr2(hStables2,j);       
        
        StableData c1(fit.getTrackP4(1).px(),fit.getTrackP4(1).py(), fit.getTrackP4(1).pz(), ptr1);
        StableData c2(fit.getTrackP4(2).px(),fit.getTrackP4(2).py(), fit.getTrackP4(2).pz(), ptr2);
        
        const ThreeVector vtxPos = fit.getVertex(MultiVertexFitterD::VERTEX_1);
        const ThreeVector trkMom1(fit.getTrackP4(1).px(),
                                  fit.getTrackP4(1).py(), 
                                  fit.getTrackP4(1).pz());
        const ThreeVector trkMom2(fit.getTrackP4(2).px(),
                                  fit.getTrackP4(2).py(), 
                                  fit.getTrackP4(2).pz());
        
        // Build corrected HitPattern for StableData, removing hits before the fit vertex
        reco::HitPattern hits1 = dropper->CorrectedHits(s1.track(), vtxPos, trkMom1, dlErr, dlzErr);
        reco::HitPattern hits2 = dropper->CorrectedHits(s2.track(), vtxPos, trkMom2, dlErr, dlzErr);
        
        c1.SetHits(hits1);
        c2.SetHits(hits2);
        c1.SetHitsFilled();
        c2.SetHitsFilled();
        
        d->addStableChild(c1);
        d->addStableChild(c2);
               
               
        d->setFittedMass     (mass);
        d->setFittedMassError(massErr);
        
        d->setLxy(dl);
        d->setLxyError(dlErr);
        d->setLxyToPv(dl);
        d->setLxyToPvError(dlErr);
        
        d->setLz(dlz);
        d->setLzError(dlzErr);
        d->setLzToPv(dlz);
        d->setLzToPvError(dlzErr);
        
        d->setDxy(dxy);
        d->setDxyError(dxyErr);
        d->setDxyToPv(dxy);
        d->setDxyToPvError(dxyErr);
        
        if (usePVertex_)
          d->setPrimaryVertex(vPtr);
        
        // Put the result into our collection
        pD->push_back(*d);

        delete d;
      }
    }
  }

  // Write the collection even if it is empty
  if (0) {
    cout << " ProducerConversions::produce - " << pD->size() << " entries collection created -"
         << " (Pid: " << oPid_ << ")\n";
  }
  evt.put(pD);
}

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