Producers/src/ProducerConversions.cc

// $Id: ProducerConversions.cc,v 1.6 2008/09/30 12:59:14 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"

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))
    return;  
  const StablePartCol *pS1 = hStables1.product();
  // Second input collection
  Handle<StablePartCol> hStables2;
  if (!GetProduct(iStables2_, hStables2, evt))
    return;
  const StablePartCol *pS2 = hStables2.product();

  ThreeVector pvPos;
  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);
      }
    }
    pvPos.SetXYZ(vertex->x(),
                 vertex->y(),
                 vertex->z());
  }
  else
    pvPos.SetXYZ(0.0,0.0,0.0);
  
  //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);

      // Vertex fit now, possibly with conversion constraint
      MultiVertexFitter fit;
      fit.init(3.8); // Reset to the MC magnetic field of 3.8 Tesla
      MvfInterface fitInt(&fit);
      fitInt.addTrack(s1.track(),1,s1.mass(),MultiVertexFitter::VERTEX_1);
      fitInt.addTrack(s2.track(),2,s2.mass(),MultiVertexFitter::VERTEX_1);
      if (convConstraint3D_) {
        fit.conversion_3d(MultiVertexFitter::VERTEX_1);
        //printf("applying 3d conversion constraint\n");
      }
      else if (convConstraint_) {
        fit.conversion_2d(MultiVertexFitter::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);
      }
      
      //only perform fit for oppositely-charged tracks
      int trackCharge = s1.track()->charge() + s2.track()->charge();
      int fitStatus = 0;
      if (trackCharge==0)
        fitStatus = fit.fit();
        
      if (fitStatus) {
        DecayPart *d = new DecayPart(oPid_,DecayPart::Fast);
        
        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(MultiVertexFitter::VERTEX_1);
        //build corrected HitPattern for StableData, removing hits before the fit vertex
        reco::HitPattern hits1 = dropper->CorrectedHits(s1.track(), vtxPos, pvPos);
        reco::HitPattern hits2 = dropper->CorrectedHits(s2.track(), vtxPos, pvPos);
        
        c1.SetHits(hits1);
        c2.SetHits(hits2);
        
        d->addStableChild(c1);
        d->addStableChild(c2);
        
        // 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     (MultiVertexFitter::VERTEX_1));
        d->setError   (fit.getErrorMatrix(MultiVertexFitter::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(MultiVertexFitter::PRIMARY_VERTEX, MultiVertexFitter::VERTEX_1,
               p3Fitted, dlErr);
               
        //Get Z decay length               
        float dlz, dlzErr;
        dlz = fit.getZDecayLength(MultiVertexFitter::PRIMARY_VERTEX, MultiVertexFitter::VERTEX_1,
               p3Fitted, dlzErr);
               
        //get impact parameter               
        float dxy, dxyErr;
        dxy = fit.getImpactPar(MultiVertexFitter::PRIMARY_VERTEX, MultiVertexFitter::VERTEX_1,
               p3Fitted, dxyErr);

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