Producers/src/ProducerV2SS.cc

// $Id: ProducerV2SS.cc,v 1.21 2010/01/18 14:41:33 bendavid Exp $

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

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

//--------------------------------------------------------------------------------------------------
ProducerV2SS::ProducerV2SS(const ParameterSet& cfg) :
  ProducerD2SS(cfg),
  rhoMin_     (cfg.getUntrackedParameter<double>("minRadius",   0.0)),
  massMin_    (cfg.getUntrackedParameter<double>("minMass",     0.0)),
  massMax_    (cfg.getUntrackedParameter<double>("maxMass",     3.0)),
  dZMax_      (cfg.getUntrackedParameter<double>("maxZDistance",5.0)),
  useHitDropper_(cfg.getUntrackedParameter<bool>("useHitDropper",true)),
  applyChargeConstraint_(cfg.getUntrackedParameter<bool>  ("applyChargeConstraint",false))
{
  // Constructor.
}

//--------------------------------------------------------------------------------------------------
ProducerV2SS::~ProducerV2SS()
{
  // Destructor.
}

//--------------------------------------------------------------------------------------------------
void ProducerV2SS::produce(Event &evt, const EventSetup &setup)
{
  // Produce the output collection.

  auto_ptr<DecayPartCol> pD(new DecayPartCol());

  // First input collection
  Handle<StablePartCol> hStables1;
  Handle<StablePartCol> hStables2;
  if (!GetProduct(iStables1_, hStables1, evt) ) {
    cout << "Couldn't get in collection in Producer V2SS" << endl;
    evt.put(pD);
    return;  
  }
  const StablePartCol *pS1 = hStables1.product();

  // Second input collection
  if(!GetProduct(iStables2_, hStables2, evt)) {
    cout << "Couldn't get in collection in Producer V2SS" << endl;
    evt.put(pD);
    return;  
  }
  const StablePartCol *pS2 = hStables2.product();
  
  //get hit dropper
  ESHandle<HitDropper> hDropper;
  const HitDropper *dropper = 0;
  if (useHitDropper_) {
    setup.get<HitDropperRecord>().get("HitDropper",hDropper);
    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);
  }
  
  std::vector<TrackParameters> trkPars2;
  if (iStables1_ == iStables2_)
    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);
  }

  // -----------------------------------------------------------------------------------------------
  // Simple double loop
  // -----------------------------------------------------------------------------------------------
  if (0)
    cout << "Starting V finder loop" << endl;

  ClosestApproachInRPhi helixIntersector;
  
  
  //sX_y: X= pion or proton collection.  
  //i, j = 2 loop particles. 
  //ex.: s1_i and s2_i are same particle as pion and proton
  for (UInt_t i=0; i<pS1->size(); ++i) {
    const StablePart &s1 =  pS1->at(i);
    const TrackParameters &trkPar1 = trkPars1.at(i);
   
    //const reco::Track * t1 = s1.track();
    
    UInt_t j;
    if (iStables1_ == iStables2_)
      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);
      const TrackParameters &trkPar2 = trkPars2.at(j);

      if( applyChargeConstraint_ && (s1.charge() + s2.charge() != 0) ) continue;

      // do fast helix fit to check if there's any hope
      //
      //const reco::Track * t2 = s2.track();
         
      double dZ0 = -999;
      double dR0 = -999;
      double mass0 = 0.0;
      
      FreeTrajectoryState initialState2 = tsTransform.initialFreeState(*s2.track(),&*magneticField);
      helixIntersector.calculate(initialState1, initialState2);
      if (helixIntersector.status()) {
        dZ0 = fabs(helixIntersector.points().first.z() - helixIntersector.points().second.z());
        dR0 = helixIntersector.crossingPoint().perp();
        
        GlobalVector     v1, v2;
        v1 = helixIntersector.trajectoryParameters().first.momentum();
        v2 = helixIntersector.trajectoryParameters().second.momentum();

        double e1 = sqrt(v1.mag2()+s1.mass()*s1.mass());
        double x1 = v1.x();
        double y1 = v1.y();
        double z1 = v1.z();

        double e2 = sqrt(v2.mag2()+s2.mass()*s2.mass());
        double x2 = v2.x();
        double y2 = v2.y();
        double z2 = v2.z();

        FourVector sum(x1+x2, y1+y2, z1+z2, e1+e2);

        mass0 = sqrt(sum.M2());
      }
      
      
      // Basic cuts on helix intersection
      if(mass0 > massMax_ || mass0<massMin_ || fabs(dZ0) > dZMax_ || dR0 < rhoMin_) continue;

      // -------------------------------------------------------------------------------------------
      // Do vertex fit for all pairs
      // -------------------------------------------------------------------------------------------
      mithep::MultiVertexFitterD fit;
      fit.init(bfield); // Reset to the magnetic field from the event setup
      fit.setChisqMax(100);
      
      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 (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
        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);

        // put the result into our collection
        pD->push_back(*d);
        delete d;
      }  //done processing fit
    } //end j loop
  } //end i loop

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

