85 #include "TDirectory.h"
95 class ParticlePtGreater {
98 return p1->momentum().perp() >
p2->momentum().perp();
102 class ParticlePGreater {
105 return p1->momentum().rho() >
p2->momentum().rho();
321 : ptMin_(iConfig.getUntrackedParameter<double>(
"PTMin", 1.0)),
322 etaMax_(iConfig.getUntrackedParameter<double>(
"MaxChargedHadronEta", 2.5)),
323 pCutIsolate_(iConfig.getUntrackedParameter<double>(
"PMaxIsolation", 20.0)),
324 a_Isolation_(iConfig.getUntrackedParameter<
bool>(
"UseConeIsolation",
false)),
325 genSrc_(iConfig.getUntrackedParameter(
"GenSrc",
std::
string(
"generatorSmeared"))),
326 useHepMC_(iConfig.getUntrackedParameter<
bool>(
"UseHepMC",
false)),
327 a_coneR_(iConfig.getUntrackedParameter<double>(
"ConeRadius", 34.98)),
328 a_mipR_(iConfig.getUntrackedParameter<double>(
"ConeRadiusMIP", 14.0)),
329 debugL1Info_(iConfig.getUntrackedParameter<
bool>(
"DebugL1Info",
false)),
330 verbosity_(iConfig.getUntrackedParameter<
int>(
"Verbosity", 0)) {
347 tok_L1GTrorsrc_ = consumes<L1GlobalTriggerReadoutRecord>(L1GTReadoutRcdSource_);
348 tok_L1GTobjmap_ = consumes<L1GlobalTriggerObjectMapRecord>(L1GTObjectMapRcdSource_);
349 tok_L1extMusrc_ = consumes<l1extra::L1MuonParticleCollection>(L1extraMuonSource_);
350 tok_L1Em_ = consumes<l1extra::L1EmParticleCollection>(L1extraIsoEmSource_);
351 tok_L1extNonIsoEm_ = consumes<l1extra::L1EmParticleCollection>(L1extraNonIsoEmSource_);
352 tok_L1extTauJet_ = consumes<l1extra::L1JetParticleCollection>(L1extraTauJetSource_);
353 tok_L1extCenJet_ = consumes<l1extra::L1JetParticleCollection>(L1extraCenJetSource_);
354 tok_L1extFwdJet_ = consumes<l1extra::L1JetParticleCollection>(L1extraFwdJetSource_);
356 if (!strcmp(
"Dummy",
genSrc_.c_str())) {
368 tok_geom_ = esConsumes<CaloGeometry, CaloGeometryRecord>();
370 tok_topo_ = esConsumes<HcalTopology, HcalRecNumberingRecord>();
371 tok_magField_ = esConsumes<MagneticField, IdealMagneticFieldRecord>();
372 tok_pdt_ = esConsumes<HepPDT::ParticleDataTable, PDTRecord>();
378 desc.addUntracked<
bool>(
"UseHepMC",
false);
379 desc.addUntracked<
double>(
"ChargedHadronSeedP", 1.0);
380 desc.addUntracked<
double>(
"PTMin", 1.0);
381 desc.addUntracked<
double>(
"MaxChargedHadronEta", 2.5);
382 desc.addUntracked<
double>(
"ConeRadius", 34.98);
383 desc.addUntracked<
double>(
"ConeRadiusMIP", 14.0);
384 desc.addUntracked<
bool>(
"UseConeIsolation",
true);
385 desc.addUntracked<
double>(
"PMaxIsolation", 5.0);
386 desc.addUntracked<
int>(
"Verbosity", 0);
387 desc.addUntracked<
bool>(
"DebugL1Info",
false);
396 descriptions.
add(
"isolatedGenParticles",
desc);
431 <<
"not found\n --> returning false by "
435 edm::LogVerbatim(
"IsoTrack") <<
"\nL1GlobalTriggerObjectMapRecord with \n\n"
436 <<
"not found\n --> returning false by "
442 unsigned int numberTriggerBits = dWord.size();
447 edm::LogVerbatim(
"IsoTrack") <<
"\nNumber of Trigger bits " << numberTriggerBits <<
"\n";
451 const std::vector<L1GlobalTriggerObjectMap> &objMapVec = gtOMRec->
gtObjectMap();
452 for (std::vector<L1GlobalTriggerObjectMap>::const_iterator itMap = objMapVec.begin(); itMap != objMapVec.end();
455 int itrig = (*itMap).algoBitNumber();
469 for (
unsigned int iBit = 0; iBit < numberTriggerBits; ++iBit) {
470 bool accept = dWord[iBit];
484 l1extra::L1JetParticleCollection::const_iterator itr;
485 for (itr = l1TauHandle->begin(); itr != l1TauHandle->end(); ++itr) {
490 edm::LogVerbatim(
"IsoTrack") <<
"tauJ p/pt " << itr->momentum() <<
" " << itr->pt() <<
" eta/phi " << itr->eta()
491 <<
" " << itr->phi();
498 for (itr = l1CenJetHandle->begin(); itr != l1CenJetHandle->end(); ++itr) {
503 edm::LogVerbatim(
"IsoTrack") <<
"cenJ p/pt " << itr->momentum() <<
" " << itr->pt() <<
" eta/phi "
504 << itr->eta() <<
" " << itr->phi();
510 for (itr = l1FwdJetHandle->begin(); itr != l1FwdJetHandle->end(); ++itr) {
515 edm::LogVerbatim(
"IsoTrack") <<
"fwdJ p/pt " << itr->momentum() <<
" " << itr->pt() <<
" eta/phi "
516 << itr->eta() <<
" " << itr->phi();
520 l1extra::L1EmParticleCollection::const_iterator itrEm;
523 for (itrEm = l1IsoEmHandle->begin(); itrEm != l1IsoEmHandle->end(); ++itrEm) {
528 edm::LogVerbatim(
"IsoTrack") <<
"isoEm p/pt " << itrEm->momentum() <<
" " << itrEm->pt() <<
" eta/phi "
529 << itrEm->eta() <<
" " << itrEm->phi();
535 for (itrEm = l1NonIsoEmHandle->begin(); itrEm != l1NonIsoEmHandle->end(); ++itrEm) {
540 edm::LogVerbatim(
"IsoTrack") <<
"nonIsoEm p/pt " << itrEm->momentum() <<
" " << itrEm->pt() <<
" eta/phi "
541 << itrEm->eta() <<
" " << itrEm->phi();
546 l1extra::L1MuonParticleCollection::const_iterator itrMu;
549 for (itrMu = l1MuHandle->begin(); itrMu != l1MuHandle->end(); ++itrMu) {
554 edm::LogVerbatim(
"IsoTrack") <<
"l1muon p/pt " << itrMu->momentum() <<
" " << itrMu->pt() <<
" eta/phi "
555 << itrMu->eta() <<
" " << itrMu->phi();
568 for (
unsigned int indx = 0; indx < trackIDs.size(); ++indx) {
569 int charge = trackIDs[indx].charge;
570 HepMC::GenEvent::particle_const_iterator
p = trackIDs[indx].trkItr;
572 (*p)->momentum().px(), (*p)->momentum().py(), (*p)->momentum().pz(), (*p)->momentum().e());
574 edm::LogVerbatim(
"IsoTrack") <<
"trkIndx " << indx <<
" pdgid " << trackIDs[indx].pdgId <<
" charge " <<
charge
575 <<
" momVec " << momVec;
580 posVec =
GlobalPoint(0.1 * (*p)->production_vertex()->position().x(),
581 0.1 * (*p)->production_vertex()->position().y(),
582 0.1 * (*p)->production_vertex()->position().z());
583 posECAL = trackIDs[indx].pointECAL;
584 fillTrack(posVec, momVec, posECAL, trackIDs[indx].
pdgId, trackIDs[indx].okECAL,
true);
586 edm::LogVerbatim(
"IsoTrack") <<
"posECAL " << posECAL <<
" okECAL " << trackIDs[indx].okECAL <<
"okHCAL "
587 << trackIDs[indx].okHCAL;
588 if (trackIDs[indx].okECAL) {
605 trackIDs[indx].directionECAL,
614 trackIDs[indx].directionECAL,
617 if (trackIDs[indx].okHCAL) {
628 trackIDs[indx].directionHCAL,
637 trackIDs[indx].directionHCAL,
642 bool saveTrack =
true;
652 fillTrack(posVec, momVec, posECAL, 0,
false,
false);
658 HepMC::GenEvent::particle_const_iterator
p;
659 for (
p = myGenEvent->particles_begin(), indx = 0;
p != myGenEvent->particles_end(); ++
p, ++indx) {
660 int pdgId = ((*p)->pdg_id());
663 double pp = (*p)->momentum().rho();
664 double eta = (*p)->momentum().eta();
669 std::vector<spr::propagatedGenParticleID> trackIDs =
672 for (
unsigned int indx = 0; indx < trackIDs.size(); ++indx) {
673 int charge = trackIDs[indx].charge;
674 reco::GenParticleCollection::const_iterator
p = trackIDs[indx].trkItr;
678 edm::LogVerbatim(
"IsoTrack") <<
"trkIndx " << indx <<
" pdgid " << trackIDs[indx].pdgId <<
" charge " <<
charge
679 <<
" momVec " << momVec;
684 edm::LogVerbatim(
"IsoTrack") <<
" pt " << momVec.Pt() <<
" eta " << momVec.eta();
686 posVec =
GlobalPoint(
p->vertex().x(),
p->vertex().y(),
p->vertex().z());
687 posECAL = trackIDs[indx].pointECAL;
689 edm::LogVerbatim(
"IsoTrack") <<
"posECAL " << posECAL <<
" okECAL " << trackIDs[indx].okECAL <<
"okHCAL "
690 << trackIDs[indx].okHCAL;
691 fillTrack(posVec, momVec, posECAL, trackIDs[indx].
pdgId, trackIDs[indx].okECAL,
true);
692 if (trackIDs[indx].okECAL) {
718 trackIDs[indx].directionECAL,
727 trackIDs[indx].directionECAL,
730 if (trackIDs[indx].okHCAL) {
745 trackIDs[indx].directionHCAL,
754 trackIDs[indx].directionHCAL,
759 bool saveTrack =
true;
769 fillTrack(posVec, momVec, posECAL, 0,
false,
false);
775 reco::GenParticleCollection::const_iterator
p;
780 double pp = (
p->momentum()).
R();
781 double eta = (
p->momentum()).Eta();
797 double tempgen_TH[
NPBins_ + 1] = {0.0, 5.0, 12.0, 300.0};
801 double tempgen_Eta[
NEtaBins_ + 1] = {0.0, 0.5, 1.1, 1.7, 2.3};
817 double phi1 = momVec.phi();
818 double phi2 = (posECAL - posVec).
phi();
820 double deta = momVec.eta() - (posECAL - posVec).
eta();
1018 h_NEventProc = fs->
make<TH1I>(
"h_NEventProc",
"h_NEventProc", 2, -0.5, 0.5);
1019 h_L1AlgoNames = fs->
make<TH1I>(
"h_L1AlgoNames",
"h_L1AlgoNames:Bin Labels", 128, -0.5, 127.5);
1021 double pBin[
PBins_ + 1] = {0.0, 2.0, 4.0, 6.0, 8.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0,
1022 70.0, 80.0, 90.0, 100.0, 150.0, 200.0, 250.0, 300.0, 350.0, 400.0, 450.0,
1023 500.0, 550.0, 600.0, 650.0, 700.0, 750.0, 800.0, 850.0, 900.0, 950.0, 1000.0};
1024 double etaBin[
EtaBins_ + 1] = {-3.0, -2.9, -2.8, -2.7, -2.6, -2.5, -2.4, -2.3, -2.2, -2.1, -2.0, -1.9, -1.8,
1025 -1.7, -1.6, -1.5, -1.4, -1.3, -1.2, -1.1, -1.0, -0.9, -0.8, -0.7, -0.6, -0.5,
1026 -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8,
1027 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1,
1028 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3.0};
1030 "electron",
"positron",
"#gamma",
"#pi^+",
"#pi^-",
"K^+",
"K^-",
"p",
"n",
"pbar",
"nbar",
"K^0_L"};
1034 sprintf(
name,
"pEta%d",
i);
1035 sprintf(
title,
"#eta vs momentum for %s", particle[
i].c_str());
1040 tree_ = fs->
make<TTree>(
"tree_",
"tree");
1197 t_muEneR =
new std::vector<double>();
1410 tree_->Branch(
"t_muEneR",
"std::vector<double>", &
t_muEneR);
1682 int partID[
Particles] = {11, -11, 21, 211, -211, 321, -321, 2212, 2112, -2212, -2112, 130};
1684 for (
int ik = 0; ik <
Particles; ++ik) {
1685 if (
pdgId == partID[ik]) {