Producers/src/ProducerConversions.cc

// $Id: ProducerConversions.cc,v 1.4 2008/09/24 09:00:54 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 "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/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","")),
  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();

  // 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");
      }
      if (fit.fit()) {
        DecayPart *d = new DecayPart(oPid_,DecayPart::Fast);
        
        BasePartPtr ptr1(hStables1,i);
        BasePartPtr ptr2(hStables2,j);
        d->addChild(ptr1);
        d->addChild(ptr2);
        d->addChildMom(fit.getTrackP4(1));
        d->addChildMom(fit.getTrackP4(2));

        // 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);
        
        if(0) {
          const reco::Track *p1 = s1.track();
          const reco::Track *p2 = s2.track();
          // create the dimuon system
          FourVector mu1(p1->px(),p1->py(),p1->pz(),sqrt(p1->p()*p1->p()+0.105658357*0.105658357));
          FourVector mu2(p2->px(),p2->py(),p2->pz(),sqrt(p2->p()*p2->p()+0.105658357*0.105658357));
          FourVector diMu = mu1+mu2;

          // for convenience and economy
          double mass4Vec = sqrt(diMu.M2());
        
          printf(" Generated mass:   ....\n");
          printf(" Four vector mass: %14.6f\n",mass4Vec);
          printf(" Fitted mass:      %14.6f +- %14.6f\n",mass,massErr);
        }

        d->setFittedMass(mass);
        d->setFittedMassError(massErr);

        // 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.5
Committed:	Sat Sep 27 05:48:25 2008 UTC (16 years, 7 months ago) by loizides
Content type:	text/plain
Branch:	MAIN
Changes since 1.4:	+27 -38 lines
Log Message:	Cleanup
#	User	Rev	Content
1	loizides	1.5	// $Id: ProducerConversions.cc,v 1.4 2008/09/24 09:00:54 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			#include "MitEdm/DataFormats/interface/Types.h"
8	loizides	1.5	#include "MitEdm/DataFormats/interface/Collections.h"
9	bendavid	1.1	#include "MitEdm/DataFormats/interface/DecayPart.h"
10			#include "MitEdm/DataFormats/interface/StablePart.h"
11			#include "MitEdm/VertexFitInterface/interface/MvfInterface.h"
12
13			using namespace std;
14			using namespace edm;
15			using namespace reco;
16			using namespace mitedm;
17			using namespace mithep;
18
19			//--------------------------------------------------------------------------------------------------
20			ProducerConversions::ProducerConversions(const ParameterSet& cfg) :
21	loizides	1.5	BaseCandProducer(cfg),
22			iStables1_(cfg.getUntrackedParameter<string>("iStables1","")),
23			iStables2_(cfg.getUntrackedParameter<string>("iStables2","")),
24			convConstraint_(cfg.getUntrackedParameter<bool>("convConstraint",false)),
25			convConstraint3D_(cfg.getUntrackedParameter<bool>("convConstraint3D",true)),
26			rhoMin_(cfg.getUntrackedParameter<double>("rhoMin",0.0))
27	bendavid	1.1	{
28	loizides	1.5	// Constructor.
29
30	bendavid	1.1	produces<DecayPartCol>();
31			}
32
33			//--------------------------------------------------------------------------------------------------
34			ProducerConversions::~ProducerConversions()
35			{
36	loizides	1.5	// Destructor.
37	bendavid	1.1	}
38
39			//--------------------------------------------------------------------------------------------------
40			void ProducerConversions::produce(Event &evt, const EventSetup &setup)
41			{
42	loizides	1.5	// Produce our DecayPartCol.
43
44	bendavid	1.1	// First input collection
45			Handle<StablePartCol> hStables1;
46			if (!GetProduct(iStables1_, hStables1, evt))
47			return;
48			const StablePartCol *pS1 = hStables1.product();
49			// Second input collection
50			Handle<StablePartCol> hStables2;
51			if (!GetProduct(iStables2_, hStables2, evt))
52			return;
53			const StablePartCol *pS2 = hStables2.product();
54
55			// Create the output collection
56			auto_ptr<DecayPartCol> pD(new DecayPartCol());
57
58			// Simple double loop
59			for (UInt_t i = 0; i<pS1->size(); ++i) {
60			const StablePart &s1 = pS1->at(i);
61
62			UInt_t j;
63			if (iStables1_ == iStables2_)
64			j = i+1;
65			else
66			j = 0;
67
68			for (; j<pS2->size(); ++j) {
69			const StablePart &s2 = pS2->at(j);
70
71			// Vertex fit now, possibly with conversion constraint
72			MultiVertexFitter fit;
73	loizides	1.5	fit.init(3.8); // Reset to the MC magnetic field of 3.8 Tesla
74	bendavid	1.1	MvfInterface fitInt(&fit);
75			fitInt.addTrack(s1.track(),1,s1.mass(),MultiVertexFitter::VERTEX_1);
76			fitInt.addTrack(s2.track(),2,s2.mass(),MultiVertexFitter::VERTEX_1);
77	bendavid	1.3	if (convConstraint3D_) {
78			fit.conversion_3d(MultiVertexFitter::VERTEX_1);
79			//printf("applying 3d conversion constraint\n");
80			}
81			else if (convConstraint_) {
82			fit.conversion_2d(MultiVertexFitter::VERTEX_1);
83			//printf("applying 2d conversion constraint\n");
84			}
85	bendavid	1.1	if (fit.fit()) {
86	bendavid	1.3	DecayPart *d = new DecayPart(oPid_,DecayPart::Fast);
87	bendavid	1.1
88	bendavid	1.4	BasePartPtr ptr1(hStables1,i);
89			BasePartPtr ptr2(hStables2,j);
90			d->addChild(ptr1);
91			d->addChild(ptr2);
92	bendavid	1.2	d->addChildMom(fit.getTrackP4(1));
93			d->addChildMom(fit.getTrackP4(2));
94	loizides	1.5
95	bendavid	1.3	// Update temporarily some of the quantities (prob, chi2, nDoF, mass, lxy, pt, fourMomentum)
96			d->setProb(fit.prob());
97			d->setChi2(fit.chisq());
98			d->setNdof(fit.ndof());
99
100			FourVector p4Fitted(0.,0.,0.,0.);
101			p4Fitted += fit.getTrackP4(1);
102			p4Fitted += fit.getTrackP4(2);
103			d->setFourMomentum(p4Fitted);
104			d->setPosition(fit.getVertex (MultiVertexFitter::VERTEX_1));
105			d->setError (fit.getErrorMatrix(MultiVertexFitter::VERTEX_1));
106			float mass, massErr;
107			const int trksIds[2] = { 1, 2 };
108			mass = fit.getMass(2,trksIds,massErr);
109
110	bendavid	1.1	if(0) {
111			const reco::Track *p1 = s1.track();
112			const reco::Track *p2 = s2.track();
113	loizides	1.5	// create the dimuon system
114	bendavid	1.1	FourVector mu1(p1->px(),p1->py(),p1->pz(),sqrt(p1->p()p1->p()+0.1056583570.105658357));
115			FourVector mu2(p2->px(),p2->py(),p2->pz(),sqrt(p2->p()p2->p()+0.1056583570.105658357));
116			FourVector diMu = mu1+mu2;
117
118	loizides	1.5	// for convenience and economy
119	bendavid	1.1	double mass4Vec = sqrt(diMu.M2());
120	bendavid	1.3
121	bendavid	1.1	printf(" Generated mass: ....\n");
122			printf(" Four vector mass: %14.6f\n",mass4Vec);
123			printf(" Fitted mass: %14.6f +- %14.6f\n",mass,massErr);
124			}
125
126	loizides	1.5	d->setFittedMass(mass);
127			d->setFittedMassError(massErr);
128
129			// Put the result into our collection
130	bendavid	1.1	if (d->position().rho() > rhoMin_)
131	bendavid	1.3	pD->push_back(*d);
132	loizides	1.5
133	bendavid	1.1	delete d;
134			}
135			}
136			}
137
138			// Write the collection even if it is empty
139	loizides	1.5	if (0) {
140			cout << " ProducerConversions::produce - " << pD->size() << " entries collection created -"
141			<< " (Pid: " << oPid_ << ")\n";
142			}
143	bendavid	1.1	evt.put(pD);
144			}
145
146			//define this as a plug-in
147			DEFINE_FWK_MODULE(ProducerConversions);