Producers/src/VertexZProducer.cc

// $Id: VertexZProducer.cc,v 1.12 2009/11/29 10:19:06 loizides Exp $

#include "MitEdm/Producers/interface/VertexZProducer.h"
#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
#include "DataFormats/GeometryVector/interface/LocalPoint.h"
#include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
#include "DataFormats/SiPixelDetId/interface/PXFDetId.h"
#include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
#include "DataFormats/VertexReco/interface/VertexFwd.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "Geometry/CommonDetUnit/interface/GeomDet.h"
#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
#include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetType.h"
#include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetUnit.h"
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
#include "Geometry/TrackerTopology/interface/RectangularPixelTopology.h" 
#include <TString.h>

using namespace std;
using namespace edm;

#define CP(level) \
  if (level>=verbose_)

//--------------------------------------------------------------------------------------------------
VertexZProducer::VertexZProducer(const ParameterSet &parameters) :
  pixelName_(parameters.getUntrackedParameter<string>("pixelRecHits","siPixelRecHits")),
  dPhiVc_(parameters.getUntrackedParameter<double>("dPhiVertexCut",0.08)),
  dZVc_(parameters.getUntrackedParameter<double>("dZVertexCut",0.25)),
  verbose_(parameters.getUntrackedParameter<int>("verbose",3)),
  pixLayers_(parameters.getUntrackedParameter<int>("pixLayerCombinations",12)),
  doClusVertex_(parameters.getUntrackedParameter<bool>("doClusVertex",false)),
  useRecHitQ_(parameters.getUntrackedParameter<bool>("useRecHitQualityWord",false)),
  usePixelQ_(parameters.getUntrackedParameter<bool>("usePixelQualityWord",true)),
  tgeo_(0)
{
  // Constructor.

  produces<reco::VertexCollection>();

  if ((pixLayers_!=12) && (pixLayers_!=13) && (pixLayers_!=23)) {
    print(2,Form("Value for pixLayerCombinations must be one of 12, 13, or 23. "
                 "Got %d, set value to 12", pixLayers_));
    pixLayers_ = 12;
  }
}

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

//--------------------------------------------------------------------------------------------------
void VertexZProducer::produce(Event &iEvent, const EventSetup &iSetup) 
{
  // Produce vertex collection.

  // get tracker geometry
  ESHandle<TrackerGeometry> trackerHandle;
  iSetup.get<TrackerDigiGeometryRecord>().get(trackerHandle);
  tgeo_ = trackerHandle.product();
  if (!tgeo_)
    print(3,"Could not obtain pointer to TrackerGeometry");

  // get pixels
  fillPixels(iEvent);

  // do vertex finding
  trackletV_.setinvalid();
  if (doClusVertex_) {
    vertexFromClusters(iEvent, pixLayers_);
  } else {
    trackletVertexUnbinned(iEvent, pixLayers_);
  }

  // publish result
  std::auto_ptr<reco::VertexCollection> vertices(new reco::VertexCollection);
  if (trackletV_.n()>0) {
    reco::Vertex::Error err;
    err(2,2) = trackletV_.zs()*trackletV_.zs();
    reco::Vertex ver(reco::Vertex::Point(0,0,trackletV_.z()), err, 
                     trackletV_.chi2(), trackletV_.n(), 0);
    vertices->push_back(ver);
  }
  iEvent.put(vertices);
}

