Producers/src/ProducerEvtSelData.cc

// $Id: ProducerEvtSelData.cc,v 1.9 2009/12/11 21:30:37 edwenger Exp $

#include "MitEdm/Producers/interface/ProducerEvtSelData.h"
#include "MitEdm/DataFormats/interface/EvtSelData.h"
#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
#include "DataFormats/GeometryVector/interface/LocalPoint.h"
#include "DataFormats/HcalDetId/interface/HcalDetId.h"
#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
#include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
#include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/VertexReco/interface/VertexFwd.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
#include "Geometry/CommonDetUnit/interface/GeomDet.h"
#include "Geometry/Records/interface/CaloGeometryRecord.h"
#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
#include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetUnit.h"
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
#include "Geometry/TrackerTopology/interface/RectangularPixelTopology.h"
#include <TMath.h>

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

//--------------------------------------------------------------------------------------------------
ProducerEvtSelData::ProducerEvtSelData(const ParameterSet& cfg)
  : srcHF_(cfg.getUntrackedParameter<string>("hfRecHits","hfreco")),
    srcHBHE_(cfg.getUntrackedParameter<string>("hbheRecHits","hbhereco")),
    srcCastor_(cfg.getUntrackedParameter<string>("castorRecHits","castorreco")),
    srcZDC_(cfg.getUntrackedParameter<string>("zdcRecHits","zdcreco")),
    srcPixels_(cfg.getUntrackedParameter<string>("pixelRecHits","siPixelRecHits")),
    srcVertex_(cfg.getUntrackedParameter<string>("vertex",""))
{
  // Constructor.

  produces<EvtSelData>();
}

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

