194 typedef recHitContainer::const_iterator const_iterator;
195 typedef std::pair<const_iterator, const_iterator>
range;
205 edm::LogInfo(
"") <<
"\n\n =================> Treating event " << e.
id() <<
" Number of seeds " 209 float mass = .000511;
213 float dphi1 = 0., dphi2 = 0., drz1 = 0., drz2 = 0.;
214 float phi1 = 0., phi2 = 0., rz1 = 0., rz2 = 0.;
216 for (ElectronSeedCollection::const_iterator MyS = (*elSeeds).begin(); MyS != (*elSeeds).end(); ++MyS) {
217 LogDebug(
"") <<
"\nSeed nr " << is <<
": ";
218 range
r = (*MyS).recHits();
219 LogDebug(
"") <<
" Number of RecHits= " << (*MyS).nHits();
226 LogDebug(
"") <<
" First hit local x,y,z " << (*it).localPosition() <<
" det " << id1.
det() <<
" subdet " 228 LogDebug(
"") <<
" First hit global " << det1->
toGlobal((*it).localPosition());
234 LogDebug(
"") <<
" Second hit local x,y,z " << (*it).localPosition() <<
" det " << id2.
det() <<
" subdet " 236 LogDebug(
"") <<
" Second hit global " << det2->
toGlobal((*it).localPosition());
242 for (TrackingRecHitCollection::const_iterator rhits = r.first; rhits != r.second; rhits++)
243 det = pDD->
idToDet(((*rhits)).geographicalId());
255 LogDebug(
"") <<
" ElectronSeed superCluster energy: " << theClus->energy() <<
", position: " << theClus->position();
257 LogDebug(
"") <<
" ElectronSeed supercluster Et: " 258 << theClus->energy() *
sin(2. * atan(
exp(-theClus->position().eta())));
260 LogDebug(
"") <<
" ElectronSeed supercluster eta: " << theClus->position().eta();
261 LogDebug(
"") <<
" ElectronSeed seed charge: " << (*MyS).getCharge();
269 GlobalPoint xmeas(theClus->position().x(), theClus->position().y(), theClus->position().z());
271 float energy = theClus->energy();
285 DetId id = (*it).geographicalId();
295 float SCl_phi = xmeas.phi();
296 float localDphi = SCl_phi - hitPos.
phi();
308 if (
id.subdetId() % 2 == 1) {
315 DetId id2 = (*it).geographicalId();
320 double pxHit1z = hitPos.
z();
321 double pxHit1x = hitPos.
x();
322 double pxHit1y = hitPos.
y();
323 double r1diff = (pxHit1x - vprim.x()) * (pxHit1x - vprim.x()) + (pxHit1y - vprim.y()) * (pxHit1y - vprim.y());
324 r1diff =
sqrt(r1diff);
325 double r2diff = (xmeas.x() - pxHit1x) * (xmeas.x() - pxHit1x) + (xmeas.y() - pxHit1y) * (xmeas.y() - pxHit1y);
326 r2diff =
sqrt(r2diff);
327 double zVertexPred = pxHit1z - r1diff * (xmeas.z() - pxHit1z) / r2diff;
329 GlobalPoint vertexPred(vprim.x(), vprim.y(), zVertexPred);
337 phi2 = hitPos2.
phi();
339 rz2 = hitPos2.
perp();
351 histetclu_->Fill(theClus->energy() *
sin(2. * atan(
exp(-theClus->position().eta()))));
354 histq_->Fill((*MyS).getCharge());
366 seedQ[is] = (*MyS).getCharge();
392 e.
getByLabel(
"correctedHybridSuperClusters", clusters);
394 if (!clusters.
product()->empty())
402 e.
getByLabel(
"generatorSmeared",
"", HepMCEvt);
406 HepMC::FourVector pAssSim;
408 for (HepMC::GenEvent::particle_const_iterator partIter = MCEvt->particles_begin(); partIter != MCEvt->particles_end();
410 for (HepMC::GenEvent::vertex_const_iterator vertIter = MCEvt->vertices_begin(); vertIter != MCEvt->vertices_end();
413 HepMC::GenVertex *creation = (*partIter)->production_vertex();
415 HepMC::FourVector momentum = (*partIter)->momentum();
417 int id = (*partIter)->pdg_id();
418 LogDebug(
"") <<
"MC particle id " <<
id <<
", creationVertex " << (*creation) <<
" cm, initialMomentum " 419 << momentum.rho() <<
" GeV/c" << std::endl;
421 if (
id == 11 ||
id == -11) {
424 if ((*partIter)->production_vertex()) {
425 if ((*partIter)->production_vertex()->particles_begin(
HepMC::parents) !=
427 mother = *((*partIter)->production_vertex()->particles_begin(
HepMC::parents));
429 if (((mother ==
nullptr) || ((mother !=
nullptr) && (mother->pdg_id() == 23)) ||
430 ((mother !=
nullptr) && (mother->pdg_id() == 32)) ||
431 ((mother !=
nullptr) && (
std::abs(mother->pdg_id()) == 24)))) {
433 pAssSim = genPc->momentum();
448 bool okSeedFound =
false;
449 double seedOkRatio = 999999.;
453 for (ElectronSeedCollection::const_iterator gsfIter = (*elSeeds).begin(); gsfIter != (*elSeeds).end();
457 for (TrackingRecHitCollection::const_iterator rhits = r.first; rhits != r.second; rhits++)
458 det = pDD->
idToDet(((*rhits)).geographicalId());
465 double dphi = phi - pAssSim.phi();
467 dphi = dphi < 0 ? (CLHEP::twopi) + dphi : dphi - CLHEP::twopi;
473 double tmpSeedRatio = p / pAssSim.t();
475 seedOkRatio = tmpSeedRatio;
476 bestElectronSeed = *gsfIter;
490 mcEta[ip] = pAssSim.eta();
491 mcPhi[ip] = pAssSim.phi();
492 mcPt[ip] = pAssSim.perp();
493 mcQ[ip] = ((
id == 11) ? -1. : +1.);
499 seedOkRatio = 999999.;
502 for (ElectronSeedCollection::const_iterator gsfIter = (*elSeeds).begin(); gsfIter != (*elSeeds).end();
504 range r = gsfIter->recHits();
506 for (TrackingRecHitCollection::const_iterator rhits = r.first; rhits != r.second; rhits++)
507 det = pDD->
idToDet(((*rhits)).geographicalId());
514 double dphi = phi - pAssSim.phi();
516 dphi = dphi < 0 ? (CLHEP::twopi) + dphi : dphi - CLHEP::twopi;
518 if (gsfIter->isEcalDriven()) {
523 double tmpSeedRatio = p / pAssSim.t();
525 seedOkRatio = tmpSeedRatio;
526 bestElectronSeed = *gsfIter;
543 seedOkRatio = 999999.;
546 for (ElectronSeedCollection::const_iterator gsfIter = (*elSeeds).begin(); gsfIter != (*elSeeds).end();
548 range r = gsfIter->recHits();
550 for (TrackingRecHitCollection::const_iterator rhits = r.first; rhits != r.second; rhits++)
551 det = pDD->
idToDet(((*rhits)).geographicalId());
558 double dphi = phi - pAssSim.phi();
560 dphi = dphi < 0 ? (CLHEP::twopi) + dphi : dphi - CLHEP::twopi;
562 if (gsfIter->isTrackerDriven()) {
567 double tmpSeedRatio = p / pAssSim.t();
569 seedOkRatio = tmpSeedRatio;
570 bestElectronSeed = *gsfIter;
GlobalPoint toGlobal(const Point2DBase< Scalar, LocalTag > lp) const
static FreeTrajectoryState get(MagneticField const &magField, GlobalPoint const &xmeas, GlobalPoint const &xvert, float momentum, TrackCharge charge)
float superclusterEnergy[10]
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Sin< T >::type sin(const T &t)
float superclusterPhi[10]
Geom::Phi< T > phi() const
unsigned int side(const DetId &id) const
GlobalPoint globalPosition() const
TH1F * histtrackerdrivenetaMCmatched_
const Plane & surface() const
The nominal surface of the GeomDet.
edm::InputTag inputCollection_
recHitContainer::const_iterator const_iterator
TH1F * histecaldrivenptMCmatched_
constexpr int subdetId() const
get the contents of the subdetector field (not cast into any detector's numbering enum) ...
Abs< T >::type abs(const T &t)
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
GlobalVector momentum() const
TH1F * histecaldriveneMCmatched_
GlobalPoint position() const
float superclusterEta[10]
const HepMC::GenEvent * GetEvent() const
TH1F * histecaldrivenetaMCmatched_
T const * product() const
bool isNull() const
Checks for null.
TrajectoryStateOnSurface propagate(STA const &state, SUR const &surface) const
TH1F * histtrackerdriveneMCmatched_
unsigned int layer(const DetId &id) const
GlobalVector globalMomentum() const
const TrackerGeomDet * idToDet(DetId) const override
const Point & position() const
position
Power< A, B >::type pow(const A &a, const B &b)
TH1F * histtrackerdrivenptMCmatched_
constexpr Detector det() const
get the detector field from this detid