Producers/src/ProducerConversions.cc

// $Id: ProducerConversions.cc,v 1.23 2010/01/18 14:41:33 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 "MagneticField/Engine/interface/MagneticField.h"
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
#include "TrackingTools/PatternTools/interface/ClosestApproachInRPhi.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/TrackParameters.h"
#include <TMath.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)),
  useRhoMin_       (cfg.getUntrackedParameter<bool>  ("useRhoMin",true)),
  useHitDropper_   (cfg.getUntrackedParameter<bool>  ("useHitDropper",true)),
  applyChargeConstraint_(cfg.getUntrackedParameter<bool>  ("applyChargeConstraint",false))
{
  // 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();
  
  //Get Magnetic Field from event setup, taking value at (0,0,0)
  edm::ESHandle<MagneticField> magneticField;
  setup.get<IdealMagneticFieldRecord>().get(magneticField);
  const double bfield = magneticField->inTesla(GlobalPoint(0.,0.,0.)).z();
  
  edm::ESHandle<TransientTrackBuilder> hTransientTrackBuilder;
  setup.get<TransientTrackRecord>().get("TransientTrackBuilder",hTransientTrackBuilder);
  const TransientTrackBuilder *transientTrackBuilder = hTransientTrackBuilder.product();
  
  //construct intermediate collection of TrackParameters in mvf format for vertex fit
  std::vector<TrackParameters> trkPars1;
  for (UInt_t i = 0; i<pS1->size(); ++i) {
    const reco::Track *t = pS1->at(i).track();
    const reco::TransientTrack ttrk = transientTrackBuilder->build(t);
    TrackParameters cmsTrk(ttrk);
    TrackParameters mvfTrk = cmsTrk.mvfTrack(); 
    trkPars1.push_back(mvfTrk);
  }
  
  bool sameCollection = iStables1_ == iStables2_;
  
  std::vector<TrackParameters> trkPars2;
  if (sameCollection)
    trkPars2 = trkPars1;
  else for (UInt_t i = 0; i<pS2->size(); ++i) {
    const reco::Track *t = pS2->at(i).track();
    const reco::TransientTrack ttrk = transientTrackBuilder->build(t);
    TrackParameters cmsTrk(ttrk);
    TrackParameters mvfTrk = cmsTrk.mvfTrack(); 
    trkPars2.push_back(mvfTrk);
  }

  
  // Create the output collection
  auto_ptr<DecayPartCol> pD(new DecayPartCol());
  
  ClosestApproachInRPhi helixIntersector;
  
  int nFits = 0;
  
  //printf("S1 size = %i\n", pS1->size());
  
  std::set<std::pair<edm::Ptr<reco::Track>,edm::Ptr<reco::Track> > > pairSet;
  // Simple double loop
  for (UInt_t i = 0; i<pS1->size(); ++i) {
    const StablePart &s1 =  pS1->at(i);
    
    const reco::Track * t1 = s1.track();
    const TrackParameters &trkPar1 = trkPars1.at(i);
    
    UInt_t j;
    if (sameCollection)
      j = i+1; 
    else
      j = 0;
    
    TrajectoryStateTransform tsTransform;
    
    FreeTrajectoryState initialState1 = tsTransform.initialFreeState(*s1.track(),&*magneticField);
    
    for (; j<pS2->size(); ++j) {
      const StablePart &s2 = pS2->at(j);
      
      if (!sameCollection) {
        std::pair<edm::Ptr<reco::Track>,edm::Ptr<reco::Track> > pair1(s1.trackPtr(),s2.trackPtr());
        std::pair<edm::Ptr<reco::Track>,edm::Ptr<reco::Track> > pair2(s2.trackPtr(),s1.trackPtr());
        if (pairSet.find(pair1)!=pairSet.end() || pairSet.find(pair2)!=pairSet.end()) {
          continue;
        }
        
        pairSet.insert(pair1);
        pairSet.insert(pair2);
      }
      
      //Do fast helix fit to check if there's any hope
      const reco::Track * t2 = s2.track();
      const TrackParameters &trkPar2 = trkPars2.at(j);
      
      int trackCharge = t1->charge() + t2->charge();
      
      double dR0 = 0.0;
      
      if (!applyChargeConstraint_ || trackCharge==0) {
        FreeTrajectoryState initialState2 = tsTransform.initialFreeState(*s2.track(),&*magneticField);
        helixIntersector.calculate(initialState1, initialState2);
        if (helixIntersector.status())
          dR0 = helixIntersector.crossingPoint().perp();
      }

      int fitStatus = 0;
      MultiVertexFitterD fit;      
      if ( (!applyChargeConstraint_ || trackCharge==0) && (!useRhoMin_ || dR0 > rhoMin_) ) {
      
        // Vertex fit now, possibly with conversion constraint
        nFits++;
  
        fit.init(bfield); // Reset to the magnetic field from the event setup
        
        fit.addTrack(*trkPar1.pars(),*trkPar1.cMat(),1,s1.mass(),MultiVertexFitterD::VERTEX_1);
        fit.addTrack(*trkPar2.pars(),*trkPar2.cMat(),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
        if (useHitDropper_) {
          std::pair<reco::HitPattern,uint> hits1 = dropper->CorrectedHitsAOD(s1.track(), vtxPos, trkMom1, dlErr, dlzErr);
          std::pair<reco::HitPattern,uint> hits2 = dropper->CorrectedHitsAOD(s2.track(), vtxPos, trkMom2, dlErr, dlzErr);                 
   
          c1.SetHits(hits1.first);
          c2.SetHits(hits2.first);
          c1.SetHitsFilled();
          c2.SetHitsFilled();
          c1.SetNWrongHits(hits1.second);
          c2.SetNWrongHits(hits2.second);
          
          reco::HitPattern sharedHits = dropper->SharedHits(s1.track(),s2.track());
          d->setSharedHits(sharedHits);
        }
        
        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;
      }
    }
  }
  
  //printf("nConversionFits = %i\n",nFits);

  // 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.24
Committed:	Tue Jun 8 20:17:30 2010 UTC (14 years, 11 months ago) by bendavid
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_025c_branch2, Mit_025c_branch1, Mit_025c_branch0, 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
Branch point for:	Mit_025c_branch
Changes since 1.23:	+24 -8 lines
Log Message:	Add simple tool to do conversion selection and matching
#	User	Rev	Content
1	bendavid	1.24	// $Id: ProducerConversions.cc,v 1.23 2010/01/18 14:41:33 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.18	#include "MagneticField/Engine/interface/MagneticField.h"
10			#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
11	bendavid	1.19	#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
12			#include "TrackingTools/PatternTools/interface/ClosestApproachInRPhi.h"
13	bendavid	1.7	#include "MitEdm/Producers/interface/HitDropperRecord.h"
14			#include "MitEdm/Producers/interface/HitDropper.h"
15	bendavid	1.1	#include "MitEdm/DataFormats/interface/Types.h"
16	loizides	1.5	#include "MitEdm/DataFormats/interface/Collections.h"
17	bendavid	1.1	#include "MitEdm/DataFormats/interface/DecayPart.h"
18			#include "MitEdm/DataFormats/interface/StablePart.h"
19	bendavid	1.6	#include "MitEdm/DataFormats/interface/StableData.h"
20	bendavid	1.1	#include "MitEdm/VertexFitInterface/interface/MvfInterface.h"
21	bendavid	1.23	#include "MitEdm/VertexFitInterface/interface/TrackParameters.h"
22	bendavid	1.20	#include <TMath.h>
23	bendavid	1.1
24			using namespace std;
25			using namespace edm;
26			using namespace reco;
27			using namespace mitedm;
28			using namespace mithep;
29
30			//--------------------------------------------------------------------------------------------------
31			ProducerConversions::ProducerConversions(const ParameterSet& cfg) :
32	paus	1.11	BaseCandProducer (cfg),
33			iStables1_ (cfg.getUntrackedParameter<string>("iStables1","")),
34			iStables2_ (cfg.getUntrackedParameter<string>("iStables2","")),
35			iPVertexes_ (cfg.getUntrackedParameter<string>("iPVertexes","offlinePrimaryVerticesWithBS")),
36			usePVertex_ (cfg.getUntrackedParameter<bool> ("usePVertex",true)),
37			convConstraint_ (cfg.getUntrackedParameter<bool> ("convConstraint",false)),
38			convConstraint3D_(cfg.getUntrackedParameter<bool> ("convConstraint3D",true)),
39	bendavid	1.15	rhoMin_ (cfg.getUntrackedParameter<double>("rhoMin",0.0)),
40	bendavid	1.23	useRhoMin_ (cfg.getUntrackedParameter<bool> ("useRhoMin",true)),
41	bendavid	1.20	useHitDropper_ (cfg.getUntrackedParameter<bool> ("useHitDropper",true)),
42	bendavid	1.21	applyChargeConstraint_(cfg.getUntrackedParameter<bool> ("applyChargeConstraint",false))
43	bendavid	1.1	{
44	loizides	1.5	// Constructor.
45
46	bendavid	1.1	produces<DecayPartCol>();
47			}
48
49			//--------------------------------------------------------------------------------------------------
50			ProducerConversions::~ProducerConversions()
51			{
52	loizides	1.5	// Destructor.
53	bendavid	1.1	}
54
55			//--------------------------------------------------------------------------------------------------
56			void ProducerConversions::produce(Event &evt, const EventSetup &setup)
57			{
58	loizides	1.5	// Produce our DecayPartCol.
59
60	bendavid	1.1	// First input collection
61			Handle<StablePartCol> hStables1;
62	bendavid	1.9	if (!GetProduct(iStables1_, hStables1, evt)) {
63			printf("Stable collection 1 not found in ProducerConversions\n");
64	bendavid	1.1	return;
65	bendavid	1.9	}
66	bendavid	1.1	const StablePartCol *pS1 = hStables1.product();
67			// Second input collection
68			Handle<StablePartCol> hStables2;
69	bendavid	1.9	if (!GetProduct(iStables2_, hStables2, evt)) {
70			printf("Stable collection 2 not found in ProducerConversions\n");
71	bendavid	1.1	return;
72	bendavid	1.9	}
73	bendavid	1.1	const StablePartCol *pS2 = hStables2.product();
74
75	bendavid	1.6	const reco::Vertex *vertex = 0;
76			mitedm::VertexPtr vPtr;
77			if (usePVertex_) {
78			// Primary vertex collection
79			Handle<reco::VertexCollection> hVertexes;
80			if (!GetProduct(iPVertexes_, hVertexes, evt))
81			return;
82			const reco::VertexCollection *pV = hVertexes.product();
83
84			//Choose the primary vertex with the largest number of tracks
85			UInt_t maxTracks=0;
86			for (UInt_t i=0; i<pV->size(); ++i) {
87			const reco::Vertex &v = pV->at(i);
88			UInt_t nTracks = v.tracksSize();
89			if (nTracks >= maxTracks) {
90			maxTracks = nTracks;
91			vertex = &v;
92			vPtr = mitedm::VertexPtr(hVertexes,i);
93			}
94			}
95			}
96	bendavid	1.7
97	paus	1.11	// Get hit dropper
98	bendavid	1.7	ESHandle<HitDropper> hDropper;
99			setup.get<HitDropperRecord>().get("HitDropper",hDropper);
100			const HitDropper *dropper = hDropper.product();
101	bendavid	1.6
102	bendavid	1.18	//Get Magnetic Field from event setup, taking value at (0,0,0)
103			edm::ESHandle<MagneticField> magneticField;
104			setup.get<IdealMagneticFieldRecord>().get(magneticField);
105			const double bfield = magneticField->inTesla(GlobalPoint(0.,0.,0.)).z();
106
107	bendavid	1.23	edm::ESHandle<TransientTrackBuilder> hTransientTrackBuilder;
108			setup.get<TransientTrackRecord>().get("TransientTrackBuilder",hTransientTrackBuilder);
109			const TransientTrackBuilder *transientTrackBuilder = hTransientTrackBuilder.product();
110
111			//construct intermediate collection of TrackParameters in mvf format for vertex fit
112			std::vector<TrackParameters> trkPars1;
113			for (UInt_t i = 0; i<pS1->size(); ++i) {
114			const reco::Track *t = pS1->at(i).track();
115			const reco::TransientTrack ttrk = transientTrackBuilder->build(t);
116			TrackParameters cmsTrk(ttrk);
117			TrackParameters mvfTrk = cmsTrk.mvfTrack();
118			trkPars1.push_back(mvfTrk);
119			}
120
121	bendavid	1.24	bool sameCollection = iStables1_ == iStables2_;
122
123	bendavid	1.23	std::vector<TrackParameters> trkPars2;
124	bendavid	1.24	if (sameCollection)
125	bendavid	1.23	trkPars2 = trkPars1;
126			else for (UInt_t i = 0; i<pS2->size(); ++i) {
127			const reco::Track *t = pS2->at(i).track();
128			const reco::TransientTrack ttrk = transientTrackBuilder->build(t);
129			TrackParameters cmsTrk(ttrk);
130			TrackParameters mvfTrk = cmsTrk.mvfTrack();
131			trkPars2.push_back(mvfTrk);
132			}
133
134
135	bendavid	1.1	// Create the output collection
136			auto_ptr<DecayPartCol> pD(new DecayPartCol());
137	bendavid	1.6
138	bendavid	1.19	ClosestApproachInRPhi helixIntersector;
139
140	bendavid	1.20	int nFits = 0;
141
142			//printf("S1 size = %i\n", pS1->size());
143
144	bendavid	1.24	std::set<std::pair<edm::Ptr<reco::Track>,edm::Ptr<reco::Track> > > pairSet;
145	bendavid	1.1	// Simple double loop
146			for (UInt_t i = 0; i<pS1->size(); ++i) {
147			const StablePart &s1 = pS1->at(i);
148
149	bendavid	1.20	const reco::Track * t1 = s1.track();
150	bendavid	1.23	const TrackParameters &trkPar1 = trkPars1.at(i);
151
152	bendavid	1.1	UInt_t j;
153	bendavid	1.24	if (sameCollection)
154	bendavid	1.1	j = i+1;
155			else
156			j = 0;
157
158	bendavid	1.19	TrajectoryStateTransform tsTransform;
159
160			FreeTrajectoryState initialState1 = tsTransform.initialFreeState(s1.track(),&magneticField);
161
162	bendavid	1.1	for (; j<pS2->size(); ++j) {
163			const StablePart &s2 = pS2->at(j);
164	bendavid	1.24
165			if (!sameCollection) {
166			std::pair<edm::Ptr<reco::Track>,edm::Ptr<reco::Track> > pair1(s1.trackPtr(),s2.trackPtr());
167			std::pair<edm::Ptr<reco::Track>,edm::Ptr<reco::Track> > pair2(s2.trackPtr(),s1.trackPtr());
168			if (pairSet.find(pair1)!=pairSet.end() \|\| pairSet.find(pair2)!=pairSet.end()) {
169			continue;
170			}
171
172			pairSet.insert(pair1);
173			pairSet.insert(pair2);
174			}
175
176	bendavid	1.21	//Do fast helix fit to check if there's any hope
177	bendavid	1.13	const reco::Track * t2 = s2.track();
178	bendavid	1.23	const TrackParameters &trkPar2 = trkPars2.at(j);
179	bendavid	1.13
180			int trackCharge = t1->charge() + t2->charge();
181	bendavid	1.19
182	bendavid	1.13	double dR0 = 0.0;
183
184	bendavid	1.20	if (!applyChargeConstraint_ \|\| trackCharge==0) {
185	bendavid	1.19	FreeTrajectoryState initialState2 = tsTransform.initialFreeState(s2.track(),&magneticField);
186			helixIntersector.calculate(initialState1, initialState2);
187			if (helixIntersector.status())
188			dR0 = helixIntersector.crossingPoint().perp();
189	bendavid	1.3	}
190	bendavid	1.13
191			int fitStatus = 0;
192			MultiVertexFitterD fit;
193	bendavid	1.23	if ( (!applyChargeConstraint_ \|\| trackCharge==0) && (!useRhoMin_ \|\| dR0 > rhoMin_) ) {
194	bendavid	1.13
195			// Vertex fit now, possibly with conversion constraint
196	bendavid	1.20	nFits++;
197	bendavid	1.13
198	bendavid	1.18	fit.init(bfield); // Reset to the magnetic field from the event setup
199	bendavid	1.23
200			fit.addTrack(trkPar1.pars(),trkPar1.cMat(),1,s1.mass(),MultiVertexFitterD::VERTEX_1);
201			fit.addTrack(trkPar2.pars(),trkPar2.cMat(),2,s2.mass(),MultiVertexFitterD::VERTEX_1);
202	bendavid	1.13	if (convConstraint3D_) {
203			fit.conversion_3d(MultiVertexFitterD::VERTEX_1);
204			//printf("applying 3d conversion constraint\n");
205			}
206			else if (convConstraint_) {
207			fit.conversion_2d(MultiVertexFitterD::VERTEX_1);
208			//printf("applying 2d conversion constraint\n");
209			}
210			//initialize primary vertex parameters in the fitter
211			if (usePVertex_) {
212			float vErr[3][3];
213			for (UInt_t vi=0; vi<3; ++vi)
214			for (UInt_t vj=0; vj<3; ++vj)
215			vErr[vi][vj] = vertex->covariance(vi,vj);
216
217			fit.setPrimaryVertex(vertex->x(),vertex->y(),vertex->z());
218			fit.setPrimaryVertexError(vErr);
219			}
220
221			fitStatus = fit.fit();
222	bendavid	1.23
223	bendavid	1.6	}
224
225			if (fitStatus) {
226	bendavid	1.3	DecayPart *d = new DecayPart(oPid_,DecayPart::Fast);
227	bendavid	1.1
228	bendavid	1.3	// Update temporarily some of the quantities (prob, chi2, nDoF, mass, lxy, pt, fourMomentum)
229			d->setProb(fit.prob());
230			d->setChi2(fit.chisq());
231			d->setNdof(fit.ndof());
232
233			FourVector p4Fitted(0.,0.,0.,0.);
234			p4Fitted += fit.getTrackP4(1);
235			p4Fitted += fit.getTrackP4(2);
236			d->setFourMomentum(p4Fitted);
237	paus	1.12	d->setPosition(fit.getVertex (MultiVertexFitterD::VERTEX_1));
238			d->setError (fit.getErrorMatrix(MultiVertexFitterD::VERTEX_1));
239	bendavid	1.3	float mass, massErr;
240			const int trksIds[2] = { 1, 2 };
241			mass = fit.getMass(2,trksIds,massErr);
242
243	bendavid	1.6	ThreeVector p3Fitted(p4Fitted.px(), p4Fitted.py(), p4Fitted.pz());
244
245	paus	1.11	// Get decay length in xy plane
246	bendavid	1.6	float dl, dlErr;
247	paus	1.12	dl = fit.getDecayLength(MultiVertexFitterD::PRIMARY_VERTEX, MultiVertexFitterD::VERTEX_1,
248	bendavid	1.6	p3Fitted, dlErr);
249
250	paus	1.11	// Get Z decay length
251	bendavid	1.6	float dlz, dlzErr;
252	paus	1.12	dlz = fit.getZDecayLength(MultiVertexFitterD::PRIMARY_VERTEX, MultiVertexFitterD::VERTEX_1,
253	bendavid	1.6	p3Fitted, dlzErr);
254
255	paus	1.11	// Get impact parameter
256	bendavid	1.6	float dxy, dxyErr;
257	paus	1.12	dxy = fit.getImpactPar(MultiVertexFitterD::PRIMARY_VERTEX, MultiVertexFitterD::VERTEX_1,
258	bendavid	1.6	p3Fitted, dxyErr);
259	bendavid	1.1
260	bendavid	1.8	BasePartPtr ptr1(hStables1,i);
261			BasePartPtr ptr2(hStables2,j);
262
263			StableData c1(fit.getTrackP4(1).px(),fit.getTrackP4(1).py(), fit.getTrackP4(1).pz(), ptr1);
264			StableData c2(fit.getTrackP4(2).px(),fit.getTrackP4(2).py(), fit.getTrackP4(2).pz(), ptr2);
265
266	paus	1.12	const ThreeVector vtxPos = fit.getVertex(MultiVertexFitterD::VERTEX_1);
267	loizides	1.14	const ThreeVector trkMom1(fit.getTrackP4(1).px(),
268			fit.getTrackP4(1).py(),
269			fit.getTrackP4(1).pz());
270			const ThreeVector trkMom2(fit.getTrackP4(2).px(),
271			fit.getTrackP4(2).py(),
272			fit.getTrackP4(2).pz());
273	bendavid	1.8
274	paus	1.11	// Build corrected HitPattern for StableData, removing hits before the fit vertex
275	bendavid	1.15	if (useHitDropper_) {
276	bendavid	1.24	std::pair<reco::HitPattern,uint> hits1 = dropper->CorrectedHitsAOD(s1.track(), vtxPos, trkMom1, dlErr, dlzErr);
277			std::pair<reco::HitPattern,uint> hits2 = dropper->CorrectedHitsAOD(s2.track(), vtxPos, trkMom2, dlErr, dlzErr);
278	bendavid	1.17
279	bendavid	1.24	c1.SetHits(hits1.first);
280			c2.SetHits(hits2.first);
281	bendavid	1.15	c1.SetHitsFilled();
282			c2.SetHitsFilled();
283	bendavid	1.24	c1.SetNWrongHits(hits1.second);
284			c2.SetNWrongHits(hits2.second);
285	bendavid	1.22
286			reco::HitPattern sharedHits = dropper->SharedHits(s1.track(),s2.track());
287			d->setSharedHits(sharedHits);
288	bendavid	1.15	}
289	bendavid	1.8
290			d->addStableChild(c1);
291			d->addStableChild(c2);
292
293
294	bendavid	1.6	d->setFittedMass (mass);
295			d->setFittedMassError(massErr);
296
297			d->setLxy(dl);
298			d->setLxyError(dlErr);
299			d->setLxyToPv(dl);
300			d->setLxyToPvError(dlErr);
301
302			d->setLz(dlz);
303			d->setLzError(dlzErr);
304			d->setLzToPv(dlz);
305			d->setLzToPvError(dlzErr);
306
307			d->setDxy(dxy);
308			d->setDxyError(dxyErr);
309			d->setDxyToPv(dxy);
310			d->setDxyToPvError(dxyErr);
311
312			if (usePVertex_)
313			d->setPrimaryVertex(vPtr);
314
315			// Put the result into our collection
316	bendavid	1.13	pD->push_back(*d);
317	loizides	1.5
318	bendavid	1.1	delete d;
319			}
320			}
321			}
322	bendavid	1.20
323			//printf("nConversionFits = %i\n",nFits);
324	bendavid	1.1
325			// Write the collection even if it is empty
326	loizides	1.5	if (0) {
327			cout << " ProducerConversions::produce - " << pD->size() << " entries collection created -"
328			<< " (Pid: " << oPid_ << ")\n";
329			}
330	bendavid	1.1	evt.put(pD);
331			}
332
333			//define this as a plug-in
334			DEFINE_FWK_MODULE(ProducerConversions);