//--------------------------------------------------------------------------------------------------
void ProducerEvtSelData::produce(Event &evt, const EventSetup &setup)
{
  // Produce event selection data for this event.

  double eHfNeg       = 0;
  double eHfPos       = 0;
  double eHfPosTime   = 0;
  double eHfNegTime   = 0;
  double eHcalNeg     = 0;
  double eHcalPos     = 0;
  double eCaNeg       = 0;
  double eCaPos       = 0;
  double eCaPosTime   = 0;
  double eCaNegTime   = 0;
  double eZdNeg       = 0;
  double eZdPos       = 0;
  double eZdPosTime   = 0;
  double eZdNegTime   = 0;
  int    ePxbHits     = 0;
  int    ePxHits      = 0;
  double eClusVtxQual = 0;
  double eClusVtxDiff = 0;

  Handle<HFRecHitCollection> hfhits;
  try {
    evt.getByLabel(edm::InputTag(srcHF_),hfhits);
  } catch (...) {}  
  if (hfhits.isValid()) {
    for(size_t ihit = 0; ihit<hfhits->size(); ++ihit){
      const HFRecHit h = (*hfhits)[ihit];
      double energy = h.energy();
      double time = h.time();
      const HcalDetId id(h.id());
      if(energy<0) continue; 
      if (id.zside()<0) {
        eHfNeg     += energy;
        eHfNegTime += energy * time;
      } else {
        eHfPos     += energy;
        eHfPosTime += energy * time;
      }
    }
  }

  Handle<HBHERecHitCollection> hbhehits;
  try {
    evt.getByLabel(edm::InputTag(srcHBHE_),hbhehits);
  } catch (...) {}  
  if (hbhehits.isValid()) {
    for(size_t ihit = 0; ihit<hbhehits->size(); ++ihit){
      const HBHERecHit h = (*hbhehits)[ihit];
      double energy = h.energy();
      const HcalDetId id(h.id()); 
      if(energy<0) continue;
      if (id.zside()<0) {
        eHcalNeg   += energy;
      } else {
        eHcalPos   += energy;
      }
    }
  }

  Handle<CastorRecHitCollection> castorhits;
  try {
    evt.getByLabel(edm::InputTag(srcCastor_),castorhits);
  } catch (...) {}  
  if (castorhits.isValid()) {
    for(size_t ihit = 0; ihit<castorhits->size(); ++ihit){
      const CastorRecHit h = (*castorhits)[ihit];
      double energy = h.energy();
      double time = h.time();
      const HcalCastorDetId id(h.id()); 
      if(energy<0) continue;
      if (id.zside()<0) {
        eCaNeg     += energy;
        eCaNegTime += energy * time;
      } else {
        eCaPos     += energy;
        eCaPosTime += energy * time;
      }
    }
  }

  Handle<ZDCRecHitCollection> zdchits;
  try {
    evt.getByLabel(edm::InputTag(srcZDC_),zdchits);
  } catch (...) {}  
  if (zdchits.isValid()) {
    for(size_t ihit = 0; ihit<zdchits->size(); ++ihit){
      const ZDCRecHit h = (*zdchits)[ihit];
      double energy = h.energy();
      double time = h.time();
      const HcalZDCDetId id(h.id()); 
      if(energy<0) continue;
      if (id.zside()<0) {
        eZdNeg     += energy;
        eZdNegTime += energy * time;
      } else {
        eZdPos     += energy;
        eZdPosTime += energy * time;
      }
    }
  }

  if(eHfPos>0)
    eHfPosTime /= eHfPos;
  if(eHfNeg>0)
    eHfNegTime /= eHfNeg;
  if(eCaPos>0)
    eCaPosTime /= eCaPos;
  if(eCaNeg>0)
    eCaNegTime /= eCaNeg;
  if(eZdPos>0)
    eZdPosTime /= eZdPos;
  if(eZdNeg>0)
    eZdNegTime /= eZdNeg;
  
  Handle<SiPixelRecHitCollection> hRecHits;
  try {
    evt.getByLabel(edm::InputTag(srcPixels_),hRecHits);
  } catch (...) {}

  if (hRecHits.isValid()) {
    ESHandle<TrackerGeometry> trackerHandle;
    setup.get<TrackerDigiGeometryRecord>().get(trackerHandle);
    const TrackerGeometry *tgeo = trackerHandle.product();
    const SiPixelRecHitCollection *hits = hRecHits.product();

    vector<VertexHit> vhits(hits->dataSize());
    for(SiPixelRecHitCollection::DataContainer::const_iterator hit = hits->data().begin(), 
          end = hits->data().end(); hit != end; ++hit) {
      if (!hit->isValid())
        continue;
      ++ePxHits;
      DetId id(hit->geographicalId());
      if(id.subdetId() != int(PixelSubdetector::PixelBarrel))
        continue;
      ++ePxbHits;
      const PixelGeomDetUnit *pgdu = static_cast<const PixelGeomDetUnit*>(tgeo->idToDet(id));
      if (1) {
        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);
      VertexHit vh;
      vh.z = gpos.z(); 
      vh.r = gpos.perp(); 
      vh.w = hit->cluster()->sizeY();
      vhits.push_back(vh);
    }

    double zest = 0.0;
    if (srcVertex_.length()!=0) {
      edm::Handle<reco::VertexCollection> vertexCol;
      evt.getByLabel(srcVertex_,vertexCol);
      const reco::VertexCollection *vertices = vertexCol.product();
      unsigned int vtracks = 0;
      for(reco::VertexCollection::const_iterator
            vertex = vertices->begin(); vertex!= vertices->end(); ++vertex) {
        if(vertex->tracksSize()>vtracks) {
          vtracks = vertex->tracksSize();
          zest = vertex->z();
        }
      }
    } else {
      int nhits = 0, nhits_max = 0;
      double chi = 0, chi_max = 1e+9;
      for(double z0 = -15.9; z0 <= 15.95; z0 += 0.1) {
        nhits = getContainedHits(vhits, z0, chi);
        if(nhits > 0) {
          if(nhits >  nhits_max) { 
            chi_max = 1e+9; 
            nhits_max = nhits; 
          }
          if(nhits >= nhits_max) {
            if(chi < chi_max) { 
              chi_max = chi; 
              zest = z0; 
            }
          }
        }
      }
    }

    double chi = 0;
    int nbest=0, nminus=0, nplus=0;
    nbest = getContainedHits(vhits,zest,chi);
    nminus = getContainedHits(vhits,zest-10.,chi);
    nplus = getContainedHits(vhits,zest+10.,chi);

    eClusVtxDiff = nbest - (nminus+nplus)/2.;
    if ((nminus+nplus)> 0)
      eClusVtxQual = (2.0*nbest)/(nminus+nplus);  // A/B
    else if (nbest>0)
      eClusVtxQual = 1000.0;                      // A/0 (set to arbitrarily large number)
    else
      eClusVtxQual = 0;                           // 0/0 (already the default)

  }

  std::auto_ptr<EvtSelData> output(new EvtSelData(eHcalNeg, eHcalPos,
                                                  eHfNeg,eHfPos,eHfNegTime,eHfPosTime,
                                                  eCaNeg,eCaPos,eCaNegTime,eCaPosTime,
                                                  eZdNeg,eZdPos,eZdNegTime,eZdPosTime,
                                                  ePxbHits,ePxHits,eClusVtxQual,eClusVtxDiff));
  evt.put(output);
}

