9 #include "TDirectory.h" 11 #include "TLorentzVector.h" 12 #include "TInterpreter.h" 100 double deltaR(
double eta1,
double eta2,
double phi1,
double phi2);
142 std::map< std::pair<unsigned int,std::string>,
int>
l1AlgoMap_;
148 verbosity_(iConfig.getUntrackedParameter<
int>(
"Verbosity",0)),
151 a_coneR_(iConfig.getUntrackedParameter<double>(
"ConeRadius",34.98)),
153 a_mipR_(iConfig.getUntrackedParameter<double>(
"ConeRadiusMIP",14.0)),
190 tok_L1extTauJet_ = consumes<l1extra::L1JetParticleCollection>(L1extraTauJetSource_);
191 tok_L1extCenJet_ = consumes<l1extra::L1JetParticleCollection>(L1extraCenJetSource_);
192 tok_L1extFwdJet_ = consumes<l1extra::L1JetParticleCollection>(L1extraFwdJetSource_);
195 <<
"Parameters read from config file \n" 210 <<
"\t isItAOD " << isItAOD;
227 std::vector<std::string> seeds = {
"L1_SingleJet36",
"L1_SingleJet52",
228 "L1_SingleJet68",
"L1_SingleJet92",
232 desc.
addUntracked<std::vector<std::string> >(
"L1Seed", seeds);
246 descriptions.
add(
"isoTrackCalib",desc);
256 <<
"Run " <<
t_Run <<
" Event " <<
t_Event <<
" Luminosity " 258 <<
" starts ==========";
287 for (
unsigned iGenJet = 0; iGenJet < genJets->size(); ++iGenJet) {
289 double genJetPt = genJet.
pt();
290 double genJetEta = genJet.
eta();
293 if (genJetEta>-2.5 && genJetEta<2.5) {
311 if (!recVtxs->empty() && !((*recVtxs)[0].isFake())) {
312 leadPV =
math::XYZPoint( (*recVtxs)[0].
x(),(*recVtxs)[0].
y(), (*recVtxs)[0].
z() );
313 }
else if (beamSpotH.
isValid()) {
318 edm::LogInfo(
"IsoTrack") <<
"Primary Vertex " << leadPV;
332 for (rhitItr=hbhe->
begin();rhitItr!=hbhe->
end();rhitItr++) {
333 double rec_energy = rhitItr->energy();
334 int rec_ieta = rhitItr->id().ieta();
335 int rec_depth = rhitItr->id().depth();
336 int rec_zside = rhitItr->id().zside();
337 double num1_1 = rec_zside*(rec_ieta+0.2*(rec_depth-1));
339 edm::LogInfo(
"IsoTrack") <<
"detid/rechit/ieta/zside/depth/num " <<
" = " 340 << rhitItr->id() <<
"/" << rec_energy <<
"/" 341 << rec_ieta <<
"/" << rec_zside <<
"/" 342 << rec_depth <<
"/" << num1_1;
350 std::vector<spr::propagatedTrackDirection> trkCaloDirections;
353 std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
354 for (trkDetItr = trkCaloDirections.begin();
355 trkDetItr != trkCaloDirections.end(); trkDetItr++) {
356 if (trkDetItr->okHCAL) {
358 int tk_ieta = detId.
ieta();
359 const reco::Track* pTrack = &(*(trkDetItr->trkItr));
360 double tk_p = pTrack->
p();
362 for (
unsigned int k=1;
k<
pbin.size(); ++
k) {
381 edm::LogInfo(
"IsoTrack") <<
"\nL1 configuration code:" << l1ConfCode
382 <<
"\nNo valid L1 trigger configuration available." 383 <<
"\nSee text above for error code interpretation" 384 <<
"\nNo return here, in order to test each method" 385 <<
", protected against configuration error.";
391 edm::LogInfo(
"IsoTrack") <<
"menuName " << menuName << std::endl;
393 std::vector<int> algbits;
394 for (
CItAlgo itAlgo = algorithmMap.begin(); itAlgo != algorithmMap.end();
397 int algBitNumber = ( itAlgo->second ).algoBitNumber();
402 edm::LogInfo(
"IsoTrack") << algName <<
" "<< algBitNumber <<
" " 405 if (algName.find(
l1Names_[
i])!=std::string::npos){
407 edm::LogInfo(
"IsoTrack") <<
"match found" <<
" " << algName <<
" " 410 if (decision > 0) l1ok =
true;
418 l1extra::L1JetParticleCollection::const_iterator itr;
419 double ptTriggered = -10;
420 double etaTriggered = -100;
421 double phiTriggered = -100;
423 for(itr = l1TauHandle->begin(); itr != l1TauHandle->end(); ++itr) {
424 if(itr->pt()>ptTriggered){
425 ptTriggered = itr->pt();
426 etaTriggered = itr->eta();
427 phiTriggered = itr->phi() ;
430 edm::LogInfo(
"IsoTrack") <<
"tauJ pt " << itr->pt() <<
" eta/phi " 431 << itr->eta() <<
" " << itr->phi();
435 for (itr = l1CenJetHandle->begin(); itr != l1CenJetHandle->end(); ++itr){
436 if (itr->pt()>ptTriggered) {
437 ptTriggered = itr->pt();
438 etaTriggered = itr->eta();
439 phiTriggered = itr->phi() ;
443 <<
" eta/phi " << itr->eta() <<
" " 448 for (itr = l1FwdJetHandle->begin(); itr != l1FwdJetHandle->end(); ++itr) {
449 if (itr->pt()>ptTriggered) {
450 ptTriggered = itr->pt();
451 etaTriggered = itr->eta();
452 phiTriggered = itr->phi() ;
455 edm::LogInfo(
"IsoTrack") <<
"forJ pt " << itr->pt() <<
" eta/phi " 456 << itr->eta() <<
" " << itr->phi();
459 edm::LogInfo(
"IsoTrack") <<
"jets pt/eta/phi = " << ptTriggered <<
"/" 460 << etaTriggered <<
"/" <<phiTriggered;
462 unsigned int nTracks(0),nselTracks(0);
463 for (trkDetItr = trkCaloDirections.begin(),
nTracks=0;
464 trkDetItr != trkCaloDirections.end(); trkDetItr++,
nTracks++) {
465 const reco::Track* pTrack = &(*(trkDetItr->trkItr));
467 pTrack->
pz(), pTrack->
p());
473 << pTrack->
pt() <<
"/" << pTrack->
eta() <<
"/" 474 << pTrack->
phi() <<
"/" << pTrack->
p();
489 oneCutParameters.
maxDzPV = 100;
495 oneCutParameters.
maxDzPV = 100;
502 if (trkDetItr->okHCAL) {
507 edm::LogInfo(
"IsoTrack") <<
"qltyFlag|okECAL|okHCAL : " << qltyFlag <<
"|" 508 << trkDetItr->okECAL <<
"/" << trkDetItr->okHCAL;
509 t_qltyFlag = (qltyFlag && trkDetItr->okECAL && trkDetItr->okHCAL);
512 for (
unsigned int k=1;
k<
pbin.size(); ++
k) {
521 for (
unsigned int k=1;
k<
pbin.size(); ++
k) {
527 int nRH_eMipDR(0), nNearTRKs(0), nRecHits(-999);
530 trkDetItr->pointHCAL,
531 trkDetItr->pointECAL,
532 a_mipR_, trkDetItr->directionECAL,
535 std::vector<DetId>
ids;
538 trkDetItr->directionHCAL,nRecHits,
540 for (
unsigned int k=0;
k<ids.size(); ++
k) {
548 for (
unsigned int k=1;
k<
pbin.size(); ++
k) {
556 for (
unsigned int k=1;
k<
pbin.size(); ++
k) {
564 for (
unsigned int k=1;
k<
pbin.size(); ++
k) {
575 << pTrack->
pt() <<
"/" << pTrack->
eta() <<
"/" 576 << pTrack->
phi() <<
"/" <<
t_p;
581 for (
unsigned int lll=0;lll<
t_DetIds->size();lll++) {
601 double prange[5] = {20,30,40,60,100};
602 for (
int k=0;
k<5; ++
k)
pbin.push_back(prange[
k]);
604 "Charge Isolation",
"MIP Cut",
"L1 Cut"};
605 for (
unsigned int k=0; k<
pbin.size(); ++
k) {
607 if (k == 0) sprintf (namp,
"all momentum");
608 else sprintf (namp,
"p = %4.0f:%4.0f GeV",
pbin[k-1],
pbin[k]);
609 sprintf (name,
"TrackEta0%d", k);
610 sprintf (title,
"Track #eta for tracks with %s (%s)",namp,type[0].c_str());
612 sprintf (name,
"TrackEta1%d", k);
613 sprintf (title,
"Track #eta for tracks with %s (%s)",namp,type[1].c_str());
615 sprintf (name,
"TrackEta2%d", k);
616 sprintf (title,
"Track #eta for tracks with %s (%s)",namp,type[2].c_str());
618 sprintf (name,
"TrackEta3%d", k);
619 sprintf (title,
"Track #eta for tracks with %s (%s)",namp,type[3].c_str());
621 sprintf (name,
"TrackEta4%d", k);
622 sprintf (title,
"Track #eta for tracks with %s (%s)",namp,type[4].c_str());
624 sprintf (name,
"TrackEta5%d", k);
625 sprintf (title,
"Track #eta for tracks with %s (%s)",namp,type[5].c_str());
628 h_jetpt[0] =
fs_->
make<TH1F>(
"Jetpt0",
"Jet p_T (All)", 500,0.,2500.);
629 h_jetpt[1] =
fs_->
make<TH1F>(
"Jetpt1",
"Jet p_T (All Weighted)", 500,0.,2500.);
630 h_jetpt[2] =
fs_->
make<TH1F>(
"Jetpt2",
"Jet p_T (|#eta| < 2.5)", 500,0.,2500.);
631 h_jetpt[3] =
fs_->
make<TH1F>(
"Jetpt3",
"Jet p_T (|#eta| < 2.5 Weighted)", 500,0.,2500.);
635 tree->Branch(
"t_Run", &
t_Run,
"t_Run/I");
645 tree->Branch(
"t_p", &
t_p,
"t_p/D");
656 t_DetIds =
new std::vector<unsigned int>();
659 tree->Branch(
"t_DetIds",
"std::vector<unsigned int>", &
t_DetIds);
660 tree->Branch(
"t_HitEnergies",
"std::vector<double>", &t_HitEnergies);
661 tree->Branch(
"t_l1bits",
"std::vector<bool>", &t_l1bits);
666 bool changed =
false;
676 for (
CItAlgo itAlgo = algorithmMap.begin(); itAlgo != algorithmMap.end();
679 int algBitNumber = ( itAlgo->second ).algoBitNumber();
680 l1AlgoMap_.insert( std::pair<std::pair<unsigned int,std::string>,
int>( std::pair<unsigned int,std::string>(algBitNumber, algName) , 0) ) ;
682 std::map< std::pair<unsigned int,std::string>,
int>::iterator itr;
686 << (itr->first).second <<
" "<<itr->second;
696 return (vec1.eta()-vec2.eta());
701 double phi1 = vec1.phi();
702 if (phi1 < 0) phi1 += 2.0*
M_PI;
703 double phi2 = vec2.phi();
704 if (phi2 < 0) phi2 += 2.0*
M_PI;
705 double dphi = phi1-phi2;
707 else if (dphi < -
M_PI) dphi += 2.*
M_PI;
712 double deta =
dEta(vec1,vec2);
713 double dphi =
dPhi(vec1,vec2);
718 double deta = eta1-eta2;
719 if (phi1 < 0) phi1 += 2.0*
M_PI;
720 if (phi2 < 0) phi2 += 2.0*
M_PI;
721 double dphi = phi1-phi2;
723 else if (dphi < -
M_PI) dphi += 2.*
M_PI;
static const std::string kSharedResource
double p() const
momentum vector magnitude
void beginRun(edm::Run const &, edm::EventSetup const &) override
double eCone_hcal(const CaloGeometry *geo, edm::Handle< T > &hits, const GlobalPoint &hpoint1, const GlobalPoint &point1, double dR, const GlobalVector &trackMom, int &nRecHits, double hbThr=-100, double heThr=-100, double hfThr=-100, double hoThr=-100, double tMin=-500, double tMax=500, int detOnly=-1, int useRaw=0, bool debug=false)
T getParameter(std::string const &) const
EventNumber_t event() const
spr::trackSelectionParameters selectionParameters_
std::vector< double > * t_HitEnergies
T getUntrackedParameter(std::string const &, T const &) const
const edm::EDGetTokenT< reco::TrackCollection > tok_genTrack_
std::vector< spr::propagatedTrackID > propagateCALO(edm::Handle< reco::TrackCollection > &trkCollection, const CaloGeometry *geo, const MagneticField *bField, const std::string &theTrackQuality, bool debug=false)
const unsigned int nTracks(const reco::Vertex &sv)
const L1GtTriggerMenu * m_l1GtMenu
const edm::EDGetTokenT< GenEventInfoProduct > tok_ew_
std::vector< unsigned int > * t_DetIds
double eta() const final
momentum pseudorapidity
edm::EDGetTokenT< EcalRecHitCollection > tok_EE_
ParameterDescriptionBase * addUntracked(U const &iLabel, T const &value)
HLTConfigProvider hltConfig_
edm::EDGetTokenT< l1extra::L1JetParticleCollection > tok_L1extCenJet_
std::map< std::pair< unsigned int, std::string >, int > l1AlgoMap_
bool getByToken(EDGetToken token, Handle< PROD > &result) const
TrackQuality
track quality
const bool availableL1Configuration(int &errorCode, int &l1ConfCode) const
const edm::EDGetTokenT< reco::BeamSpot > tok_bs_
std::vector< GenJet > GenJetCollection
collection of GenJet objects
std::vector< Track > TrackCollection
collection of Tracks
std::vector< T >::const_iterator const_iterator
int bunchCrossing() const
edm::LuminosityBlockNumber_t luminosityBlock() const
double phi() const
azimuthal angle of momentum vector
double pt() const final
transverse momentum
const bool decision(const edm::Event &iEvent, const std::string &nameAlgoTechTrig, int &errorCode) const
std::vector< Vertex > VertexCollection
collection of Vertex objects
T * make(const Args &...args) const
make new ROOT object
double chargeIsolationCone(unsigned int trkIndex, std::vector< spr::propagatedTrackDirection > &trkDirs, double dR, int &nNearTRKs, bool debug=false)
double px() const
x coordinate of momentum vector
edm::EDGetTokenT< l1extra::L1JetParticleCollection > tok_L1extTauJet_
edm::EDGetTokenT< HBHERecHitCollection > tok_hbhe_
const std::string theTrackQuality_
EDGetTokenT< ProductType > consumes(edm::InputTag const &tag)
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
#define DEFINE_FWK_MODULE(type)
double eta() const
pseudorapidity of momentum vector
std::vector< double > pbin
static const bool useL1GtTriggerMenuLite(true)
edm::Service< TFileService > fs_
double deltaR(double eta1, double eta2, double phi1, double phi2)
bool goodTrack(const reco::Track *pTrack, math::XYZPoint leadPV, trackSelectionParameters parameters, bool debug=false)
std::vector< bool > * t_l1bits
std::vector< double > vec1
ConsumesCollector consumesCollector()
Use a ConsumesCollector to gather consumes information from helper functions.
double pt() const
track transverse momentum
int ieta() const
get the cell ieta
edm::EDGetTokenT< EcalRecHitCollection > tok_EB_
Jets made from MC generator particles.
const edm::EDGetTokenT< reco::VertexCollection > tok_recVtx_
double eCone_ecal(const CaloGeometry *geo, edm::Handle< T > &barrelhits, edm::Handle< T > &endcaphits, const GlobalPoint &hpoint1, const GlobalPoint &point1, double dR, const GlobalVector &trackMom, int &nRecHits, double ebThr=-100, double eeThr=-100, double tMin=-500, double tMax=500, bool debug=false)
~IsoTrackCalib() override
edm::EDGetTokenT< l1extra::L1JetParticleCollection > tok_L1extFwdJet_
const edm::EDGetTokenT< reco::GenJetCollection > tok_jets_
double pz() const
z coordinate of momentum vector
const_iterator end() const
double dEta(math::XYZTLorentzVector &, math::XYZTLorentzVector &)
static TrackQuality qualityByName(const std::string &name)
XYZPointD XYZPoint
point in space with cartesian internal representation
bool init(const edm::Run &iRun, const edm::EventSetup &iSetup, const std::string &processName, bool &changed)
d'tor
void add(std::string const &label, ParameterSetDescription const &psetDescription)
std::vector< PFJet > PFJetCollection
collection of PFJet objects
void getL1GtRunCache(const edm::Run &, const edm::EventSetup &, const bool, const bool)
get all the run-constant quantities for L1 trigger and cache them
reco::TrackBase::TrackQuality minQuality
void endRun(edm::Run const &, edm::EventSetup const &) override
void analyze(const edm::Event &, const edm::EventSetup &) override
const Point & position() const
position
double dR(math::XYZTLorentzVector &, math::XYZTLorentzVector &)
double dPhi(math::XYZTLorentzVector &, math::XYZTLorentzVector &)
const edm::EDGetTokenT< reco::PFJetCollection > tok_pfjets_
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
const L1GtTriggerMenu * ptrL1TriggerMenuEventSetup(int &errorCode)
return a pointer to the L1 trigger menu from event setup
T const * product() const
const std::vector< std::string > l1Names_
const_iterator begin() const
double py() const
y coordinate of momentum vector
static const bool useL1EventSetup(true)
IsoTrackCalib(const edm::ParameterSet &)