//--------------------------------------------------------------------------------------------------
void VertexZProducer::fillPixels(const Event &iEvent) 
{
  // Fill pixel hit collections.

  bpix1_.clear();
  bpix2_.clear();
  bpix3_.clear();

  Handle<SiPixelRecHitCollection> hRecHits;
  if (!getProductSafe(pixelName_, hRecHits, iEvent)) {
    CP(2) print(2,Form("Can not obtain pixel hit collection with name %s", pixelName_.c_str()));
    return;
  }

  const SiPixelRecHitCollection *hits = hRecHits.product();
  for(SiPixelRecHitCollection::DataContainer::const_iterator hit = hits->data().begin(), 
        end = hits->data().end(); hit != end; ++hit) {

    if (!hit->isValid())
      continue;

    if (useRecHitQ_) {
      if (hit->isOnEdge() || hit->hasBadPixels())
        continue;
    }

    DetId id(hit->geographicalId());
    if(id.subdetId() != int(PixelSubdetector::PixelBarrel))
      continue;

    const PixelGeomDetUnit *pgdu = static_cast<const PixelGeomDetUnit*>(tgeo_->idToDet(id));

    if (usePixelQ_) {
      const RectangularPixelTopology *pixTopo = 
        static_cast<const RectangularPixelTopology*>(&(pgdu->specificTopology()));
      vector<SiPixelCluster::Pixel> pixels(hit->cluster()->pixels());
      bool pixelOnEdge = false;
      for(std::vector<SiPixelCluster::Pixel>::const_iterator pixel = pixels.begin(); 
          pixel != pixels.end(); ++pixel) {
        int pixelX = pixel->x;
        int pixelY = pixel->y;
        if(pixTopo->isItEdgePixelInX(pixelX) || pixTopo->isItEdgePixelInY(pixelY)) {
          pixelOnEdge = true;
          break;
        }
      }
      if (pixelOnEdge)
        continue;
    }

    LocalPoint lpos = LocalPoint(hit->localPosition().x(),
                                 hit->localPosition().y(),
                                 hit->localPosition().z());
    GlobalPoint gpos = pgdu->toGlobal(lpos);
    double adc   = hit->cluster()->charge()/135.;
    double sizex = hit->cluster()->sizeX();
    double sizey = hit->cluster()->sizeY();

    Pixel pix(gpos, adc, sizex, sizey);

    PXBDetId pid(id);
    int layer = pid.layer();

    if (layer==1) {
      bpix1_.push_back(pix);     
    } else if (layer==2) {
      bpix2_.push_back(pix);     
    } else {
      bpix3_.push_back(pix);     
    }
  }
}

//--------------------------------------------------------------------------------------------------
void VertexZProducer::reallyPrint(int level, const char *msg)
{
  // Print out information dependent on level and verbosity.

  if (level==0) {
    printf("VertexZProducer: %s\n", msg); 
  } else if (level==1) {
    LogWarning("VertexZProducer") << msg << std::endl;
  } else if (level==2) {
    LogError("VertexZProducer") << msg << std::endl;
  } else if (level==3) {
    LogError("VertexZProducer") << msg << std::endl;
    throw edm::Exception(errors::Configuration, "VertexZProducer\n") << msg << std::endl;
  }
}

//--------------------------------------------------------------------------------------------------
void VertexZProducer::trackletVertexUnbinned(const Event &iEvent, int which)
{
  // Estimate tracklet based z vertex.

  if (which==12) {
    trackletVertexUnbinned(bpix1_,bpix2_,trackletV_);
  } else if (which==13) {
    trackletVertexUnbinned(bpix1_,bpix3_,trackletV_);
  } else {
    trackletVertexUnbinned(bpix2_,bpix3_,trackletV_);
  }
}