//--------------------------------------------------------------------------------------------------
int ProducerEvtSelData::getContainedHits(const std::vector<VertexHit> &hits, double z0, double &chi) 
{
  // Calculate number of hits contained in v-shaped window in cluster y-width vs. z-position.

  int n = 0;
  chi   = 0.;

  for(vector<VertexHit>::const_iterator hit = hits.begin(); hit!= hits.end(); hit++) {
    double p = 2 * fabs(hit->z - z0)/hit->r + 0.5; // FIXME
    if(TMath::Abs(p - hit->w) <= 1.) { 
      chi += fabs(p - hit->w);
      n++;
    }
  }
  return n;
}

//define this as a plug-in
DEFINE_FWK_MODULE(ProducerEvtSelData);
Revision:	1.10
Committed:	Sat Dec 12 13:37:30 2009 UTC (15 years, 5 months ago) by edwenger
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_012g, Mit_012f, Mit_012e
Changes since 1.9:	+2 -2 lines
Log Message:	truncate cluster shape cut at 150 pixel hits to avoid any interference with vertexing inefficiency at 0,0
#	User	Rev	Content
1	edwenger	1.10	// $Id: ProducerEvtSelData.cc,v 1.9 2009/12/11 21:30:37 edwenger Exp $
2	loizides	1.1
3			#include "MitEdm/Producers/interface/ProducerEvtSelData.h"
4			#include "MitEdm/DataFormats/interface/EvtSelData.h"
5			#include "DataFormats/Common/interface/Handle.h"
6	loizides	1.4	#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
7			#include "DataFormats/GeometryVector/interface/LocalPoint.h"
8			#include "DataFormats/HcalDetId/interface/HcalDetId.h"
9	loizides	1.1	#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
10	loizides	1.4	#include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
11			#include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
12			#include "DataFormats/VertexReco/interface/Vertex.h"
13			#include "DataFormats/VertexReco/interface/VertexFwd.h"
14	loizides	1.1	#include "FWCore/Framework/interface/ESHandle.h"
15			#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
16	loizides	1.4	#include "Geometry/CommonDetUnit/interface/GeomDet.h"
17	loizides	1.1	#include "Geometry/Records/interface/CaloGeometryRecord.h"
18	loizides	1.4	#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
19			#include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetUnit.h"
20	edwenger	1.2	#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
21			#include "Geometry/TrackerTopology/interface/RectangularPixelTopology.h"
22	loizides	1.4	#include <TMath.h>
23	edwenger	1.2
24	loizides	1.1	using namespace std;
25			using namespace edm;
26			using namespace mitedm;
27
28			//--------------------------------------------------------------------------------------------------
29			ProducerEvtSelData::ProducerEvtSelData(const ParameterSet& cfg)
30			: srcHF_(cfg.getUntrackedParameter<string>("hfRecHits","hfreco")),
31			srcHBHE_(cfg.getUntrackedParameter<string>("hbheRecHits","hbhereco")),
32			srcCastor_(cfg.getUntrackedParameter<string>("castorRecHits","castorreco")),
33	edwenger	1.2	srcZDC_(cfg.getUntrackedParameter<string>("zdcRecHits","zdcreco")),
34	loizides	1.4	srcPixels_(cfg.getUntrackedParameter<string>("pixelRecHits","siPixelRecHits")),
35			srcVertex_(cfg.getUntrackedParameter<string>("vertex",""))
36	loizides	1.1	{
37			// Constructor.
38
39			produces<EvtSelData>();
40			}
41
42			//--------------------------------------------------------------------------------------------------
43			ProducerEvtSelData::~ProducerEvtSelData()
44			{
45			// Destructor.
46			}
47
48			//--------------------------------------------------------------------------------------------------
49			void ProducerEvtSelData::produce(Event &evt, const EventSetup &setup)
50			{
51			// Produce event selection data for this event.
52
53	loizides	1.4	double eHfNeg = 0;
54			double eHfPos = 0;
55			double eHfPosTime = 0;
56			double eHfNegTime = 0;
57			double eHcalNeg = 0;
58			double eHcalPos = 0;
59			double eCaNeg = 0;
60			double eCaPos = 0;
61			double eCaPosTime = 0;
62			double eCaNegTime = 0;
63			double eZdNeg = 0;
64			double eZdPos = 0;
65			double eZdPosTime = 0;
66			double eZdNegTime = 0;
67			int ePxbHits = 0;
68	edwenger	1.5	int ePxHits = 0;
69	edwenger	1.2	double eClusVtxQual = 0;
70	edwenger	1.5	double eClusVtxDiff = 0;
71	loizides	1.1
72			Handle<HFRecHitCollection> hfhits;
73			try {
74			evt.getByLabel(edm::InputTag(srcHF_),hfhits);
75			} catch (...) {}
76			if (hfhits.isValid()) {
77			for(size_t ihit = 0; ihit<hfhits->size(); ++ihit){
78			const HFRecHit h = (*hfhits)[ihit];
79			double energy = h.energy();
80			double time = h.time();
81	edwenger	1.5	const HcalDetId id(h.id());
82			if(energy<0) continue;
83	loizides	1.1	if (id.zside()<0) {
84			eHfNeg += energy;
85			eHfNegTime += energy * time;
86			} else {
87			eHfPos += energy;
88			eHfPosTime += energy * time;
89			}
90			}
91			}
92
93			Handle<HBHERecHitCollection> hbhehits;
94			try {
95			evt.getByLabel(edm::InputTag(srcHBHE_),hbhehits);
96			} catch (...) {}
97			if (hbhehits.isValid()) {
98			for(size_t ihit = 0; ihit<hbhehits->size(); ++ihit){
99			const HBHERecHit h = (*hbhehits)[ihit];
100			double energy = h.energy();
101			const HcalDetId id(h.id());
102	edwenger	1.5	if(energy<0) continue;
103	loizides	1.1	if (id.zside()<0) {
104			eHcalNeg += energy;
105			} else {
106			eHcalPos += energy;
107			}
108			}
109			}
110
111			Handle<CastorRecHitCollection> castorhits;
112			try {
113			evt.getByLabel(edm::InputTag(srcCastor_),castorhits);
114			} catch (...) {}
115			if (castorhits.isValid()) {
116			for(size_t ihit = 0; ihit<castorhits->size(); ++ihit){
117			const CastorRecHit h = (*castorhits)[ihit];
118			double energy = h.energy();
119			double time = h.time();
120			const HcalCastorDetId id(h.id());
121	edwenger	1.5	if(energy<0) continue;
122	loizides	1.1	if (id.zside()<0) {
123			eCaNeg += energy;
124			eCaNegTime += energy * time;
125			} else {
126			eCaPos += energy;
127			eCaPosTime += energy * time;
128			}
129			}
130			}
131
132			Handle<ZDCRecHitCollection> zdchits;
133			try {
134			evt.getByLabel(edm::InputTag(srcZDC_),zdchits);
135			} catch (...) {}
136	edwenger	1.2	if (zdchits.isValid()) {
137	loizides	1.1	for(size_t ihit = 0; ihit<zdchits->size(); ++ihit){
138			const ZDCRecHit h = (*zdchits)[ihit];
139			double energy = h.energy();
140			double time = h.time();
141			const HcalZDCDetId id(h.id());
142	edwenger	1.5	if(energy<0) continue;
143	loizides	1.1	if (id.zside()<0) {
144			eZdNeg += energy;
145			eZdNegTime += energy * time;
146			} else {
147			eZdPos += energy;
148			eZdPosTime += energy * time;
149			}
150			}
151			}
152
153			if(eHfPos>0)
154			eHfPosTime /= eHfPos;
155			if(eHfNeg>0)
156			eHfNegTime /= eHfNeg;
157			if(eCaPos>0)
158			eCaPosTime /= eCaPos;
159			if(eCaNeg>0)
160			eCaNegTime /= eCaNeg;
161			if(eZdPos>0)
162			eZdPosTime /= eZdPos;
163			if(eZdNeg>0)
164			eZdNegTime /= eZdNeg;
165
166	loizides	1.4	Handle<SiPixelRecHitCollection> hRecHits;
167	edwenger	1.2	try {
168	loizides	1.4	evt.getByLabel(edm::InputTag(srcPixels_),hRecHits);
169	edwenger	1.2	} catch (...) {}
170
171	loizides	1.4	if (hRecHits.isValid()) {
172			ESHandle<TrackerGeometry> trackerHandle;
173			setup.get<TrackerDigiGeometryRecord>().get(trackerHandle);
174			const TrackerGeometry *tgeo = trackerHandle.product();
175			const SiPixelRecHitCollection *hits = hRecHits.product();
176	edwenger	1.5
177	loizides	1.8	vector<VertexHit> vhits(hits->dataSize());
178	loizides	1.4	for(SiPixelRecHitCollection::DataContainer::const_iterator hit = hits->data().begin(),
179			end = hits->data().end(); hit != end; ++hit) {
180			if (!hit->isValid())
181			continue;
182	loizides	1.8	++ePxHits;
183	loizides	1.4	DetId id(hit->geographicalId());
184			if(id.subdetId() != int(PixelSubdetector::PixelBarrel))
185			continue;
186	loizides	1.6	++ePxbHits;
187	loizides	1.4	const PixelGeomDetUnit pgdu = static_cast<const PixelGeomDetUnit>(tgeo->idToDet(id));
188			if (1) {
189			const RectangularPixelTopology *pixTopo =
190			static_cast<const RectangularPixelTopology*>(&(pgdu->specificTopology()));
191			vector<SiPixelCluster::Pixel> pixels(hit->cluster()->pixels());
192			bool pixelOnEdge = false;
193			for(std::vector<SiPixelCluster::Pixel>::const_iterator pixel = pixels.begin();
194			pixel != pixels.end(); ++pixel) {
195			int pixelX = pixel->x;
196			int pixelY = pixel->y;
197			if(pixTopo->isItEdgePixelInX(pixelX) \|\| pixTopo->isItEdgePixelInY(pixelY)) {
198			pixelOnEdge = true;
199			break;
200			}
201			}
202			if (pixelOnEdge)
203			continue;
204			}
205
206			LocalPoint lpos = LocalPoint(hit->localPosition().x(),
207			hit->localPosition().y(),
208			hit->localPosition().z());
209			GlobalPoint gpos = pgdu->toGlobal(lpos);
210			VertexHit vh;
211			vh.z = gpos.z();
212			vh.r = gpos.perp();
213			vh.w = hit->cluster()->sizeY();
214			vhits.push_back(vh);
215			}
216	edwenger	1.2
217	loizides	1.4	double zest = 0.0;
218			if (srcVertex_.length()!=0) {
219			edm::Handle<reco::VertexCollection> vertexCol;
220			evt.getByLabel(srcVertex_,vertexCol);
221			const reco::VertexCollection *vertices = vertexCol.product();
222			unsigned int vtracks = 0;
223			for(reco::VertexCollection::const_iterator
224			vertex = vertices->begin(); vertex!= vertices->end(); ++vertex) {
225			if(vertex->tracksSize()>vtracks) {
226			vtracks = vertex->tracksSize();
227			zest = vertex->z();
228			}
229			}
230			} else {
231			int nhits = 0, nhits_max = 0;
232			double chi = 0, chi_max = 1e+9;
233			for(double z0 = -15.9; z0 <= 15.95; z0 += 0.1) {
234			nhits = getContainedHits(vhits, z0, chi);
235			if(nhits > 0) {
236			if(nhits > nhits_max) {
237			chi_max = 1e+9;
238			nhits_max = nhits;
239			}
240			if(nhits >= nhits_max) {
241			if(chi < chi_max) {
242			chi_max = chi;
243			zest = z0;
244			}
245			}
246			}
247			}
248			}
249	edwenger	1.2
250	loizides	1.4	double chi = 0;
251			int nbest=0, nminus=0, nplus=0;
252			nbest = getContainedHits(vhits,zest,chi);
253			nminus = getContainedHits(vhits,zest-10.,chi);
254			nplus = getContainedHits(vhits,zest+10.,chi);
255	edwenger	1.5
256			eClusVtxDiff = nbest - (nminus+nplus)/2.;
257	loizides	1.7	if ((nminus+nplus)> 0)
258	edwenger	1.9	eClusVtxQual = (2.0*nbest)/(nminus+nplus); // A/B
259			else if (nbest>0)
260			eClusVtxQual = 1000.0; // A/0 (set to arbitrarily large number)
261			else
262	edwenger	1.10	eClusVtxQual = 0; // 0/0 (already the default)
263	edwenger	1.9
264	edwenger	1.2	}
265	loizides	1.4
266			std::auto_ptr<EvtSelData> output(new EvtSelData(eHcalNeg, eHcalPos,
267			eHfNeg,eHfPos,eHfNegTime,eHfPosTime,
268	loizides	1.1	eCaNeg,eCaPos,eCaNegTime,eCaPosTime,
269	edwenger	1.2	eZdNeg,eZdPos,eZdNegTime,eZdPosTime,
270	edwenger	1.5	ePxbHits,ePxHits,eClusVtxQual,eClusVtxDiff));
271	loizides	1.1	evt.put(output);
272			}
273
274	edwenger	1.2	//--------------------------------------------------------------------------------------------------
275	loizides	1.4	int ProducerEvtSelData::getContainedHits(const std::vector<VertexHit> &hits, double z0, double &chi)
276			{
277			// Calculate number of hits contained in v-shaped window in cluster y-width vs. z-position.
278	edwenger	1.2
279			int n = 0;
280	loizides	1.4	chi = 0.;
281	edwenger	1.2
282	loizides	1.4	for(vector<VertexHit>::const_iterator hit = hits.begin(); hit!= hits.end(); hit++) {
283			double p = 2 * fabs(hit->z - z0)/hit->r + 0.5; // FIXME
284			if(TMath::Abs(p - hit->w) <= 1.) {
285	edwenger	1.2	chi += fabs(p - hit->w);
286			n++;
287			}
288			}
289			return n;
290			}
291
292	loizides	1.1	//define this as a plug-in
293			DEFINE_FWK_MODULE(ProducerEvtSelData);