Producers/src/ProducerV2SS.cc

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

#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"
#include "MitEdm/VertexFitInterface/interface/HisInterface.h"
#include "MitEdm/Producers/interface/ProducerV2SS.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))
{
}

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

//--------------------------------------------------------------------------------------------------
void ProducerV2SS::produce(Event &evt, const EventSetup &setup)
{
  // Create 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();
  }

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

  //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);
   
    UInt_t j;
    if (iStables1_ == iStables2_)
      j = i+1; 
    else
      j = 0;
    
    for (; j<pS2->size(); ++j) {
      const StablePart &s2 = pS2->at(j);

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

      //Do fast helix fit to check if there's any hope
      const reco::Track * t1 = s1.track();
      const reco::Track * t2 = s2.track();
      HisInterface            hisInt(t1,t2);
      const mithep::HelixIntersector *his = hisInt.hISector();
      double dZ0 = -999;
      double dR0 = -999;
      dZ0 = his->IntersectionI(0).DeltaZ();
      dR0 = sqrt(his->IntersectionI(0).Location().X()*his->IntersectionI(0).Location().X()+
                 his->IntersectionI(0).Location().Y()*his->IntersectionI(0).Location().Y());

      TVector3     v1, v2;
      v1 = his->IntersectionI(0).TrackParamsI(0).Momentum();
      v2 = his->IntersectionI(0).TrackParamsI(1).Momentum();

      double e1 = sqrt(v1.Mag()*v1.Mag() +s1.mass()*s1.mass());
      double x1 = v1.X();
      double y1 = v1.Y();
      double z1 = v1.Z();

      double e2 = sqrt(v2.Mag()*v2.Mag()+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);

      double 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(3.8); // Reset to the MC magnetic field of 3.8 Tesla
      fit.setChisqMax(100);
      MvfInterface fitInt(&fit);
      fitInt.addTrack(s1.track(),1,s1.mass(),mithep::MultiVertexFitterD::VERTEX_1);
      fitInt.addTrack(s2.track(),2,s2.mass(),mithep::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_) {
          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);
      }  //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.11