//--------------------------------------------------------------------------------------------------
void VertexZProducer::trackletVertexUnbinned(std::vector<Pixel> &pix1, 
                                             std::vector<Pixel> &pix2,
                                             Vertex &vtx)
{
  // Calculate tracklet based z vertex position. 
  // At first build zvertex candidates from tracklet prototypes,
  // then group zvertex candidates and calculate mean position
  // from most likely cluster.

  vector<double> zvCands;
  zvCands.reserve(pix1.size()*pix2.size());

  // build candidates
  for(size_t i = 0; i<pix1.size(); ++i) {
    const Pixel &p1(pix1.at(i));
    const double r12 = p1.x()*p1.x()+p1.y()*p1.y();
    for(size_t j = 0; j<pix2.size(); ++j) {
      const Pixel &p2(pix2.at(j));
      if (TMath::Abs(Geom::deltaPhi(p1.phi(),p2.phi()))>dPhiVc_)
        continue;
      const double r22 = p2.x()*p2.x()+p2.y()*p2.y();
      const double vz = p1.z() - (p2.z()-p1.z())/(TMath::Sqrt(r22/r12)-1);
      if (TMath::IsNaN(vz))
        continue;
      if (TMath::Abs(vz)>25)
        continue;
      zvCands.push_back(vz);
    }
  }

  // sort cluster candidates
  sort(zvCands.begin(),zvCands.end());

  int    mcl=0;
  double ms2=1e10;
  double mzv=1e10;

  // cluster candidates and calculate mean z position
  for(size_t i = 0; i<zvCands.size(); ++i) {
    double z1 = zvCands.at(i);
    int    ncl   = 0;
    double mean  = 0;
    double mean2 = 0;
    for(size_t j = i; j<zvCands.size(); ++j) {
      double z2 = zvCands.at(j);
      if (TMath::Abs(z1-z2)>dZVc_)
        break;
      ++ncl;
      mean += z2;
      mean2 += z2*z2;
    }      
    if (ncl>0) {
      mean /= ncl;
      mean2 /= ncl;
    }
    double_t s2 = mean*mean - mean2;

    if ((ncl<mcl) || (ncl==mcl && s2>ms2))
      continue;

    mzv = mean;
    ms2 = s2;
    mcl = ncl;
  }

  // set the vertex
  vtx.set(mcl, mzv, ms2);
}

//--------------------------------------------------------------------------------------------------
void VertexZProducer::vertexFromClusters(const edm::Event &iEvent, int which)
{
  // Vertex from cluster positions.

  if (which==12) {
    vertexFromClusters(bpix1_,trackletV_);
  } else if (which==13) {
    vertexFromClusters(bpix3_,trackletV_);
  } else {
    vertexFromClusters(bpix2_,trackletV_);
  }
}

//--------------------------------------------------------------------------------------------------
void VertexZProducer::vertexFromClusters(const std::vector<Pixel> &pix, Vertex &vtx)
{
  // Estimate z vertex position from clusters.

  double chi_max = 1e+9;
  double z_best  = -999;
  int nhits_max  = 0;

  for(double z0 = -15.9; z0 <= 15.95; z0 += 0.1) {
    int nhits  = 0;
    double chi = 0;
    for(size_t i=0; i<pix.size(); ++i) {
      const Pixel &p = pix.at(i);

      // predicted cluster width in y direction
      double pval = 2*TMath::Abs(p.z()-z0)/p.rho() + 0.5; // FIXME
      double chitest = TMath::Abs(pval - p.sizey());
      if(chitest <= 1.) { 
        chi += chitest;
        ++nhits;
      }
    }

    if(nhits <= 0)
      continue;

    if(nhits < nhits_max)
      continue;

    if ((nhits > nhits_max) || (chi < chi_max)) { 
      z_best    = z0; 
      nhits_max = nhits; 
      chi_max   = chi; 
    }
  }

  // set the vertex
  vtx.set(nhits_max, z_best, 0.6);
  vtx.setchi2(chi_max);
}

// define this as a plug-in
DEFINE_FWK_MODULE(VertexZProducer);

