258 edm::LogInfo(
"") <<
"\n\n =================> Treating event " << e.
id() <<
" Number of seeds "
262 float mass = .000511;
266 float dphi1 = 0., dphi2 = 0., drz1 = 0., drz2 = 0.;
267 float phi1 = 0., phi2 = 0., rz1 = 0., rz2 = 0.;
269 for (ElectronSeedCollection::const_iterator MyS = (*elSeeds).begin(); MyS != (*elSeeds).end(); ++MyS) {
270 LogDebug(
"") <<
"\nSeed nr " << is <<
": ";
272 LogDebug(
"") <<
" Number of RecHits= " << (*MyS).nHits();
277 DetId id1 = (*it).geographicalId();
278 det1 = pDD->idToDet(id1);
279 LogDebug(
"") <<
" First hit local x,y,z " << (*it).localPosition() <<
" det " << id1.
det() <<
" subdet "
281 LogDebug(
"") <<
" First hit global " << det1->
toGlobal((*it).localPosition());
285 DetId id2 = (*it).geographicalId();
286 det2 = pDD->idToDet(id2);
287 LogDebug(
"") <<
" Second hit local x,y,z " << (*it).localPosition() <<
" det " << id2.
det() <<
" subdet "
289 LogDebug(
"") <<
" Second hit global " << det2->
toGlobal((*it).localPosition());
295 for (
auto const &recHit : r) {
296 det = pDD->idToDet(recHit.geographicalId());
309 LogDebug(
"") <<
" ElectronSeed superCluster energy: " << theClus->energy() <<
", position: " << theClus->position();
311 LogDebug(
"") <<
" ElectronSeed supercluster Et: "
312 << theClus->energy() *
sin(2. * atan(
exp(-theClus->position().eta())));
314 LogDebug(
"") <<
" ElectronSeed supercluster eta: " << theClus->position().eta();
315 LogDebug(
"") <<
" ElectronSeed seed charge: " << (*MyS).getCharge();
322 int charge = int((*MyS).getCharge());
323 GlobalPoint xmeas(theClus->position().x(), theClus->position().y(), theClus->position().z());
324 GlobalPoint vprim(theBeamSpot->position().x(), theBeamSpot->position().y(), theBeamSpot->position().z());
325 float energy = theClus->energy();
339 DetId id = (*it).geographicalId();
340 const GeomDet *geomdet = pDD->idToDet((*it).geographicalId());
349 float SCl_phi = xmeas.phi();
350 float localDphi = SCl_phi - hitPos.
phi();
362 if (
id.subdetId() % 2 == 1) {
369 DetId id2 = (*it).geographicalId();
370 const GeomDet *geomdet2 = pDD->idToDet((*it).geographicalId());
374 double pxHit1z = hitPos.
z();
375 double pxHit1x = hitPos.
x();
376 double pxHit1y = hitPos.
y();
377 double r1diff = (pxHit1x - vprim.x()) * (pxHit1x - vprim.x()) + (pxHit1y - vprim.y()) * (pxHit1y - vprim.y());
378 r1diff =
sqrt(r1diff);
379 double r2diff = (xmeas.x() - pxHit1x) * (xmeas.x() - pxHit1x) + (xmeas.y() - pxHit1y) * (xmeas.y() - pxHit1y);
380 r2diff =
sqrt(r2diff);
381 double zVertexPred = pxHit1z - r1diff * (xmeas.z() - pxHit1z) / r2diff;
383 GlobalPoint vertexPred(vprim.x(), vprim.y(), zVertexPred);
391 phi2 = hitPos2.
phi();
393 rz2 = hitPos2.
perp();
405 histetclu_->Fill(theClus->energy() *
sin(2. * atan(
exp(-theClus->position().eta()))));
408 histq_->Fill((*MyS).getCharge());
420 seedQ[is] = (*MyS).getCharge();
446 e.
getByLabel(
"correctedHybridSuperClusters", clusters);
448 if (!clusters.
product()->empty())
456 e.
getByLabel(
"generatorSmeared",
"", HepMCEvt);
460 HepMC::FourVector pAssSim;
462 for (HepMC::GenEvent::particle_const_iterator partIter = MCEvt->particles_begin(); partIter != MCEvt->particles_end();
464 for (HepMC::GenEvent::vertex_const_iterator vertIter = MCEvt->vertices_begin(); vertIter != MCEvt->vertices_end();
467 HepMC::GenVertex *creation = (*partIter)->production_vertex();
469 HepMC::FourVector momentum = (*partIter)->momentum();
471 int id = (*partIter)->pdg_id();
472 LogDebug(
"") <<
"MC particle id " <<
id <<
", creationVertex " << (*creation) <<
" cm, initialMomentum "
473 << momentum.rho() <<
" GeV/c" << std::endl;
475 if (
id == 11 ||
id == -11) {
478 if ((*partIter)->production_vertex()) {
479 if ((*partIter)->production_vertex()->particles_begin(
HepMC::parents) !=
481 mother = *((*partIter)->production_vertex()->particles_begin(
HepMC::parents));
483 if (((mother ==
nullptr) || ((mother !=
nullptr) && (mother->pdg_id() == 23)) ||
484 ((mother !=
nullptr) && (mother->pdg_id() == 32)) ||
485 ((mother !=
nullptr) && (
std::abs(mother->pdg_id()) == 24)))) {
487 pAssSim = genPc->momentum();
502 bool okSeedFound =
false;
503 double seedOkRatio = 999999.;
507 for (ElectronSeedCollection::const_iterator gsfIter = (*elSeeds).begin(); gsfIter != (*elSeeds).end();
510 for (
auto const &recHit : gsfIter->recHits()) {
511 det = pDD->idToDet(recHit.geographicalId());
519 double dphi = phi - pAssSim.phi();
521 dphi = dphi < 0 ? (CLHEP::twopi) + dphi : dphi - CLHEP::twopi;
527 double tmpSeedRatio = p / pAssSim.t();
529 seedOkRatio = tmpSeedRatio;
530 bestElectronSeed = *gsfIter;
544 mcEta[ip] = pAssSim.eta();
545 mcPhi[ip] = pAssSim.phi();
546 mcPt[ip] = pAssSim.perp();
547 mcQ[ip] = ((
id == 11) ? -1. : +1.);
553 seedOkRatio = 999999.;
556 for (ElectronSeedCollection::const_iterator gsfIter = (*elSeeds).begin(); gsfIter != (*elSeeds).end();
559 for (
auto const &recHit : gsfIter->recHits()) {
560 det = pDD->idToDet(recHit.geographicalId());
568 double dphi = phi - pAssSim.phi();
570 dphi = dphi < 0 ? (CLHEP::twopi) + dphi : dphi - CLHEP::twopi;
572 if (gsfIter->isEcalDriven()) {
577 double tmpSeedRatio = p / pAssSim.t();
579 seedOkRatio = tmpSeedRatio;
580 bestElectronSeed = *gsfIter;
597 seedOkRatio = 999999.;
600 for (ElectronSeedCollection::const_iterator gsfIter = (*elSeeds).begin(); gsfIter != (*elSeeds).end();
603 for (
auto const &recHit : gsfIter->recHits()) {
604 det = pDD->idToDet(recHit.geographicalId());
612 double dphi = phi - pAssSim.phi();
614 dphi = dphi < 0 ? (CLHEP::twopi) + dphi : dphi - CLHEP::twopi;
616 if (gsfIter->isTrackerDriven()) {
621 double tmpSeedRatio = p / pAssSim.t();
623 seedOkRatio = tmpSeedRatio;
624 bestElectronSeed = *gsfIter;
GlobalPoint toGlobal(const Point2DBase< Scalar, LocalTag > lp) const
float superclusterEnergy[10]
edm::ESGetToken< TrackerTopology, TrackerTopologyRcd > topoToken_
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
GlobalPoint globalPosition() const
TH1F * histtrackerdrivenetaMCmatched_
Exp< T >::type exp(const T &t)
const Plane & surface() const
The nominal surface of the GeomDet.
edm::InputTag inputCollection_
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
edm::ESGetToken< MagneticField, IdealMagneticFieldRecord > magFieldToken_
TH1F * histecaldriveneMCmatched_
Log< level::Info, false > LogInfo
float superclusterEta[10]
TH1F * histecaldrivenetaMCmatched_
T const * product() const
bool isNull() const
Checks for null.
TrajectoryStateOnSurface propagate(STA const &state, SUR const &surface) const
TH1F * histtrackerdriveneMCmatched_
GlobalVector globalMomentum() const
ESHandle< T > getHandle(const ESGetToken< T, R > &iToken) const
edm::ESGetToken< TrackerGeometry, TrackerDigiGeometryRecord > trackerToken_
Power< A, B >::type pow(const A &a, const B &b)
TH1F * histtrackerdrivenptMCmatched_
constexpr Detector det() const
get the detector field from this detid