Committed:	Tue Mar 3 21:31:08 2009 UTC (16 years, 2 months ago) by bendavid
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_008pre2, Mit_008pre1
Changes since 1.10:	+17 -11 lines
Log Message:	Make V producer runnable in Aod by making HitDropper usage configurable, and changing the nofit python file a bit
#	User	Rev	Content
1	bendavid	1.11	// $Id: ProducerV2SS.cc,v 1.10 2008/11/13 17:08:31 paus 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.8	#include "MitEdm/Producers/interface/HitDropperRecord.h"
7			#include "MitEdm/Producers/interface/HitDropper.h"
8	mrudolph	1.3	#include "MitEdm/DataFormats/interface/Types.h"
9	mrudolph	1.4	#include "MitEdm/DataFormats/interface/Collections.h"
10	mrudolph	1.3	#include "MitEdm/DataFormats/interface/DecayPart.h"
11			#include "MitEdm/DataFormats/interface/StablePart.h"
12	bendavid	1.5	#include "MitEdm/DataFormats/interface/StableData.h"
13	mrudolph	1.3	#include "MitEdm/VertexFitInterface/interface/MvfInterface.h"
14			#include "MitEdm/VertexFitInterface/interface/HisInterface.h"
15			#include "MitEdm/Producers/interface/ProducerV2SS.h"
16
17			using namespace std;
18			using namespace edm;
19			using namespace mitedm;
20	bendavid	1.5	using namespace mithep;
21	mrudolph	1.3
22			//--------------------------------------------------------------------------------------------------
23			ProducerV2SS::ProducerV2SS(const ParameterSet& cfg) :
24			ProducerD2SS(cfg),
25	paus	1.10	rhoMin_ (cfg.getUntrackedParameter<double>("minRadius", 0.0)),
26			massMin_ (cfg.getUntrackedParameter<double>("minMass", 0.0)),
27			massMax_ (cfg.getUntrackedParameter<double>("maxMass", 3.0)),
28	bendavid	1.11	dZMax_ (cfg.getUntrackedParameter<double>("maxZDistance",5.0)),
29			useHitDropper_(cfg.getUntrackedParameter<bool>("useHitDropper",true))
30	mrudolph	1.3	{
31			}
32
33			//--------------------------------------------------------------------------------------------------
34			ProducerV2SS::~ProducerV2SS()
35			{
36			}
37
38			//--------------------------------------------------------------------------------------------------
39			void ProducerV2SS::produce(Event &evt, const EventSetup &setup)
40			{
41			// Create the output collection
42			auto_ptr<DecayPartCol> pD(new DecayPartCol());
43
44			// First input collection
45			Handle<StablePartCol> hStables1;
46			Handle<StablePartCol> hStables2;
47			if (!GetProduct(iStables1_, hStables1, evt) ) {
48			cout << "Couldn't get in collection in Producer V2SS" << endl;
49			evt.put(pD);
50			return;
51			}
52	loizides	1.7	const StablePartCol *pS1 = hStables1.product();
53	mrudolph	1.3
54	loizides	1.7	// Second input collection
55	mrudolph	1.3	if(!GetProduct(iStables2_, hStables2, evt)) {
56			cout << "Couldn't get in collection in Producer V2SS" << endl;
57			evt.put(pD);
58			return;
59			}
60			const StablePartCol *pS2 = hStables2.product();
61	bendavid	1.8
62			//get hit dropper
63			ESHandle<HitDropper> hDropper;
64	bendavid	1.11	const HitDropper *dropper = 0;
65			if (useHitDropper_) {
66			setup.get<HitDropperRecord>().get("HitDropper",hDropper);
67			dropper = hDropper.product();
68			}
69	mrudolph	1.3
70			// -----------------------------------------------------------------------------------------------
71			// Simple double loop
72			// -----------------------------------------------------------------------------------------------
73	loizides	1.7	if (0)
74			cout << "Starting V finder loop" << endl;
75	mrudolph	1.3
76	loizides	1.7	//sX_y: X= pion or proton collection.
77			//i, j = 2 loop particles.
78			//ex.: s1_i and s2_i are same particle as pion and proton
79	mrudolph	1.3	for (UInt_t i=0; i<pS1->size(); ++i) {
80			const StablePart &s1 = pS1->at(i);
81
82			UInt_t j;
83			if (iStables1_ == iStables2_)
84			j = i+1;
85			else
86			j = 0;
87
88			for (; j<pS2->size(); ++j) {
89			const StablePart &s2 = pS2->at(j);
90
91			if(s1.charge() + s2.charge() != 0) continue;
92
93			//Do fast helix fit to check if there's any hope
94			const reco::Track * t1 = s1.track();
95			const reco::Track * t2 = s2.track();
96			HisInterface hisInt(t1,t2);
97			const mithep::HelixIntersector *his = hisInt.hISector();
98			double dZ0 = -999;
99			double dR0 = -999;
100			dZ0 = his->IntersectionI(0).DeltaZ();
101			dR0 = sqrt(his->IntersectionI(0).Location().X()*his->IntersectionI(0).Location().X()+
102			his->IntersectionI(0).Location().Y()*his->IntersectionI(0).Location().Y());
103
104			TVector3 v1, v2;
105			v1 = his->IntersectionI(0).TrackParamsI(0).Momentum();
106			v2 = his->IntersectionI(0).TrackParamsI(1).Momentum();
107
108			double e1 = sqrt(v1.Mag()v1.Mag() +s1.mass()s1.mass());
109			double x1 = v1.X();
110			double y1 = v1.Y();
111			double z1 = v1.Z();
112
113			double e2 = sqrt(v2.Mag()v2.Mag()+s2.mass()s2.mass());
114			double x2 = v2.X();
115			double y2 = v2.Y();
116			double z2 = v2.Z();
117
118			FourVector sum(x1+x2, y1+y2, z1+z2, e1+e2);
119
120			double mass0 = sqrt(sum.M2());
121
122			//Basic cuts on helix intersection
123			if(mass0 > massMax_ \|\| mass0<massMin_ \|\| fabs(dZ0) > dZMax_ \|\| dR0 < rhoMin_) continue;
124
125			// -------------------------------------------------------------------------------------------
126			// Do vertex fit for all pairs
127	loizides	1.7	// -------------------------------------------------------------------------------------------
128	paus	1.10	mithep::MultiVertexFitterD fit;
129	loizides	1.7	fit.init(3.8); // Reset to the MC magnetic field of 3.8 Tesla
130	mrudolph	1.3	fit.setChisqMax(100);
131			MvfInterface fitInt(&fit);
132	paus	1.10	fitInt.addTrack(s1.track(),1,s1.mass(),mithep::MultiVertexFitterD::VERTEX_1);
133			fitInt.addTrack(s2.track(),2,s2.mass(),mithep::MultiVertexFitterD::VERTEX_1);
134	mrudolph	1.3
135
136			if (fit.fit()){
137
138			DecayPart *d = new DecayPart(oPid_,DecayPart::Fast);
139	bendavid	1.5
140			// Update temporarily some of the quantities (prob, chi2, nDoF, mass, lxy, pt, fourMomentum)
141			d->setProb(fit.prob());
142			d->setChi2(fit.chisq());
143			d->setNdof(fit.ndof());
144
145			FourVector p4Fitted(0.,0.,0.,0.);
146			p4Fitted += fit.getTrackP4(1);
147			p4Fitted += fit.getTrackP4(2);
148			d->setFourMomentum(p4Fitted);
149	paus	1.10	d->setPosition (fit.getVertex (MultiVertexFitterD::VERTEX_1));
150			d->setError (fit.getErrorMatrix(MultiVertexFitterD::VERTEX_1));
151	bendavid	1.5	float mass, massErr;
152			const int trksIds[2] = { 1, 2 };
153			mass = fit.getMass(2,trksIds,massErr);
154
155			ThreeVector p3Fitted(p4Fitted.px(), p4Fitted.py(), p4Fitted.pz());
156
157			//Get decay length in xy plane
158			float dl, dlErr;
159	paus	1.10	dl = fit.getDecayLength (MultiVertexFitterD::PRIMARY_VERTEX, MultiVertexFitterD::VERTEX_1,
160			p3Fitted, dlErr);
161	bendavid	1.5
162			//Get Z decay length
163			float dlz, dlzErr;
164	paus	1.10	dlz = fit.getZDecayLength(MultiVertexFitterD::PRIMARY_VERTEX, MultiVertexFitterD::VERTEX_1,
165			p3Fitted, dlzErr);
166	bendavid	1.5
167			//get impact parameter
168			float dxy, dxyErr;
169	paus	1.10	dxy = fit.getImpactPar (MultiVertexFitterD::PRIMARY_VERTEX, MultiVertexFitterD::VERTEX_1,
170			p3Fitted, dxyErr);
171	bendavid	1.5
172	bendavid	1.8	BasePartPtr ptr1(hStables1,i);
173			BasePartPtr ptr2(hStables2,j);
174
175	paus	1.10	StableData c1(fit.getTrackP4(1).px(),fit.getTrackP4(1).py(), fit.getTrackP4(1).pz(), ptr1);
176			StableData c2(fit.getTrackP4(2).px(),fit.getTrackP4(2).py(), fit.getTrackP4(2).pz(), ptr2);
177	bendavid	1.8
178	paus	1.10	const ThreeVector vtxPos = fit.getVertex(MultiVertexFitterD::VERTEX_1);
179			const ThreeVector trkMom1(fit.getTrackP4(1).px(),fit.getTrackP4(1).py(),
180			fit.getTrackP4(1).pz());
181			const ThreeVector trkMom2(fit.getTrackP4(2).px(),fit.getTrackP4(2).py(),
182			fit.getTrackP4(2).pz());
183	bendavid	1.8
184			//build corrected HitPattern for StableData, removing hits before the fit vertex
185	bendavid	1.11	if (useHitDropper_) {
186			reco::HitPattern hits1 = dropper->CorrectedHits(s1.track(), vtxPos, trkMom1, dlErr, dlzErr);
187			reco::HitPattern hits2 = dropper->CorrectedHits(s2.track(), vtxPos, trkMom2, dlErr, dlzErr);
188
189			c1.SetHits(hits1);
190			c2.SetHits(hits2);
191			c1.SetHitsFilled();
192			c2.SetHitsFilled();
193			}
194	bendavid	1.8
195	paus	1.10	d->addStableChild (c1);
196			d->addStableChild (c2);
197			d->setFittedMass (mass);
198	bendavid	1.5	d->setFittedMassError(massErr);
199	paus	1.10	d->setLxy (dl);
200			d->setLxyError (dlErr);
201			d->setLxyToPv (dl);
202			d->setLxyToPvError (dlErr);
203			d->setLz (dlz);
204			d->setLzError (dlzErr);
205			d->setLzToPv (dlz);
206			d->setLzToPvError (dlzErr);
207			d->setDxy (dxy);
208			d->setDxyError (dxyErr);
209			d->setDxyToPv (dxy);
210			d->setDxyToPvError (dxyErr);
211
212	mrudolph	1.3	// Put the result into our collection
213			pD->push_back(*d);
214			} //done processing fit
215	loizides	1.7	} //end j loop
216	mrudolph	1.3	} //end i loop
217
218			// -----------------------------------------------------------------------------------------------
219			// Write the collection even if it is empty
220			// -----------------------------------------------------------------------------------------------
221	loizides	1.7	if (0)
222			cout << " V2SS::produce - " << pD->size() << " entries collection created -"
223			<< " (Pid: " << oPid_ << ")\n";
224	mrudolph	1.3	evt.put(pD);
225			}
226
227			//define this as a plug-in
228			DEFINE_FWK_MODULE(ProducerV2SS);
229