Producers/src/ProducerD2SS.cc

// $Id: ProducerD2SS.cc,v 1.13 2008/11/13 17:08:31 paus Exp $

#include "MitEdm/Producers/interface/ProducerD2SS.h"
#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.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;

//--------------------------------------------------------------------------------------------------
ProducerD2SS::ProducerD2SS(const ParameterSet& cfg) :
  BaseCandProducer(cfg),
  iStables1_(cfg.getUntrackedParameter<string>("iStables1","")),
  iStables2_(cfg.getUntrackedParameter<string>("iStables2",""))
{
  // Constructor
  produces<DecayPartCol>();
}

//--------------------------------------------------------------------------------------------------
ProducerD2SS::~ProducerD2SS()
{
  // Destructor
}

//--------------------------------------------------------------------------------------------------
void ProducerD2SS::produce(Event &evt, const EventSetup &setup)
{
  // Produce the DecayPartCol for a topology: D -> S S

  // 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();

  // 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 (temporary screwing around)
      MultiVertexFitterD fit;
      fit.init(3.8); // Reset to the MC magnetic field of 4 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 (fit.fit()) {
        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);

        // Put the result into our collection
        pD->push_back(*d);
        delete d;
      }
    }
  }

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

// Define this as a plug-in
DEFINE_FWK_MODULE(ProducerD2SS);
Revision:	1.14
Committed:	Sun Oct 4 12:49:26 2009 UTC (15 years, 7 months ago) by bendavid
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, Mit_028, Mit_027a, Mit_027, Mit_026, Mit_025e, Mit_025d, Mit_025c, Mit_025b, Mit_025a, Mit_025, Mit_025pre2, Mit_024b, Mit_025pre1, Mit_024a, Mit_024, Mit_023, Mit_022a, Mit_022, Mit_020d, TMit_020d, Mit_020c, Mit_021, Mit_021pre2, Mit_021pre1, Mit_020b, Mit_020a, Mit_020, Mit_020pre1, Mit_018, Mit_017, Mit_017pre3, Mit_017pre2, Mit_017pre1, V07-05-00, Mit_016, Mit_015b, Mit_015a, Mit_015, Mit_014e, Mit_014d, Mit_014c, Mit_014b, ConvRejection-10-06-09, Mit_014a, Mit_014, Mit_014pre3, Mit_014pre2, Mit_014pre1, Mit_013d, Mit_013c, Mit_013b, Mit_013a, Mit_013, Mit_013pre1, Mit_012i, Mit_012h, Mit_012g, Mit_012f, Mit_012e, Mit_012d, Mit_012c, Mit_012b, Mit_012a, Mit_012, Mit_011a, HEAD
Branch point for:	Mit_025c_branch
Changes since 1.13:	+3 -1 lines
Log Message:	Make HitDropper work in AOD
#	User	Rev	Content
1	bendavid	1.14	// $Id: ProducerD2SS.cc,v 1.13 2008/11/13 17:08:31 paus Exp $
2	loizides	1.4
3	loizides	1.7	#include "MitEdm/Producers/interface/ProducerD2SS.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.11	#include "MitEdm/Producers/interface/HitDropperRecord.h"
8			#include "MitEdm/Producers/interface/HitDropper.h"
9	bendavid	1.1	#include "MitEdm/DataFormats/interface/Types.h"
10	loizides	1.7	#include "MitEdm/DataFormats/interface/Collections.h"
11	bendavid	1.1	#include "MitEdm/DataFormats/interface/DecayPart.h"
12			#include "MitEdm/DataFormats/interface/StablePart.h"
13	bendavid	1.8	#include "MitEdm/DataFormats/interface/StableData.h"
14	bendavid	1.1	#include "MitEdm/VertexFitInterface/interface/MvfInterface.h"
15
16			using namespace std;
17			using namespace edm;
18			using namespace reco;
19			using namespace mitedm;
20	loizides	1.4	using namespace mithep;
21	bendavid	1.1
22			//--------------------------------------------------------------------------------------------------
23			ProducerD2SS::ProducerD2SS(const ParameterSet& cfg) :
24	loizides	1.7	BaseCandProducer(cfg),
25			iStables1_(cfg.getUntrackedParameter<string>("iStables1","")),
26			iStables2_(cfg.getUntrackedParameter<string>("iStables2",""))
27	bendavid	1.1	{
28	paus	1.12	// Constructor
29	bendavid	1.1	produces<DecayPartCol>();
30			}
31
32			//--------------------------------------------------------------------------------------------------
33			ProducerD2SS::~ProducerD2SS()
34			{
35	paus	1.12	// Destructor
36	bendavid	1.1	}
37
38			//--------------------------------------------------------------------------------------------------
39			void ProducerD2SS::produce(Event &evt, const EventSetup &setup)
40			{
41	paus	1.12	// Produce the DecayPartCol for a topology: D -> S S
42	loizides	1.7
43	bendavid	1.1	// First input collection
44			Handle<StablePartCol> hStables1;
45			if (!GetProduct(iStables1_, hStables1, evt))
46	paus	1.13	return;
47	bendavid	1.1	const StablePartCol *pS1 = hStables1.product();
48			// Second input collection
49			Handle<StablePartCol> hStables2;
50			if (!GetProduct(iStables2_, hStables2, evt))
51			return;
52			const StablePartCol *pS2 = hStables2.product();
53
54	paus	1.12	// Get hit dropper
55	bendavid	1.11	ESHandle<HitDropper> hDropper;
56			setup.get<HitDropperRecord>().get("HitDropper",hDropper);
57			const HitDropper *dropper = hDropper.product();
58	paus	1.13
59	bendavid	1.1	// Create the output collection
60			auto_ptr<DecayPartCol> pD(new DecayPartCol());
61
62			// Simple double loop
63			for (UInt_t i = 0; i<pS1->size(); ++i) {
64			const StablePart &s1 = pS1->at(i);
65	paus	1.13
66	bendavid	1.1	UInt_t j;
67			if (iStables1_ == iStables2_)
68	paus	1.13	j = i+1;
69	bendavid	1.1	else
70			j = 0;
71	paus	1.13
72	bendavid	1.1	for (; j<pS2->size(); ++j) {
73			const StablePart &s2 = pS2->at(j);
74
75			// Vertex fit now (temporary screwing around)
76	paus	1.13	MultiVertexFitterD fit;
77	loizides	1.7	fit.init(3.8); // Reset to the MC magnetic field of 4 Tesla
78	bendavid	1.1	MvfInterface fitInt(&fit);
79	paus	1.13	fitInt.addTrack(s1.track(),1,s1.mass(),MultiVertexFitterD::VERTEX_1);
80			fitInt.addTrack(s2.track(),2,s2.mass(),MultiVertexFitterD::VERTEX_1);
81	bendavid	1.1	if (fit.fit()) {
82	paus	1.13	DecayPart *d = new DecayPart(oPid_,DecayPart::Fast);
83
84			// Update temporarily some of the quantities (prob, chi2, nDoF, mass, lxy, pt, fourMomentum)
85			d->setProb(fit.prob());
86			d->setChi2(fit.chisq());
87			d->setNdof(fit.ndof());
88
89			FourVector p4Fitted(0.,0.,0.,0.);
90			p4Fitted += fit.getTrackP4(1);
91			p4Fitted += fit.getTrackP4(2);
92			d->setFourMomentum(p4Fitted);
93			d->setPosition (fit.getVertex (MultiVertexFitterD::VERTEX_1));
94			d->setError (fit.getErrorMatrix(MultiVertexFitterD::VERTEX_1));
95			float mass, massErr;
96			const int trksIds[2] = { 1, 2 };
97			mass = fit.getMass(2,trksIds,massErr);
98
99			ThreeVector p3Fitted(p4Fitted.px(),p4Fitted.py(),p4Fitted.pz());
100	bendavid	1.8
101	paus	1.12	// Get decay length in xy plane
102	bendavid	1.8	float dl, dlErr;
103	paus	1.13	dl = fit.getDecayLength (MultiVertexFitterD::PRIMARY_VERTEX,MultiVertexFitterD::VERTEX_1,
104			p3Fitted,dlErr);
105
106			// Get Z decay length
107	bendavid	1.8	float dlz, dlzErr;
108	paus	1.13	dlz = fit.getZDecayLength(MultiVertexFitterD::PRIMARY_VERTEX,MultiVertexFitterD::VERTEX_1,
109			p3Fitted, dlzErr);
110
111			// Get impact parameter
112	bendavid	1.8	float dxy, dxyErr;
113	paus	1.13	dxy = fit.getImpactPar (MultiVertexFitterD::PRIMARY_VERTEX,MultiVertexFitterD::VERTEX_1,
114			p3Fitted, dxyErr);
115	loizides	1.4
116	bendavid	1.11	BasePartPtr ptr1(hStables1,i);
117	paus	1.13	BasePartPtr ptr2(hStables2,j);
118
119			StableData c1(fit.getTrackP4(1).px(),fit.getTrackP4(1).py(), fit.getTrackP4(1).pz(), ptr1);
120			StableData c2(fit.getTrackP4(2).px(),fit.getTrackP4(2).py(), fit.getTrackP4(2).pz(), ptr2);
121
122			const ThreeVector vtxPos = fit.getVertex(MultiVertexFitterD::VERTEX_1);
123	paus	1.12	const ThreeVector trkMom1(fit.getTrackP4(1).px(),fit.getTrackP4(1).py(),
124	paus	1.13	fit.getTrackP4(1).pz());
125	paus	1.12	const ThreeVector trkMom2(fit.getTrackP4(2).px(),fit.getTrackP4(2).py(),
126	paus	1.13	fit.getTrackP4(2).pz());
127
128	paus	1.12	// Build corrected HitPattern for StableData, removing hits before the fit vertex
129	bendavid	1.11	reco::HitPattern hits1 = dropper->CorrectedHits(s1.track(), vtxPos, trkMom1, dlErr, dlzErr);
130			reco::HitPattern hits2 = dropper->CorrectedHits(s2.track(), vtxPos, trkMom2, dlErr, dlzErr);
131	paus	1.13
132	bendavid	1.11	c1.SetHits(hits1);
133			c2.SetHits(hits2);
134	bendavid	1.14	c1.SetHitsFilled();
135			c2.SetHitsFilled();
136	paus	1.13
137			d->addStableChild (c1);
138			d->addStableChild (c2);
139	bendavid	1.8	d->setFittedMass (mass);
140			d->setFittedMassError(massErr);
141	paus	1.13	d->setLxy (dl);
142			d->setLxyError (dlErr);
143			d->setLxyToPv (dl);
144			d->setLxyToPvError (dlErr);
145			d->setLz (dlz);
146			d->setLzError (dlzErr);
147			d->setLzToPv (dlz);
148			d->setLzToPvError (dlzErr);
149			d->setDxy (dxy);
150			d->setDxyError (dxyErr);
151			d->setDxyToPv (dxy);
152			d->setDxyToPvError (dxyErr);
153
154			// Put the result into our collection
155			pD->push_back(*d);
156	bendavid	1.1	delete d;
157			}
158			}
159			}
160
161			// Write the collection even if it is empty
162	paus	1.13	if (0)
163	loizides	1.10	cout << " ProducerD2SS::produce - " << pD->size() << " entries collection created -"
164			<< " (Pid: " << oPid_ << ")\n";
165	bendavid	1.1	evt.put(pD);
166			}
167
168	paus	1.12	// Define this as a plug-in
169	bendavid	1.1	DEFINE_FWK_MODULE(ProducerD2SS);