86 #include "TDirectory.h"
96 class ParticlePtGreater {
99 return p1->momentum().perp() >
p2->momentum().perp();
103 class ParticlePGreater {
106 return p1->momentum().rho() >
p2->momentum().rho();
316 : ptMin_(iConfig.getUntrackedParameter<double>(
"PTMin", 1.0)),
317 etaMax_(iConfig.getUntrackedParameter<double>(
"MaxChargedHadronEta", 2.5)),
318 pCutIsolate_(iConfig.getUntrackedParameter<double>(
"PMaxIsolation", 20.0)),
319 a_Isolation_(iConfig.getUntrackedParameter<
bool>(
"UseConeIsolation",
false)),
320 genSrc_(iConfig.getUntrackedParameter(
"GenSrc",
std::
string(
"generatorSmeared"))),
321 useHepMC_(iConfig.getUntrackedParameter<
bool>(
"UseHepMC",
false)),
322 a_coneR_(iConfig.getUntrackedParameter<double>(
"ConeRadius", 34.98)),
323 a_mipR_(iConfig.getUntrackedParameter<double>(
"ConeRadiusMIP", 14.0)),
324 debugL1Info_(iConfig.getUntrackedParameter<
bool>(
"DebugL1Info",
false)),
325 verbosity_(iConfig.getUntrackedParameter<
int>(
"Verbosity", 0)) {
342 tok_L1GTrorsrc_ = consumes<L1GlobalTriggerReadoutRecord>(L1GTReadoutRcdSource_);
343 tok_L1GTobjmap_ = consumes<L1GlobalTriggerObjectMapRecord>(L1GTObjectMapRcdSource_);
344 tok_L1extMusrc_ = consumes<l1extra::L1MuonParticleCollection>(L1extraMuonSource_);
345 tok_L1Em_ = consumes<l1extra::L1EmParticleCollection>(L1extraIsoEmSource_);
346 tok_L1extNonIsoEm_ = consumes<l1extra::L1EmParticleCollection>(L1extraNonIsoEmSource_);
347 tok_L1extTauJet_ = consumes<l1extra::L1JetParticleCollection>(L1extraTauJetSource_);
348 tok_L1extCenJet_ = consumes<l1extra::L1JetParticleCollection>(L1extraCenJetSource_);
349 tok_L1extFwdJet_ = consumes<l1extra::L1JetParticleCollection>(L1extraFwdJetSource_);
351 if (!strcmp(
"Dummy",
genSrc_.c_str())) {
367 desc.addUntracked<
bool>(
"UseHepMC",
false);
368 desc.addUntracked<
double>(
"ChargedHadronSeedP", 1.0);
369 desc.addUntracked<
double>(
"PTMin", 1.0);
370 desc.addUntracked<
double>(
"MaxChargedHadronEta", 2.5);
371 desc.addUntracked<
double>(
"ConeRadius", 34.98);
372 desc.addUntracked<
double>(
"ConeRadiusMIP", 14.0);
373 desc.addUntracked<
bool>(
"UseConeIsolation",
true);
374 desc.addUntracked<
double>(
"PMaxIsolation", 5.0);
375 desc.addUntracked<
int>(
"Verbosity", 0);
376 desc.addUntracked<
bool>(
"DebugL1Info",
false);
385 descriptions.
add(
"isolatedGenParticles",
desc);
432 <<
"not found\n --> returning false by "
436 edm::LogVerbatim(
"IsoTrack") <<
"\nL1GlobalTriggerObjectMapRecord with \n\n"
437 <<
"not found\n --> returning false by "
443 unsigned int numberTriggerBits = dWord.size();
448 edm::LogVerbatim(
"IsoTrack") <<
"\nNumber of Trigger bits " << numberTriggerBits <<
"\n";
452 const std::vector<L1GlobalTriggerObjectMap> &objMapVec = gtOMRec->
gtObjectMap();
453 for (std::vector<L1GlobalTriggerObjectMap>::const_iterator itMap = objMapVec.begin(); itMap != objMapVec.end();
456 int itrig = (*itMap).algoBitNumber();
470 for (
unsigned int iBit = 0; iBit < numberTriggerBits; ++iBit) {
471 bool accept = dWord[iBit];
485 l1extra::L1JetParticleCollection::const_iterator
itr;
486 for (
itr = l1TauHandle->begin();
itr != l1TauHandle->end(); ++
itr) {
492 <<
" " <<
itr->phi();
499 for (
itr = l1CenJetHandle->begin();
itr != l1CenJetHandle->end(); ++
itr) {
505 <<
itr->eta() <<
" " <<
itr->phi();
511 for (
itr = l1FwdJetHandle->begin();
itr != l1FwdJetHandle->end(); ++
itr) {
517 <<
itr->eta() <<
" " <<
itr->phi();
521 l1extra::L1EmParticleCollection::const_iterator itrEm;
524 for (itrEm = l1IsoEmHandle->begin(); itrEm != l1IsoEmHandle->end(); ++itrEm) {
529 edm::LogVerbatim(
"IsoTrack") <<
"isoEm p/pt " << itrEm->momentum() <<
" " << itrEm->pt() <<
" eta/phi "
530 << itrEm->eta() <<
" " << itrEm->phi();
536 for (itrEm = l1NonIsoEmHandle->begin(); itrEm != l1NonIsoEmHandle->end(); ++itrEm) {
541 edm::LogVerbatim(
"IsoTrack") <<
"nonIsoEm p/pt " << itrEm->momentum() <<
" " << itrEm->pt() <<
" eta/phi "
542 << itrEm->eta() <<
" " << itrEm->phi();
547 l1extra::L1MuonParticleCollection::const_iterator itrMu;
550 for (itrMu = l1MuHandle->begin(); itrMu != l1MuHandle->end(); ++itrMu) {
555 edm::LogVerbatim(
"IsoTrack") <<
"l1muon p/pt " << itrMu->momentum() <<
" " << itrMu->pt() <<
" eta/phi "
556 << itrMu->eta() <<
" " << itrMu->phi();
569 for (
unsigned int indx = 0; indx < trackIDs.size(); ++indx) {
570 int charge = trackIDs[indx].charge;
571 HepMC::GenEvent::particle_const_iterator
p = trackIDs[indx].trkItr;
573 (*p)->momentum().px(), (*p)->momentum().py(), (*p)->momentum().pz(), (*p)->momentum().e());
575 edm::LogVerbatim(
"IsoTrack") <<
"trkIndx " << indx <<
" pdgid " << trackIDs[indx].pdgId <<
" charge " <<
charge
576 <<
" momVec " << momVec;
581 posVec =
GlobalPoint(0.1 * (*p)->production_vertex()->position().x(),
582 0.1 * (*p)->production_vertex()->position().y(),
583 0.1 * (*p)->production_vertex()->position().z());
584 posECAL = trackIDs[indx].pointECAL;
585 fillTrack(posVec, momVec, posECAL, trackIDs[indx].
pdgId, trackIDs[indx].okECAL,
true);
587 edm::LogVerbatim(
"IsoTrack") <<
"posECAL " << posECAL <<
" okECAL " << trackIDs[indx].okECAL <<
"okHCAL "
588 << trackIDs[indx].okHCAL;
589 if (trackIDs[indx].okECAL) {
606 trackIDs[indx].directionECAL,
615 trackIDs[indx].directionECAL,
618 if (trackIDs[indx].okHCAL) {
629 trackIDs[indx].directionHCAL,
638 trackIDs[indx].directionHCAL,
643 bool saveTrack =
true;
653 fillTrack(posVec, momVec, posECAL, 0,
false,
false);
659 HepMC::GenEvent::particle_const_iterator
p;
660 for (
p = myGenEvent->particles_begin(), indx = 0;
p != myGenEvent->particles_end(); ++
p, ++indx) {
661 int pdgId = ((*p)->pdg_id());
664 double pp = (*p)->momentum().rho();
665 double eta = (*p)->momentum().eta();
670 std::vector<spr::propagatedGenParticleID> trackIDs =
673 for (
unsigned int indx = 0; indx < trackIDs.size(); ++indx) {
674 int charge = trackIDs[indx].charge;
675 reco::GenParticleCollection::const_iterator
p = trackIDs[indx].trkItr;
679 edm::LogVerbatim(
"IsoTrack") <<
"trkIndx " << indx <<
" pdgid " << trackIDs[indx].pdgId <<
" charge " <<
charge
680 <<
" momVec " << momVec;
685 edm::LogVerbatim(
"IsoTrack") <<
" pt " << momVec.Pt() <<
" eta " << momVec.eta();
687 posVec =
GlobalPoint(
p->vertex().x(),
p->vertex().y(),
p->vertex().z());
688 posECAL = trackIDs[indx].pointECAL;
690 edm::LogVerbatim(
"IsoTrack") <<
"posECAL " << posECAL <<
" okECAL " << trackIDs[indx].okECAL <<
"okHCAL "
691 << trackIDs[indx].okHCAL;
692 fillTrack(posVec, momVec, posECAL, trackIDs[indx].
pdgId, trackIDs[indx].okECAL,
true);
693 if (trackIDs[indx].okECAL) {
719 trackIDs[indx].directionECAL,
728 trackIDs[indx].directionECAL,
731 if (trackIDs[indx].okHCAL) {
746 trackIDs[indx].directionHCAL,
755 trackIDs[indx].directionHCAL,
760 bool saveTrack =
true;
770 fillTrack(posVec, momVec, posECAL, 0,
false,
false);
776 reco::GenParticleCollection::const_iterator
p;
781 double pp = (
p->momentum()).
R();
782 double eta = (
p->momentum()).Eta();
798 double tempgen_TH[
NPBins_ + 1] = {0.0, 5.0, 12.0, 300.0};
802 double tempgen_Eta[
NEtaBins_ + 1] = {0.0, 0.5, 1.1, 1.7, 2.3};
818 double phi1 = momVec.phi();
819 double phi2 = (posECAL - posVec).
phi();
821 double deta = momVec.eta() - (posECAL - posVec).
eta();
1019 h_NEventProc = fs->
make<TH1I>(
"h_NEventProc",
"h_NEventProc", 2, -0.5, 0.5);
1020 h_L1AlgoNames = fs->
make<TH1I>(
"h_L1AlgoNames",
"h_L1AlgoNames:Bin Labels", 128, -0.5, 127.5);
1022 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,
1023 70.0, 80.0, 90.0, 100.0, 150.0, 200.0, 250.0, 300.0, 350.0, 400.0, 450.0,
1024 500.0, 550.0, 600.0, 650.0, 700.0, 750.0, 800.0, 850.0, 900.0, 950.0, 1000.0};
1025 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,
1026 -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,
1027 -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,
1028 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,
1029 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3.0};
1031 "electron",
"positron",
"#gamma",
"#pi^+",
"#pi^-",
"K^+",
"K^-",
"p",
"n",
"pbar",
"nbar",
"K^0_L"};
1035 sprintf(
name,
"pEta%d",
i);
1036 sprintf(
title,
"#eta vs momentum for %s", particle[
i].c_str());
1041 tree_ = fs->
make<TTree>(
"tree_",
"tree");
1198 t_muEneR =
new std::vector<double>();
1411 tree_->Branch(
"t_muEneR",
"std::vector<double>", &
t_muEneR);
1683 int partID[
Particles] = {11, -11, 21, 211, -211, 321, -321, 2212, 2112, -2212, -2112, 130};
1685 for (
int ik = 0; ik <
Particles; ++ik) {
1686 if (
pdgId == partID[ik]) {