Producers/src/ProducerEvtSelData.cc

// $Id: ProducerEvtSelData.cc,v 1.12 2010/01/07 17:07:55 loizides Exp $

#include "MitEdm/Producers/interface/ProducerEvtSelData.h"
#include "MitEdm/DataFormats/interface/EvtSelData.h"
#include "DataFormats/CaloTowers/interface/CaloTower.h"
#include "DataFormats/CaloTowers/interface/CaloTowerFwd.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 "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/TrackReco/interface/Track.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","")),
    srcTowers_(cfg.getUntrackedParameter<string>("caloTowers","towerMaker")),
    hfEthresh_(cfg.getUntrackedParameter<double>("hfEThreshold",3.)),
    hfETowerh_(cfg.getUntrackedParameter<double>("hfETowerThreshold",3.)),
    srcTrk_(cfg.getUntrackedParameter<string>("tracks","generalTracks"))
{
  // 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;
  int    eHfPcounts   = 0;
  int    eHfNcounts   = 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;
  int    nHfTowersP      = 0;
  int    nHfTowersN      = 0;
  double sumEsubEpPos    = 0;
  double sumEaddEpPos    = 0;
  double sumHfEsubEpPos  = 0;
  double sumHfEaddEpPos  = 0;
  double sumEsubEpNeg    = 0;
  double sumEaddEpNeg    = 0;
  double sumHfEsubEpNeg  = 0;
  double sumHfEaddEpNeg  = 0;
  double eHPTrkFrac      = 0;

  // HF
  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;
        if (energy>hfEthresh_)
          ++eHfNcounts;
      } else {
        eHfPos     += energy;
        eHfPosTime += energy * time;
        if (energy>hfEthresh_)
          ++eHfPcounts;
      }
    }
  }

  // HCAL (barrel + endcap)
  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;
      }
    }
  }

  // CASTOR
  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;
      }
    }
  }

  // ZDC
  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;
      }
    }
  }

  // time values are energy-weighted averages
  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;
  

  // pixel cluster shape compatibility with vertex
  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;
    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) { // use specified vertex
      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 { // or use pixel cluster-shape vertex
      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)

  }

  // calo based variables
  edm::Handle<CaloTowerCollection> towers; 
  try {
    evt.getByLabel(edm::InputTag(srcTowers_),towers);
  } catch (...) {}  

  if (towers.isValid()) {
    for(CaloTowerCollection::const_iterator cal = towers->begin(); cal != towers->end(); ++cal) {
      if (cal->energy()<hfETowerh_)
        continue;

      double tmp = TMath::Cos(cal->theta());
      double Esub = cal->energy()*(1-tmp);
      double Eadd = cal->energy()*(1+tmp);

      if (cal->eta()>0) {
        sumEsubEpPos += Esub;
        sumEaddEpPos += Eadd;
      } else if (cal->eta()<0) {
        sumEsubEpNeg += Esub;
        sumEaddEpNeg += Eadd;
      }
      if (cal->eta()>3) {
        sumHfEsubEpPos += Esub;
        sumHfEaddEpPos += Eadd;
      } else if (cal->eta()<-3) {
        sumHfEsubEpNeg += Esub;
        sumHfEaddEpNeg += Eadd;
      }

      for(unsigned int i = 0; i < cal->constituentsSize(); ++i) {
        const DetId id = cal->constituent(i);
        if(id.det() != DetId::Hcal)
          continue;
        HcalSubdetector subdet=(HcalSubdetector(id.subdetId()));
        if(subdet != HcalForward)
          continue;
        if (cal->eta()<-3)
          ++nHfTowersN;
        if (cal->eta()>+3)
          ++nHfTowersP;
      }    
    }
  }

  // high-purity track fraction
  Handle<reco::TrackCollection> tkRef;
  try {
    evt.getByLabel(edm::InputTag(srcTrk_),tkRef);
  } catch (...) {}  

  if (tkRef.isValid()) {

    const reco::TrackCollection* tkColl = tkRef.product();
    int numhighpurity=0;
    reco::TrackBase::TrackQuality trackQuality(reco::TrackBase::qualityByName("highPurity"));
    if(tkColl->size()>0){ 
      reco::TrackCollection::const_iterator itk = tkColl->begin();
      reco::TrackCollection::const_iterator itk_e = tkColl->end();
      for(;itk!=itk_e;++itk){
        if(itk->quality(trackQuality)) numhighpurity++;
      }
      eHPTrkFrac = (double)numhighpurity/(double)tkColl->size();
    }
  }

  // fill EvtSelData object
  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,
                                                  eHfPcounts, eHfNcounts,
                                                  nHfTowersP, nHfTowersN,
                                                  sumEsubEpPos,  sumEaddEpPos,
                                                  sumHfEsubEpPos,sumHfEaddEpPos,
                                                  sumEsubEpNeg,  sumEaddEpNeg,
                                                  sumHfEsubEpNeg,sumHfEaddEpNeg,
                                                  eHPTrkFrac));
  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.13