Revision:	1.1
Committed:	Mon Nov 30 10:11:41 2009 UTC (15 years, 5 months ago) by loizides
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_020d, TMit_020d, Mit_020c, 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, Mit_014a, Mit_014, Mit_014pre3, Mit_014pre2, Mit_014pre1, Mit_013d, Mit_013c, Mit_013b, Mit_013a, Mit_013, Mit_013pre1, Mit_012i, Mit_012h, Mit_012g, Mit_012f, Mit_012e, Mit_012d, Mit_012c
Log Message:	Added vertex z producer.
#	User	Rev	Content
1	loizides	1.1	// $Id: VertexZProducer.cc,v 1.12 2009/11/29 10:19:06 loizides Exp $
2
3			#include "MitEdm/Producers/interface/VertexZProducer.h"
4			#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
5			#include "DataFormats/GeometryVector/interface/LocalPoint.h"
6			#include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
7			#include "DataFormats/SiPixelDetId/interface/PXFDetId.h"
8			#include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
9			#include "DataFormats/VertexReco/interface/VertexFwd.h"
10			#include "DataFormats/VertexReco/interface/Vertex.h"
11			#include "FWCore/Framework/interface/Event.h"
12			#include "FWCore/Framework/interface/EventSetup.h"
13			#include "FWCore/Framework/interface/ESHandle.h"
14			#include "FWCore/Framework/interface/MakerMacros.h"
15			#include "FWCore/ParameterSet/interface/ParameterSet.h"
16			#include "FWCore/MessageLogger/interface/MessageLogger.h"
17			#include "FWCore/ServiceRegistry/interface/Service.h"
18			#include "Geometry/CommonDetUnit/interface/GeomDet.h"
19			#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
20			#include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetType.h"
21			#include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetUnit.h"
22			#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
23			#include "Geometry/TrackerTopology/interface/RectangularPixelTopology.h"
24			#include <TString.h>
25
26			using namespace std;
27			using namespace edm;
28
29			#define CP(level) \
30			if (level>=verbose_)
31
32			//--------------------------------------------------------------------------------------------------
33			VertexZProducer::VertexZProducer(const ParameterSet &parameters) :
34			pixelName_(parameters.getUntrackedParameter<string>("pixelRecHits","siPixelRecHits")),
35			dPhiVc_(parameters.getUntrackedParameter<double>("dPhiVertexCut",0.08)),
36			dZVc_(parameters.getUntrackedParameter<double>("dZVertexCut",0.25)),
37			verbose_(parameters.getUntrackedParameter<int>("verbose",3)),
38			pixLayers_(parameters.getUntrackedParameter<int>("pixLayerCombinations",12)),
39			doClusVertex_(parameters.getUntrackedParameter<bool>("doClusVertex",false)),
40			useRecHitQ_(parameters.getUntrackedParameter<bool>("useRecHitQualityWord",false)),
41			usePixelQ_(parameters.getUntrackedParameter<bool>("usePixelQualityWord",true)),
42			tgeo_(0)
43			{
44			// Constructor.
45
46			produces<reco::VertexCollection>();
47
48			if ((pixLayers_!=12) && (pixLayers_!=13) && (pixLayers_!=23)) {
49			print(2,Form("Value for pixLayerCombinations must be one of 12, 13, or 23. "
50			"Got %d, set value to 12", pixLayers_));
51			pixLayers_ = 12;
52			}
53			}
54
55			//--------------------------------------------------------------------------------------------------
56			VertexZProducer::~VertexZProducer()
57			{
58			// Destructor.
59			}
60
61			//--------------------------------------------------------------------------------------------------
62			void VertexZProducer::produce(Event &iEvent, const EventSetup &iSetup)
63			{
64			// Produce vertex collection.
65
66			// get tracker geometry
67			ESHandle<TrackerGeometry> trackerHandle;
68			iSetup.get<TrackerDigiGeometryRecord>().get(trackerHandle);
69			tgeo_ = trackerHandle.product();
70			if (!tgeo_)
71			print(3,"Could not obtain pointer to TrackerGeometry");
72
73			// get pixels
74			fillPixels(iEvent);
75
76			// do vertex finding
77			trackletV_.setinvalid();
78			if (doClusVertex_) {
79			vertexFromClusters(iEvent, pixLayers_);
80			} else {
81			trackletVertexUnbinned(iEvent, pixLayers_);
82			}
83
84			// publish result
85			std::auto_ptr<reco::VertexCollection> vertices(new reco::VertexCollection);
86			if (trackletV_.n()>0) {
87			reco::Vertex::Error err;
88			err(2,2) = trackletV_.zs()*trackletV_.zs();
89			reco::Vertex ver(reco::Vertex::Point(0,0,trackletV_.z()), err,
90			trackletV_.chi2(), trackletV_.n(), 0);
91			vertices->push_back(ver);
92			}
93			iEvent.put(vertices);
94			}
95
96			//--------------------------------------------------------------------------------------------------
97			void VertexZProducer::fillPixels(const Event &iEvent)
98			{
99			// Fill pixel hit collections.
100
101			bpix1_.clear();
102			bpix2_.clear();
103			bpix3_.clear();
104
105			Handle<SiPixelRecHitCollection> hRecHits;
106			if (!getProductSafe(pixelName_, hRecHits, iEvent)) {
107			CP(2) print(2,Form("Can not obtain pixel hit collection with name %s", pixelName_.c_str()));
108			return;
109			}
110
111			const SiPixelRecHitCollection *hits = hRecHits.product();
112			for(SiPixelRecHitCollection::DataContainer::const_iterator hit = hits->data().begin(),
113			end = hits->data().end(); hit != end; ++hit) {
114
115			if (!hit->isValid())
116			continue;
117
118			if (useRecHitQ_) {
119			if (hit->isOnEdge() \|\| hit->hasBadPixels())
120			continue;
121			}
122
123			DetId id(hit->geographicalId());
124			if(id.subdetId() != int(PixelSubdetector::PixelBarrel))
125			continue;
126
127			const PixelGeomDetUnit pgdu = static_cast<const PixelGeomDetUnit>(tgeo_->idToDet(id));
128
129			if (usePixelQ_) {
130			const RectangularPixelTopology *pixTopo =
131			static_cast<const RectangularPixelTopology*>(&(pgdu->specificTopology()));
132			vector<SiPixelCluster::Pixel> pixels(hit->cluster()->pixels());
133			bool pixelOnEdge = false;
134			for(std::vector<SiPixelCluster::Pixel>::const_iterator pixel = pixels.begin();
135			pixel != pixels.end(); ++pixel) {
136			int pixelX = pixel->x;
137			int pixelY = pixel->y;
138			if(pixTopo->isItEdgePixelInX(pixelX) \|\| pixTopo->isItEdgePixelInY(pixelY)) {
139			pixelOnEdge = true;
140			break;
141			}
142			}
143			if (pixelOnEdge)
144			continue;
145			}
146
147			LocalPoint lpos = LocalPoint(hit->localPosition().x(),
148			hit->localPosition().y(),
149			hit->localPosition().z());
150			GlobalPoint gpos = pgdu->toGlobal(lpos);
151			double adc = hit->cluster()->charge()/135.;
152			double sizex = hit->cluster()->sizeX();
153			double sizey = hit->cluster()->sizeY();
154
155			Pixel pix(gpos, adc, sizex, sizey);
156
157			PXBDetId pid(id);
158			int layer = pid.layer();
159
160			if (layer==1) {
161			bpix1_.push_back(pix);
162			} else if (layer==2) {
163			bpix2_.push_back(pix);
164			} else {
165			bpix3_.push_back(pix);
166			}
167			}
168			}
169
170			//--------------------------------------------------------------------------------------------------
171			void VertexZProducer::reallyPrint(int level, const char *msg)
172			{
173			// Print out information dependent on level and verbosity.
174
175			if (level==0) {
176			printf("VertexZProducer: %s\n", msg);
177			} else if (level==1) {
178			LogWarning("VertexZProducer") << msg << std::endl;
179			} else if (level==2) {
180			LogError("VertexZProducer") << msg << std::endl;
181			} else if (level==3) {
182			LogError("VertexZProducer") << msg << std::endl;
183			throw edm::Exception(errors::Configuration, "VertexZProducer\n") << msg << std::endl;
184			}
185			}
186
187			//--------------------------------------------------------------------------------------------------
188			void VertexZProducer::trackletVertexUnbinned(const Event &iEvent, int which)
189			{
190			// Estimate tracklet based z vertex.
191
192			if (which==12) {
193			trackletVertexUnbinned(bpix1_,bpix2_,trackletV_);
194			} else if (which==13) {
195			trackletVertexUnbinned(bpix1_,bpix3_,trackletV_);
196			} else {
197			trackletVertexUnbinned(bpix2_,bpix3_,trackletV_);
198			}
199			}
200
201			//--------------------------------------------------------------------------------------------------
202			void VertexZProducer::trackletVertexUnbinned(std::vector<Pixel> &pix1,
203			std::vector<Pixel> &pix2,
204			Vertex &vtx)
205			{
206			// Calculate tracklet based z vertex position.
207			// At first build zvertex candidates from tracklet prototypes,
208			// then group zvertex candidates and calculate mean position
209			// from most likely cluster.
210
211			vector<double> zvCands;
212			zvCands.reserve(pix1.size()*pix2.size());
213
214			// build candidates
215			for(size_t i = 0; i<pix1.size(); ++i) {
216			const Pixel &p1(pix1.at(i));
217			const double r12 = p1.x()p1.x()+p1.y()p1.y();
218			for(size_t j = 0; j<pix2.size(); ++j) {
219			const Pixel &p2(pix2.at(j));
220			if (TMath::Abs(Geom::deltaPhi(p1.phi(),p2.phi()))>dPhiVc_)
221			continue;
222			const double r22 = p2.x()p2.x()+p2.y()p2.y();
223			const double vz = p1.z() - (p2.z()-p1.z())/(TMath::Sqrt(r22/r12)-1);
224			if (TMath::IsNaN(vz))
225			continue;
226			if (TMath::Abs(vz)>25)
227			continue;
228			zvCands.push_back(vz);
229			}
230			}
231
232			// sort cluster candidates
233			sort(zvCands.begin(),zvCands.end());
234
235			int mcl=0;
236			double ms2=1e10;
237			double mzv=1e10;
238
239			// cluster candidates and calculate mean z position
240			for(size_t i = 0; i<zvCands.size(); ++i) {
241			double z1 = zvCands.at(i);
242			int ncl = 0;
243			double mean = 0;
244			double mean2 = 0;
245			for(size_t j = i; j<zvCands.size(); ++j) {
246			double z2 = zvCands.at(j);
247			if (TMath::Abs(z1-z2)>dZVc_)
248			break;
249			++ncl;
250			mean += z2;
251			mean2 += z2*z2;
252			}
253			if (ncl>0) {
254			mean /= ncl;
255			mean2 /= ncl;
256			}
257			double_t s2 = mean*mean - mean2;
258
259			if ((ncl<mcl) \|\| (ncl==mcl && s2>ms2))
260			continue;
261
262			mzv = mean;
263			ms2 = s2;
264			mcl = ncl;
265			}
266
267			// set the vertex
268			vtx.set(mcl, mzv, ms2);
269			}
270
271			//--------------------------------------------------------------------------------------------------
272			void VertexZProducer::vertexFromClusters(const edm::Event &iEvent, int which)
273			{
274			// Vertex from cluster positions.
275
276			if (which==12) {
277			vertexFromClusters(bpix1_,trackletV_);
278			} else if (which==13) {
279			vertexFromClusters(bpix3_,trackletV_);
280			} else {
281			vertexFromClusters(bpix2_,trackletV_);
282			}
283			}
284
285			//--------------------------------------------------------------------------------------------------
286			void VertexZProducer::vertexFromClusters(const std::vector<Pixel> &pix, Vertex &vtx)
287			{
288			// Estimate z vertex position from clusters.
289
290			double chi_max = 1e+9;
291			double z_best = -999;
292			int nhits_max = 0;
293
294			for(double z0 = -15.9; z0 <= 15.95; z0 += 0.1) {
295			int nhits = 0;
296			double chi = 0;
297			for(size_t i=0; i<pix.size(); ++i) {
298			const Pixel &p = pix.at(i);
299
300			// predicted cluster width in y direction
301			double pval = 2*TMath::Abs(p.z()-z0)/p.rho() + 0.5; // FIXME
302			double chitest = TMath::Abs(pval - p.sizey());
303			if(chitest <= 1.) {
304			chi += chitest;
305			++nhits;
306			}
307			}
308
309			if(nhits <= 0)
310			continue;
311
312			if(nhits < nhits_max)
313			continue;
314
315			if ((nhits > nhits_max) \|\| (chi < chi_max)) {
316			z_best = z0;
317			nhits_max = nhits;
318			chi_max = chi;
319			}
320			}
321
322			// set the vertex
323			vtx.set(nhits_max, z_best, 0.6);
324			vtx.setchi2(chi_max);
325			}
326
327			// define this as a plug-in
328			DEFINE_FWK_MODULE(VertexZProducer);
329