// define this as a plug-in
DEFINE_FWK_MODULE(ProducerV2SS);
Revision:	1.22
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.21:	+7 -5 lines
Log Message:	Add simple tool to do conversion selection and matching
#	User	Rev	Content
1	bendavid	1.22	// $Id: ProducerV2SS.cc,v 1.21 2010/01/18 14:41:33 bendavid Exp $
2	loizides	1.7
3	mrudolph	1.3	#include "DataFormats/Common/interface/Handle.h"
4			#include "DataFormats/TrackReco/interface/Track.h"
5			#include "DataFormats/TrackReco/interface/TrackFwd.h"
6	bendavid	1.15	#include "MagneticField/Engine/interface/MagneticField.h"
7			#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
8	bendavid	1.16	#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
9			#include "TrackingTools/PatternTools/interface/ClosestApproachInRPhi.h"
10	bendavid	1.8	#include "MitEdm/Producers/interface/HitDropperRecord.h"
11			#include "MitEdm/Producers/interface/HitDropper.h"
12	mrudolph	1.3	#include "MitEdm/DataFormats/interface/Types.h"
13	mrudolph	1.4	#include "MitEdm/DataFormats/interface/Collections.h"
14	mrudolph	1.3	#include "MitEdm/DataFormats/interface/DecayPart.h"
15			#include "MitEdm/DataFormats/interface/StablePart.h"
16	bendavid	1.5	#include "MitEdm/DataFormats/interface/StableData.h"
17	mrudolph	1.3	#include "MitEdm/VertexFitInterface/interface/MvfInterface.h"
18	bendavid	1.21	#include "MitEdm/VertexFitInterface/interface/TrackParameters.h"
19	mrudolph	1.3	#include "MitEdm/Producers/interface/ProducerV2SS.h"
20	bendavid	1.18	#include <TMath.h>
21	mrudolph	1.3
22			using namespace std;
23			using namespace edm;
24			using namespace mitedm;
25	bendavid	1.5	using namespace mithep;
26	mrudolph	1.3
27			//--------------------------------------------------------------------------------------------------
28			ProducerV2SS::ProducerV2SS(const ParameterSet& cfg) :
29			ProducerD2SS(cfg),
30	paus	1.10	rhoMin_ (cfg.getUntrackedParameter<double>("minRadius", 0.0)),
31			massMin_ (cfg.getUntrackedParameter<double>("minMass", 0.0)),
32			massMax_ (cfg.getUntrackedParameter<double>("maxMass", 3.0)),
33	bendavid	1.11	dZMax_ (cfg.getUntrackedParameter<double>("maxZDistance",5.0)),
34	bendavid	1.18	useHitDropper_(cfg.getUntrackedParameter<bool>("useHitDropper",true)),
35	bendavid	1.19	applyChargeConstraint_(cfg.getUntrackedParameter<bool> ("applyChargeConstraint",false))
36	mrudolph	1.3	{
37	loizides	1.12	// Constructor.
38	mrudolph	1.3	}
39
40			//--------------------------------------------------------------------------------------------------
41			ProducerV2SS::~ProducerV2SS()
42			{
43	loizides	1.12	// Destructor.
44	mrudolph	1.3	}
45
46			//--------------------------------------------------------------------------------------------------
47			void ProducerV2SS::produce(Event &evt, const EventSetup &setup)
48			{
49	loizides	1.12	// Produce the output collection.
50
51	mrudolph	1.3	auto_ptr<DecayPartCol> pD(new DecayPartCol());
52
53			// First input collection
54			Handle<StablePartCol> hStables1;
55			Handle<StablePartCol> hStables2;
56			if (!GetProduct(iStables1_, hStables1, evt) ) {
57			cout << "Couldn't get in collection in Producer V2SS" << endl;
58			evt.put(pD);
59			return;
60			}
61	loizides	1.7	const StablePartCol *pS1 = hStables1.product();
62	mrudolph	1.3
63	loizides	1.7	// Second input collection
64	mrudolph	1.3	if(!GetProduct(iStables2_, hStables2, evt)) {
65			cout << "Couldn't get in collection in Producer V2SS" << endl;
66			evt.put(pD);
67			return;
68			}
69			const StablePartCol *pS2 = hStables2.product();
70	bendavid	1.8
71			//get hit dropper
72			ESHandle<HitDropper> hDropper;
73	bendavid	1.11	const HitDropper *dropper = 0;
74			if (useHitDropper_) {
75			setup.get<HitDropperRecord>().get("HitDropper",hDropper);
76			dropper = hDropper.product();
77			}
78	bendavid	1.15
79			//Get Magnetic Field from event setup, taking value at (0,0,0)
80			edm::ESHandle<MagneticField> magneticField;
81			setup.get<IdealMagneticFieldRecord>().get(magneticField);
82			const double bfield = magneticField->inTesla(GlobalPoint(0.,0.,0.)).z();
83	bendavid	1.21
84			edm::ESHandle<TransientTrackBuilder> hTransientTrackBuilder;
85			setup.get<TransientTrackRecord>().get("TransientTrackBuilder",hTransientTrackBuilder);
86			const TransientTrackBuilder *transientTrackBuilder = hTransientTrackBuilder.product();
87
88			//construct intermediate collection of TrackParameters in mvf format for vertex fit
89			std::vector<TrackParameters> trkPars1;
90			for (UInt_t i = 0; i<pS1->size(); ++i) {
91			const reco::Track *t = pS1->at(i).track();
92			const reco::TransientTrack ttrk = transientTrackBuilder->build(t);
93			TrackParameters cmsTrk(ttrk);
94			TrackParameters mvfTrk = cmsTrk.mvfTrack();
95			trkPars1.push_back(mvfTrk);
96			}
97
98			std::vector<TrackParameters> trkPars2;
99			if (iStables1_ == iStables2_)
100			trkPars2 = trkPars1;
101			else for (UInt_t i = 0; i<pS2->size(); ++i) {
102			const reco::Track *t = pS2->at(i).track();
103			const reco::TransientTrack ttrk = transientTrackBuilder->build(t);
104			TrackParameters cmsTrk(ttrk);
105			TrackParameters mvfTrk = cmsTrk.mvfTrack();
106			trkPars2.push_back(mvfTrk);
107			}
108	mrudolph	1.3
109			// -----------------------------------------------------------------------------------------------
110			// Simple double loop
111			// -----------------------------------------------------------------------------------------------
112	loizides	1.7	if (0)
113			cout << "Starting V finder loop" << endl;
114	mrudolph	1.3
115	bendavid	1.16	ClosestApproachInRPhi helixIntersector;
116
117
118	loizides	1.7	//sX_y: X= pion or proton collection.
119			//i, j = 2 loop particles.
120			//ex.: s1_i and s2_i are same particle as pion and proton
121	mrudolph	1.3	for (UInt_t i=0; i<pS1->size(); ++i) {
122			const StablePart &s1 = pS1->at(i);
123	bendavid	1.21	const TrackParameters &trkPar1 = trkPars1.at(i);
124	mrudolph	1.3
125	bendavid	1.19	//const reco::Track * t1 = s1.track();
126	bendavid	1.18
127	mrudolph	1.3	UInt_t j;
128			if (iStables1_ == iStables2_)
129			j = i+1;
130			else
131			j = 0;
132
133	bendavid	1.16	TrajectoryStateTransform tsTransform;
134
135			FreeTrajectoryState initialState1 = tsTransform.initialFreeState(s1.track(),&magneticField);
136
137
138	mrudolph	1.3	for (; j<pS2->size(); ++j) {
139			const StablePart &s2 = pS2->at(j);
140	bendavid	1.21	const TrackParameters &trkPar2 = trkPars2.at(j);
141	mrudolph	1.3
142	bendavid	1.18	if( applyChargeConstraint_ && (s1.charge() + s2.charge() != 0) ) continue;
143	mrudolph	1.3
144	loizides	1.12	// do fast helix fit to check if there's any hope
145	bendavid	1.18	//
146	bendavid	1.19	//const reco::Track * t2 = s2.track();
147
148	mrudolph	1.3	double dZ0 = -999;
149			double dR0 = -999;
150	bendavid	1.16	double mass0 = 0.0;
151
152			FreeTrajectoryState initialState2 = tsTransform.initialFreeState(s2.track(),&magneticField);
153			helixIntersector.calculate(initialState1, initialState2);
154			if (helixIntersector.status()) {
155			dZ0 = fabs(helixIntersector.points().first.z() - helixIntersector.points().second.z());
156			dR0 = helixIntersector.crossingPoint().perp();
157
158			GlobalVector v1, v2;
159			v1 = helixIntersector.trajectoryParameters().first.momentum();
160			v2 = helixIntersector.trajectoryParameters().second.momentum();
161
162			double e1 = sqrt(v1.mag2()+s1.mass()*s1.mass());
163			double x1 = v1.x();
164			double y1 = v1.y();
165			double z1 = v1.z();
166
167			double e2 = sqrt(v2.mag2()+s2.mass()*s2.mass());
168			double x2 = v2.x();
169			double y2 = v2.y();
170			double z2 = v2.z();
171
172			FourVector sum(x1+x2, y1+y2, z1+z2, e1+e2);
173
174			mass0 = sqrt(sum.M2());
175			}
176
177
178	loizides	1.12	// Basic cuts on helix intersection
179	mrudolph	1.3	if(mass0 > massMax_ \|\| mass0<massMin_ \|\| fabs(dZ0) > dZMax_ \|\| dR0 < rhoMin_) continue;
180
181			// -------------------------------------------------------------------------------------------
182			// Do vertex fit for all pairs
183	loizides	1.7	// -------------------------------------------------------------------------------------------
184	paus	1.10	mithep::MultiVertexFitterD fit;
185	bendavid	1.15	fit.init(bfield); // Reset to the magnetic field from the event setup
186	mrudolph	1.3	fit.setChisqMax(100);
187	bendavid	1.21
188			fit.addTrack(trkPar1.pars(),trkPar1.cMat(),1,s1.mass(),MultiVertexFitterD::VERTEX_1);
189			fit.addTrack(trkPar2.pars(),trkPar2.cMat(),2,s2.mass(),MultiVertexFitterD::VERTEX_1);
190	mrudolph	1.3
191
192			if (fit.fit()){
193
194			DecayPart *d = new DecayPart(oPid_,DecayPart::Fast);
195	bendavid	1.5
196			// Update temporarily some of the quantities (prob, chi2, nDoF, mass, lxy, pt, fourMomentum)
197			d->setProb(fit.prob());
198			d->setChi2(fit.chisq());
199			d->setNdof(fit.ndof());
200
201			FourVector p4Fitted(0.,0.,0.,0.);
202			p4Fitted += fit.getTrackP4(1);
203			p4Fitted += fit.getTrackP4(2);
204			d->setFourMomentum(p4Fitted);
205	paus	1.10	d->setPosition (fit.getVertex (MultiVertexFitterD::VERTEX_1));
206			d->setError (fit.getErrorMatrix(MultiVertexFitterD::VERTEX_1));
207	bendavid	1.5	float mass, massErr;
208			const int trksIds[2] = { 1, 2 };
209			mass = fit.getMass(2,trksIds,massErr);
210
211			ThreeVector p3Fitted(p4Fitted.px(), p4Fitted.py(), p4Fitted.pz());
212
213			//Get decay length in xy plane
214			float dl, dlErr;
215	paus	1.10	dl = fit.getDecayLength (MultiVertexFitterD::PRIMARY_VERTEX, MultiVertexFitterD::VERTEX_1,
216			p3Fitted, dlErr);
217	bendavid	1.5
218			//Get Z decay length
219			float dlz, dlzErr;
220	paus	1.10	dlz = fit.getZDecayLength(MultiVertexFitterD::PRIMARY_VERTEX, MultiVertexFitterD::VERTEX_1,
221			p3Fitted, dlzErr);
222	bendavid	1.5
223			//get impact parameter
224			float dxy, dxyErr;
225	paus	1.10	dxy = fit.getImpactPar (MultiVertexFitterD::PRIMARY_VERTEX, MultiVertexFitterD::VERTEX_1,
226			p3Fitted, dxyErr);
227	bendavid	1.5
228	bendavid	1.8	BasePartPtr ptr1(hStables1,i);
229			BasePartPtr ptr2(hStables2,j);
230
231	paus	1.10	StableData c1(fit.getTrackP4(1).px(),fit.getTrackP4(1).py(), fit.getTrackP4(1).pz(), ptr1);
232			StableData c2(fit.getTrackP4(2).px(),fit.getTrackP4(2).py(), fit.getTrackP4(2).pz(), ptr2);
233	bendavid	1.8
234	paus	1.10	const ThreeVector vtxPos = fit.getVertex(MultiVertexFitterD::VERTEX_1);
235			const ThreeVector trkMom1(fit.getTrackP4(1).px(),fit.getTrackP4(1).py(),
236			fit.getTrackP4(1).pz());
237			const ThreeVector trkMom2(fit.getTrackP4(2).px(),fit.getTrackP4(2).py(),
238			fit.getTrackP4(2).pz());
239	bendavid	1.8
240			//build corrected HitPattern for StableData, removing hits before the fit vertex
241	bendavid	1.11	if (useHitDropper_) {
242	bendavid	1.22	std::pair<reco::HitPattern,uint> hits1 = dropper->CorrectedHitsAOD(s1.track(),
243	loizides	1.12	vtxPos,
244			trkMom1,
245			dlErr,
246			dlzErr);
247	bendavid	1.22	std::pair<reco::HitPattern,uint> hits2 = dropper->CorrectedHitsAOD(s2.track(),
248	loizides	1.12	vtxPos,
249			trkMom2,
250			dlErr,
251			dlzErr);
252	bendavid	1.11
253	bendavid	1.22	c1.SetHits(hits1.first);
254			c2.SetHits(hits2.first);
255	bendavid	1.11	c1.SetHitsFilled();
256			c2.SetHitsFilled();
257	bendavid	1.22	c1.SetNWrongHits(hits1.second);
258			c2.SetNWrongHits(hits2.second);
259	bendavid	1.20
260			reco::HitPattern sharedHits = dropper->SharedHits(s1.track(),s2.track());
261			d->setSharedHits(sharedHits);
262	bendavid	1.11	}
263	bendavid	1.8
264	paus	1.10	d->addStableChild (c1);
265			d->addStableChild (c2);
266			d->setFittedMass (mass);
267	bendavid	1.5	d->setFittedMassError(massErr);
268	paus	1.10	d->setLxy (dl);
269			d->setLxyError (dlErr);
270			d->setLxyToPv (dl);
271			d->setLxyToPvError (dlErr);
272			d->setLz (dlz);
273			d->setLzError (dlzErr);
274			d->setLzToPv (dlz);
275			d->setLzToPvError (dlzErr);
276			d->setDxy (dxy);
277			d->setDxyError (dxyErr);
278			d->setDxyToPv (dxy);
279			d->setDxyToPvError (dxyErr);
280
281	loizides	1.12	// put the result into our collection
282	mrudolph	1.3	pD->push_back(*d);
283	loizides	1.13	delete d;
284	mrudolph	1.3	} //done processing fit
285	loizides	1.7	} //end j loop
286	mrudolph	1.3	} //end i loop
287
288			// -----------------------------------------------------------------------------------------------
289			// Write the collection even if it is empty
290			// -----------------------------------------------------------------------------------------------
291	loizides	1.7	if (0)
292			cout << " V2SS::produce - " << pD->size() << " entries collection created -"
293			<< " (Pid: " << oPid_ << ")\n";
294	mrudolph	1.3	evt.put(pD);
295			}
296
297	loizides	1.12	// define this as a plug-in
298	mrudolph	1.3	DEFINE_FWK_MODULE(ProducerV2SS);