Committed:	Tue Aug 24 22:14:48 2010 UTC (14 years, 8 months ago) by edwenger
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
Changes since 1.12:	+2 -2 lines
Log Message:	much faster without initializing vector
#	User	Rev	Content
1	edwenger	1.13	// $Id: ProducerEvtSelData.cc,v 1.12 2010/01/07 17:07:55 loizides Exp $
2	loizides	1.1
3			#include "MitEdm/Producers/interface/ProducerEvtSelData.h"
4			#include "MitEdm/DataFormats/interface/EvtSelData.h"
5	loizides	1.12	#include "DataFormats/CaloTowers/interface/CaloTower.h"
6			#include "DataFormats/CaloTowers/interface/CaloTowerFwd.h"
7	loizides	1.1	#include "DataFormats/Common/interface/Handle.h"
8	loizides	1.4	#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
9			#include "DataFormats/GeometryVector/interface/LocalPoint.h"
10			#include "DataFormats/HcalDetId/interface/HcalDetId.h"
11	loizides	1.1	#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
12	loizides	1.4	#include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
13			#include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
14			#include "DataFormats/VertexReco/interface/Vertex.h"
15			#include "DataFormats/VertexReco/interface/VertexFwd.h"
16	edwenger	1.11	#include "DataFormats/TrackReco/interface/TrackFwd.h"
17			#include "DataFormats/TrackReco/interface/Track.h"
18	loizides	1.1	#include "FWCore/Framework/interface/ESHandle.h"
19			#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
20	loizides	1.4	#include "Geometry/CommonDetUnit/interface/GeomDet.h"
21	loizides	1.1	#include "Geometry/Records/interface/CaloGeometryRecord.h"
22	loizides	1.4	#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
23			#include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetUnit.h"
24	edwenger	1.2	#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
25			#include "Geometry/TrackerTopology/interface/RectangularPixelTopology.h"
26	loizides	1.4	#include <TMath.h>
27	edwenger	1.2
28	loizides	1.1	using namespace std;
29			using namespace edm;
30			using namespace mitedm;
31
32			//--------------------------------------------------------------------------------------------------
33			ProducerEvtSelData::ProducerEvtSelData(const ParameterSet& cfg)
34			: srcHF_(cfg.getUntrackedParameter<string>("hfRecHits","hfreco")),
35			srcHBHE_(cfg.getUntrackedParameter<string>("hbheRecHits","hbhereco")),
36			srcCastor_(cfg.getUntrackedParameter<string>("castorRecHits","castorreco")),
37	edwenger	1.2	srcZDC_(cfg.getUntrackedParameter<string>("zdcRecHits","zdcreco")),
38	loizides	1.4	srcPixels_(cfg.getUntrackedParameter<string>("pixelRecHits","siPixelRecHits")),
39	edwenger	1.11	srcVertex_(cfg.getUntrackedParameter<string>("vertex","")),
40	loizides	1.12	srcTowers_(cfg.getUntrackedParameter<string>("caloTowers","towerMaker")),
41			hfEthresh_(cfg.getUntrackedParameter<double>("hfEThreshold",3.)),
42			hfETowerh_(cfg.getUntrackedParameter<double>("hfETowerThreshold",3.)),
43	edwenger	1.11	srcTrk_(cfg.getUntrackedParameter<string>("tracks","generalTracks"))
44	loizides	1.1	{
45			// Constructor.
46
47			produces<EvtSelData>();
48			}
49
50			//--------------------------------------------------------------------------------------------------
51			ProducerEvtSelData::~ProducerEvtSelData()
52			{
53			// Destructor.
54			}
55
56			//--------------------------------------------------------------------------------------------------
57			void ProducerEvtSelData::produce(Event &evt, const EventSetup &setup)
58			{
59			// Produce event selection data for this event.
60
61	loizides	1.4	double eHfNeg = 0;
62			double eHfPos = 0;
63			double eHfPosTime = 0;
64			double eHfNegTime = 0;
65	loizides	1.12	int eHfPcounts = 0;
66			int eHfNcounts = 0;
67	loizides	1.4	double eHcalNeg = 0;
68			double eHcalPos = 0;
69			double eCaNeg = 0;
70			double eCaPos = 0;
71			double eCaPosTime = 0;
72			double eCaNegTime = 0;
73			double eZdNeg = 0;
74			double eZdPos = 0;
75			double eZdPosTime = 0;
76			double eZdNegTime = 0;
77			int ePxbHits = 0;
78	edwenger	1.5	int ePxHits = 0;
79	edwenger	1.2	double eClusVtxQual = 0;
80	edwenger	1.5	double eClusVtxDiff = 0;
81	loizides	1.12	int nHfTowersP = 0;
82			int nHfTowersN = 0;
83			double sumEsubEpPos = 0;
84			double sumEaddEpPos = 0;
85			double sumHfEsubEpPos = 0;
86			double sumHfEaddEpPos = 0;
87			double sumEsubEpNeg = 0;
88			double sumEaddEpNeg = 0;
89			double sumHfEsubEpNeg = 0;
90			double sumHfEaddEpNeg = 0;
91			double eHPTrkFrac = 0;
92	edwenger	1.11
93			// HF
94	loizides	1.1	Handle<HFRecHitCollection> hfhits;
95			try {
96			evt.getByLabel(edm::InputTag(srcHF_),hfhits);
97			} catch (...) {}
98	loizides	1.12
99	loizides	1.1	if (hfhits.isValid()) {
100			for(size_t ihit = 0; ihit<hfhits->size(); ++ihit){
101			const HFRecHit h = (*hfhits)[ihit];
102			double energy = h.energy();
103			double time = h.time();
104	edwenger	1.5	const HcalDetId id(h.id());
105			if(energy<0) continue;
106	loizides	1.1	if (id.zside()<0) {
107			eHfNeg += energy;
108			eHfNegTime += energy * time;
109	loizides	1.12	if (energy>hfEthresh_)
110			++eHfNcounts;
111	loizides	1.1	} else {
112			eHfPos += energy;
113			eHfPosTime += energy * time;
114	loizides	1.12	if (energy>hfEthresh_)
115			++eHfPcounts;
116	loizides	1.1	}
117			}
118			}
119
120	edwenger	1.11	// HCAL (barrel + endcap)
121	loizides	1.1	Handle<HBHERecHitCollection> hbhehits;
122			try {
123			evt.getByLabel(edm::InputTag(srcHBHE_),hbhehits);
124			} catch (...) {}
125	loizides	1.12
126	loizides	1.1	if (hbhehits.isValid()) {
127			for(size_t ihit = 0; ihit<hbhehits->size(); ++ihit){
128			const HBHERecHit h = (*hbhehits)[ihit];
129			double energy = h.energy();
130			const HcalDetId id(h.id());
131	edwenger	1.5	if(energy<0) continue;
132	loizides	1.1	if (id.zside()<0) {
133			eHcalNeg += energy;
134			} else {
135			eHcalPos += energy;
136			}
137			}
138			}
139
140	edwenger	1.11	// CASTOR
141	loizides	1.1	Handle<CastorRecHitCollection> castorhits;
142			try {
143			evt.getByLabel(edm::InputTag(srcCastor_),castorhits);
144			} catch (...) {}
145	loizides	1.12
146	loizides	1.1	if (castorhits.isValid()) {
147			for(size_t ihit = 0; ihit<castorhits->size(); ++ihit){
148			const CastorRecHit h = (*castorhits)[ihit];
149			double energy = h.energy();
150			double time = h.time();
151			const HcalCastorDetId id(h.id());
152	edwenger	1.5	if(energy<0) continue;
153	loizides	1.1	if (id.zside()<0) {
154			eCaNeg += energy;
155			eCaNegTime += energy * time;
156			} else {
157			eCaPos += energy;
158			eCaPosTime += energy * time;
159			}
160			}
161			}
162
163	edwenger	1.11	// ZDC
164	loizides	1.1	Handle<ZDCRecHitCollection> zdchits;
165			try {
166			evt.getByLabel(edm::InputTag(srcZDC_),zdchits);
167			} catch (...) {}
168	loizides	1.12
169	edwenger	1.2	if (zdchits.isValid()) {
170	loizides	1.1	for(size_t ihit = 0; ihit<zdchits->size(); ++ihit){
171			const ZDCRecHit h = (*zdchits)[ihit];
172			double energy = h.energy();
173			double time = h.time();
174			const HcalZDCDetId id(h.id());
175	edwenger	1.5	if(energy<0) continue;
176	loizides	1.1	if (id.zside()<0) {
177			eZdNeg += energy;
178			eZdNegTime += energy * time;
179			} else {
180			eZdPos += energy;
181			eZdPosTime += energy * time;
182			}
183			}
184			}
185
186	edwenger	1.11	// time values are energy-weighted averages
187	loizides	1.1	if(eHfPos>0)
188			eHfPosTime /= eHfPos;
189			if(eHfNeg>0)
190			eHfNegTime /= eHfNeg;
191			if(eCaPos>0)
192			eCaPosTime /= eCaPos;
193			if(eCaNeg>0)
194			eCaNegTime /= eCaNeg;
195			if(eZdPos>0)
196			eZdPosTime /= eZdPos;
197			if(eZdNeg>0)
198			eZdNegTime /= eZdNeg;
199
200	edwenger	1.11
201			// pixel cluster shape compatibility with vertex
202	loizides	1.4	Handle<SiPixelRecHitCollection> hRecHits;
203	edwenger	1.2	try {
204	loizides	1.4	evt.getByLabel(edm::InputTag(srcPixels_),hRecHits);
205	edwenger	1.2	} catch (...) {}
206
207	loizides	1.4	if (hRecHits.isValid()) {
208			ESHandle<TrackerGeometry> trackerHandle;
209			setup.get<TrackerDigiGeometryRecord>().get(trackerHandle);
210			const TrackerGeometry *tgeo = trackerHandle.product();
211			const SiPixelRecHitCollection *hits = hRecHits.product();
212	edwenger	1.5
213	edwenger	1.13	vector<VertexHit> vhits;
214	loizides	1.4	for(SiPixelRecHitCollection::DataContainer::const_iterator hit = hits->data().begin(),
215			end = hits->data().end(); hit != end; ++hit) {
216			if (!hit->isValid())
217			continue;
218	loizides	1.8	++ePxHits;
219	loizides	1.4	DetId id(hit->geographicalId());
220			if(id.subdetId() != int(PixelSubdetector::PixelBarrel))
221			continue;
222	loizides	1.6	++ePxbHits;
223	loizides	1.4	const PixelGeomDetUnit pgdu = static_cast<const PixelGeomDetUnit>(tgeo->idToDet(id));
224			if (1) {
225			const RectangularPixelTopology *pixTopo =
226			static_cast<const RectangularPixelTopology*>(&(pgdu->specificTopology()));
227			vector<SiPixelCluster::Pixel> pixels(hit->cluster()->pixels());
228			bool pixelOnEdge = false;
229			for(std::vector<SiPixelCluster::Pixel>::const_iterator pixel = pixels.begin();
230			pixel != pixels.end(); ++pixel) {
231			int pixelX = pixel->x;
232			int pixelY = pixel->y;
233			if(pixTopo->isItEdgePixelInX(pixelX) \|\| pixTopo->isItEdgePixelInY(pixelY)) {
234			pixelOnEdge = true;
235			break;
236			}
237			}
238			if (pixelOnEdge)
239			continue;
240			}
241
242			LocalPoint lpos = LocalPoint(hit->localPosition().x(),
243			hit->localPosition().y(),
244			hit->localPosition().z());
245			GlobalPoint gpos = pgdu->toGlobal(lpos);
246			VertexHit vh;
247			vh.z = gpos.z();
248			vh.r = gpos.perp();
249			vh.w = hit->cluster()->sizeY();
250			vhits.push_back(vh);
251			}
252	edwenger	1.2
253	loizides	1.4	double zest = 0.0;
254	edwenger	1.11	if (srcVertex_.length()!=0) { // use specified vertex
255	loizides	1.4	edm::Handle<reco::VertexCollection> vertexCol;
256			evt.getByLabel(srcVertex_,vertexCol);
257			const reco::VertexCollection *vertices = vertexCol.product();
258			unsigned int vtracks = 0;
259			for(reco::VertexCollection::const_iterator
260			vertex = vertices->begin(); vertex!= vertices->end(); ++vertex) {
261			if(vertex->tracksSize()>vtracks) {
262			vtracks = vertex->tracksSize();
263			zest = vertex->z();
264			}
265			}
266	edwenger	1.11	} else { // or use pixel cluster-shape vertex
267	loizides	1.4	int nhits = 0, nhits_max = 0;
268			double chi = 0, chi_max = 1e+9;
269			for(double z0 = -15.9; z0 <= 15.95; z0 += 0.1) {
270			nhits = getContainedHits(vhits, z0, chi);
271			if(nhits > 0) {
272			if(nhits > nhits_max) {
273			chi_max = 1e+9;
274			nhits_max = nhits;
275			}
276			if(nhits >= nhits_max) {
277			if(chi < chi_max) {
278			chi_max = chi;
279			zest = z0;
280			}
281			}
282			}
283			}
284			}
285	edwenger	1.2
286	loizides	1.4	double chi = 0;
287			int nbest=0, nminus=0, nplus=0;
288			nbest = getContainedHits(vhits,zest,chi);
289			nminus = getContainedHits(vhits,zest-10.,chi);
290			nplus = getContainedHits(vhits,zest+10.,chi);
291	edwenger	1.5
292			eClusVtxDiff = nbest - (nminus+nplus)/2.;
293	loizides	1.7	if ((nminus+nplus)> 0)
294	edwenger	1.9	eClusVtxQual = (2.0*nbest)/(nminus+nplus); // A/B
295			else if (nbest>0)
296			eClusVtxQual = 1000.0; // A/0 (set to arbitrarily large number)
297			else
298	edwenger	1.10	eClusVtxQual = 0; // 0/0 (already the default)
299	edwenger	1.9
300	edwenger	1.2	}
301	loizides	1.4
302	loizides	1.12	// calo based variables
303			edm::Handle<CaloTowerCollection> towers;
304			try {
305			evt.getByLabel(edm::InputTag(srcTowers_),towers);
306			} catch (...) {}
307
308			if (towers.isValid()) {
309			for(CaloTowerCollection::const_iterator cal = towers->begin(); cal != towers->end(); ++cal) {
310			if (cal->energy()<hfETowerh_)
311			continue;
312
313			double tmp = TMath::Cos(cal->theta());
314			double Esub = cal->energy()*(1-tmp);
315			double Eadd = cal->energy()*(1+tmp);
316
317			if (cal->eta()>0) {
318			sumEsubEpPos += Esub;
319			sumEaddEpPos += Eadd;
320			} else if (cal->eta()<0) {
321			sumEsubEpNeg += Esub;
322			sumEaddEpNeg += Eadd;
323			}
324			if (cal->eta()>3) {
325			sumHfEsubEpPos += Esub;
326			sumHfEaddEpPos += Eadd;
327			} else if (cal->eta()<-3) {
328			sumHfEsubEpNeg += Esub;
329			sumHfEaddEpNeg += Eadd;
330			}
331
332			for(unsigned int i = 0; i < cal->constituentsSize(); ++i) {
333			const DetId id = cal->constituent(i);
334			if(id.det() != DetId::Hcal)
335			continue;
336			HcalSubdetector subdet=(HcalSubdetector(id.subdetId()));
337			if(subdet != HcalForward)
338			continue;
339			if (cal->eta()<-3)
340			++nHfTowersN;
341			if (cal->eta()>+3)
342			++nHfTowersP;
343			}
344			}
345			}
346	edwenger	1.11
347			// high-purity track fraction
348			Handle<reco::TrackCollection> tkRef;
349			try {
350			evt.getByLabel(edm::InputTag(srcTrk_),tkRef);
351			} catch (...) {}
352
353	loizides	1.12	if (tkRef.isValid()) {
354
355			const reco::TrackCollection* tkColl = tkRef.product();
356			int numhighpurity=0;
357			reco::TrackBase::TrackQuality trackQuality(reco::TrackBase::qualityByName("highPurity"));
358			if(tkColl->size()>0){
359			reco::TrackCollection::const_iterator itk = tkColl->begin();
360			reco::TrackCollection::const_iterator itk_e = tkColl->end();
361			for(;itk!=itk_e;++itk){
362			if(itk->quality(trackQuality)) numhighpurity++;
363			}
364			eHPTrkFrac = (double)numhighpurity/(double)tkColl->size();
365	edwenger	1.11	}
366			}
367
368			// fill EvtSelData object
369	loizides	1.4	std::auto_ptr<EvtSelData> output(new EvtSelData(eHcalNeg, eHcalPos,
370			eHfNeg,eHfPos,eHfNegTime,eHfPosTime,
371	loizides	1.1	eCaNeg,eCaPos,eCaNegTime,eCaPosTime,
372	edwenger	1.2	eZdNeg,eZdPos,eZdNegTime,eZdPosTime,
373	loizides	1.12	ePxbHits,ePxHits,
374			eClusVtxQual,eClusVtxDiff,
375			eHfPcounts, eHfNcounts,
376			nHfTowersP, nHfTowersN,
377			sumEsubEpPos, sumEaddEpPos,
378			sumHfEsubEpPos,sumHfEaddEpPos,
379			sumEsubEpNeg, sumEaddEpNeg,
380			sumHfEsubEpNeg,sumHfEaddEpNeg,
381			eHPTrkFrac));
382	loizides	1.1	evt.put(output);
383			}
384
385	edwenger	1.2	//--------------------------------------------------------------------------------------------------
386	loizides	1.4	int ProducerEvtSelData::getContainedHits(const std::vector<VertexHit> &hits, double z0, double &chi)
387			{
388			// Calculate number of hits contained in v-shaped window in cluster y-width vs. z-position.
389	edwenger	1.2
390			int n = 0;
391	loizides	1.4	chi = 0.;
392	edwenger	1.2
393	loizides	1.4	for(vector<VertexHit>::const_iterator hit = hits.begin(); hit!= hits.end(); hit++) {
394			double p = 2 * fabs(hit->z - z0)/hit->r + 0.5; // FIXME
395			if(TMath::Abs(p - hit->w) <= 1.) {
396	edwenger	1.2	chi += fabs(p - hit->w);
397			n++;
398			}
399			}
400			return n;
401			}
402
403	loizides	1.1	//define this as a plug-in
404			DEFINE_FWK_MODULE(ProducerEvtSelData);