255 edm::LogInfo(
"") <<
"\n\n =================> Treating event " << e.
id() <<
" Number of seeds "
259 float mass = .000511;
263 float dphi1 = 0., dphi2 = 0., drz1 = 0., drz2 = 0.;
264 float phi1 = 0., phi2 = 0., rz1 = 0., rz2 = 0.;
266 for (ElectronSeedCollection::const_iterator MyS = (*elSeeds).begin(); MyS != (*elSeeds).end(); ++MyS) {
267 LogDebug(
"") <<
"\nSeed nr " << is <<
": ";
269 LogDebug(
"") <<
" Number of RecHits= " << (*MyS).nHits();
274 DetId id1 = (*it).geographicalId();
275 det1 = pDD->idToDet(id1);
276 LogDebug(
"") <<
" First hit local x,y,z " << (*it).localPosition() <<
" det " << id1.
det() <<
" subdet "
278 LogDebug(
"") <<
" First hit global " << det1->
toGlobal((*it).localPosition());
282 DetId id2 = (*it).geographicalId();
283 det2 = pDD->idToDet(id2);
284 LogDebug(
"") <<
" Second hit local x,y,z " << (*it).localPosition() <<
" det " << id2.
det() <<
" subdet "
286 LogDebug(
"") <<
" Second hit global " << det2->
toGlobal((*it).localPosition());
292 for (
auto const &recHit : r) {
293 det = pDD->idToDet(recHit.geographicalId());
306 LogDebug(
"") <<
" ElectronSeed superCluster energy: " << theClus->energy() <<
", position: " << theClus->position();
308 LogDebug(
"") <<
" ElectronSeed supercluster Et: "
309 << theClus->energy() *
sin(2. * atan(
exp(-theClus->position().eta())));
311 LogDebug(
"") <<
" ElectronSeed supercluster eta: " << theClus->position().eta();
312 LogDebug(
"") <<
" ElectronSeed seed charge: " << (*MyS).getCharge();
319 int charge = int((*MyS).getCharge());
320 GlobalPoint xmeas(theClus->position().x(), theClus->position().y(), theClus->position().z());
321 GlobalPoint vprim(theBeamSpot->position().x(), theBeamSpot->position().y(), theBeamSpot->position().z());
322 float energy = theClus->energy();
336 DetId id = (*it).geographicalId();
337 const GeomDet *geomdet = pDD->idToDet((*it).geographicalId());
346 float SCl_phi = xmeas.phi();
347 float localDphi = SCl_phi - hitPos.
phi();
359 if (
id.subdetId() % 2 == 1) {
366 DetId id2 = (*it).geographicalId();
367 const GeomDet *geomdet2 = pDD->idToDet((*it).geographicalId());
371 double pxHit1z = hitPos.
z();
372 double pxHit1x = hitPos.
x();
373 double pxHit1y = hitPos.
y();
374 double r1diff = (pxHit1x - vprim.x()) * (pxHit1x - vprim.x()) + (pxHit1y - vprim.y()) * (pxHit1y - vprim.y());
375 r1diff =
sqrt(r1diff);
376 double r2diff = (xmeas.x() - pxHit1x) * (xmeas.x() - pxHit1x) + (xmeas.y() - pxHit1y) * (xmeas.y() - pxHit1y);
377 r2diff =
sqrt(r2diff);
378 double zVertexPred = pxHit1z - r1diff * (xmeas.z() - pxHit1z) / r2diff;
380 GlobalPoint vertexPred(vprim.x(), vprim.y(), zVertexPred);
388 phi2 = hitPos2.
phi();
390 rz2 = hitPos2.
perp();
402 histetclu_->Fill(theClus->energy() *
sin(2. * atan(
exp(-theClus->position().eta()))));
405 histq_->Fill((*MyS).getCharge());
417 seedQ[is] = (*MyS).getCharge();
443 e.
getByLabel(
"correctedHybridSuperClusters", clusters);
445 if (!clusters.
product()->empty())
453 e.
getByLabel(
"generatorSmeared",
"", HepMCEvt);
457 HepMC::FourVector pAssSim;
459 for (HepMC::GenEvent::particle_const_iterator partIter = MCEvt->particles_begin(); partIter != MCEvt->particles_end();
461 for (HepMC::GenEvent::vertex_const_iterator vertIter = MCEvt->vertices_begin(); vertIter != MCEvt->vertices_end();
464 HepMC::GenVertex *creation = (*partIter)->production_vertex();
466 HepMC::FourVector momentum = (*partIter)->momentum();
468 int id = (*partIter)->pdg_id();
469 LogDebug(
"") <<
"MC particle id " <<
id <<
", creationVertex " << (*creation) <<
" cm, initialMomentum "
470 << momentum.rho() <<
" GeV/c" << std::endl;
472 if (
id == 11 ||
id == -11) {
475 if ((*partIter)->production_vertex()) {
476 if ((*partIter)->production_vertex()->particles_begin(
HepMC::parents) !=
478 mother = *((*partIter)->production_vertex()->particles_begin(
HepMC::parents));
480 if (((mother ==
nullptr) || ((mother !=
nullptr) && (mother->pdg_id() == 23)) ||
481 ((mother !=
nullptr) && (mother->pdg_id() == 32)) ||
482 ((mother !=
nullptr) && (
std::abs(mother->pdg_id()) == 24)))) {
484 pAssSim = genPc->momentum();
499 bool okSeedFound =
false;
500 double seedOkRatio = 999999.;
504 for (ElectronSeedCollection::const_iterator gsfIter = (*elSeeds).begin(); gsfIter != (*elSeeds).end();
507 for (
auto const &recHit : gsfIter->recHits()) {
508 det = pDD->idToDet(recHit.geographicalId());
516 double dphi = phi - pAssSim.phi();
518 dphi = dphi < 0 ? (CLHEP::twopi) + dphi : dphi - CLHEP::twopi;
524 double tmpSeedRatio = p / pAssSim.t();
526 seedOkRatio = tmpSeedRatio;
527 bestElectronSeed = *gsfIter;
541 mcEta[ip] = pAssSim.eta();
542 mcPhi[ip] = pAssSim.phi();
543 mcPt[ip] = pAssSim.perp();
544 mcQ[ip] = ((
id == 11) ? -1. : +1.);
550 seedOkRatio = 999999.;
553 for (ElectronSeedCollection::const_iterator gsfIter = (*elSeeds).begin(); gsfIter != (*elSeeds).end();
556 for (
auto const &recHit : gsfIter->recHits()) {
557 det = pDD->idToDet(recHit.geographicalId());
565 double dphi = phi - pAssSim.phi();
567 dphi = dphi < 0 ? (CLHEP::twopi) + dphi : dphi - CLHEP::twopi;
569 if (gsfIter->isEcalDriven()) {
574 double tmpSeedRatio = p / pAssSim.t();
576 seedOkRatio = tmpSeedRatio;
577 bestElectronSeed = *gsfIter;
594 seedOkRatio = 999999.;
597 for (ElectronSeedCollection::const_iterator gsfIter = (*elSeeds).begin(); gsfIter != (*elSeeds).end();
600 for (
auto const &recHit : gsfIter->recHits()) {
601 det = pDD->idToDet(recHit.geographicalId());
609 double dphi = phi - pAssSim.phi();
611 dphi = dphi < 0 ? (CLHEP::twopi) + dphi : dphi - CLHEP::twopi;
613 if (gsfIter->isTrackerDriven()) {
618 double tmpSeedRatio = p / pAssSim.t();
620 seedOkRatio = tmpSeedRatio;
621 bestElectronSeed = *gsfIter;
GlobalPoint toGlobal(const Point2DBase< Scalar, LocalTag > lp) const
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
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
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
Power< A, B >::type pow(const A &a, const B &b)
TH1F * histtrackerdrivenptMCmatched_
constexpr Detector det() const
get the detector field from this detid