Producers/src/ProducerConversions.cc

// $Id: ProducerConversions.cc,v 1.2 2008/09/19 12:00:05 bendavid Exp $

#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/CollectionsEdm.h"
#include "MitEdm/DataFormats/interface/DecayPart.h"
#include "MitEdm/DataFormats/interface/StablePart.h"
#include "MitEdm/VertexFitInterface/interface/MvfInterface.h"
#include "MitEdm/Producers/interface/ProducerConversions.h"

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

//--------------------------------------------------------------------------------------------------
ProducerConversions::ProducerConversions(const ParameterSet& cfg) :
  BaseCandidate(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 ))
{
  produces<DecayPartCol>();
}

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

//--------------------------------------------------------------------------------------------------
void ProducerConversions::produce(Event &evt, const EventSetup &setup)
{
  // -----------------------------------------------------------------------------------------------
  // Get the input
  // -----------------------------------------------------------------------------------------------
  // 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 4 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);
        
        RefToBaseProd<BasePart> baseRef1(hStables1);
        RefToBaseProd<BasePart> baseRef2(hStables2);
        BasePartBaseRef ref1(baseRef1,i);
        BasePartBaseRef ref2(baseRef2,j);
        d->addChild(ref1);
        d->addChild(ref2);
        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
  // -----------------------------------------------------------------------------------------------
 // cout << " ProducerConversions::produce - " << pD->size() << " entries collection created -"
  //     << " (Pid: " << oPid_ << ")\n";
  evt.put(pD);
}

//--------------------------------------------------------------------------------------------------
void ProducerConversions::beginJob(const EventSetup &setup)
{
}

//--------------------------------------------------------------------------------------------------
void ProducerConversions::endJob()
{
}

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