C++ JetPointer::pt方法代码示例

void JetAnalyser::analyse(const EventPtr event) {

    histMan_->setCurrentHistogramFolder(histogramFolder_);
    weight_ = event->weight() * prescale_ * scale_;
    const JetCollection jets = event->Jets();
    unsigned int numberOfBJets(0);

    for (unsigned int index = 0; index < jets.size(); ++index) {

        const JetPointer jet(jets.at(index));
        histMan_->H1D_BJetBinned("all_jet_pT")->Fill(jet->pt(), weight_);
        histMan_->H1D_BJetBinned("all_jet_phi")->Fill(jet->phi(), weight_);
        histMan_->H1D_BJetBinned("all_jet_eta")->Fill(jet->eta(), weight_);
        if (jet->isBJet(BtagAlgorithm::CombinedSecondaryVertex, BtagAlgorithm::MEDIUM))
            ++numberOfBJets;

        if (index < 7) {
            stringstream temp;
            temp << "jet" << (index + 1);
            string nthJet = temp.str();
            histMan_->H1D_BJetBinned(nthJet + "_pT")->Fill(jet->pt(), weight_);
            histMan_->H1D_BJetBinned(nthJet + "_phi")->Fill(jet->phi(), weight_);
            histMan_->H1D_BJetBinned(nthJet + "_eta")->Fill(jet->eta(), weight_);
        }

    }

    histMan_->H1D_BJetBinned("N_Jets")->Fill(jets.size(), weight_);
    histMan_->H1D("N_BJets")->Fill(numberOfBJets, weight_);
}

bool TopPairEplusJetsRefAsymJetsSelection::passesAsymmetricJetCuts(const EventPtr event) const {
	const JetCollection goodElectronCleanedJets = cleanedJets(event);
	if (goodElectronCleanedJets.size() < 3) // good jets have a cut of 30 GeV!
		return false;
	JetPointer leadingJet = goodElectronCleanedJets.front();
	JetPointer secondLeadingJet = goodElectronCleanedJets.at(1);
	JetPointer thirdLeadingJet = goodElectronCleanedJets.at(2);
	return leadingJet->pt() > 70 && secondLeadingJet->pt() > 50 && thirdLeadingJet->pt() > 50;
}

double BTagWeight::getUDSGScaleFactor(const JetPointer jet) const {
	double pt = jet->pt();
	double SF_udsg_mean(0), SF_udsg_min(0), SF_udsg_max(0);

	if (pt < 20) {
		return 0;
	} else if (pt > 670) {
		SF_udsg_mean = getMeanUDSGScaleFactor(670.);
		SF_udsg_min = getMinUDSGScaleFactor(670);
		SF_udsg_max = getMaxUDSGScaleFactor(670);
		//use twice the uncertainty
		SF_udsg_min -= (SF_udsg_mean - SF_udsg_min);
		SF_udsg_max += (SF_udsg_max - SF_udsg_mean);
	} else {
		SF_udsg_mean = getMeanUDSGScaleFactor(pt);
		SF_udsg_min = getMinUDSGScaleFactor(pt);
		SF_udsg_max = getMaxUDSGScaleFactor(pt);
	}
	if (Globals::LightJetSystematic == -1)
		return SF_udsg_min;
	else if (Globals::LightJetSystematic == 1)
		return SF_udsg_max;

	return SF_udsg_mean;
}

double BTagWeight::getBScaleFactor(const JetPointer jet, double uncertaintyFactor) const {
	const boost::array<double, 14> SFb_error = { { 0.0295675, 0.0295095, 0.0210867, 0.0219349, 0.0227033, 0.0204062,
			0.0185857, 0.0256242, 0.0383341, 0.0409675, 0.0420284, 0.0541299, 0.0578761, 0.0655432 } };

	const boost::array<double, 14> ptbins = { { 30, 40, 50, 60, 70, 80, 100, 120, 160, 210, 260, 320, 400, 500 } };

	double SFb(0);
	double sf_error(0);
	//these numbers are for CSVM only
	double pt = jet->pt();
	if (pt < 30) {
		SFb = 0.6981 * (1. + 0.414063 * 30) / (1. + 0.300155 * 30);
		sf_error = 0.12;
	} else if (pt > 670) {
		SFb = 0.6981 * (1. + 0.414063 * 670) / (1. + 0.300155 * 670);
		//use twice the uncertainty
		sf_error = 2 * SFb_error[SFb_error.size() - 1];
	} else {
		SFb = 0.6981 * (1. + 0.414063 * pt) / (1. + 0.300155 * pt);
		unsigned int ptbin(0);
		for (unsigned int bin = 0; bin < ptbins.size() + 1; ++bin) {
			double upperCut = bin + 1 < ptbins.size() ? ptbins.at(bin + 1) : 670.;
			double lowerCut = ptbins.at(bin);

			if (pt > lowerCut && pt <= upperCut) {
				ptbin = bin;
				break;
			}
		}
		sf_error = SFb_error.at(ptbin);
	}
	SFb += sf_error * Globals::BJetSystematic * uncertaintyFactor;
	return SFb;
}

bool TopPairEMuReferenceSelection::isGoodJet(const JetPointer jet) const {

	/**
	 * This function tests the jet ID and eta (and pt) range for jet
	 * The cut of 20 GeV is actually obsolete since we only store jets above that threshold.
	 * However, the jet id is only valid for jets above it.
	 */

	bool passesPtAndEta = jet->pt() > 20 && fabs(jet->eta()) < 2.5;
	bool passesJetID(false);
	JetAlgorithm::value algo = jet->getUsedAlgorithm();
	if (algo == JetAlgorithm::CA08PF || algo == JetAlgorithm::PF2PAT) { //PFJet
		bool passNOD = jet->NOD() > 1;
		bool passNHF = jet->NHF() < 0.99;
		bool passNEF = jet->NEF() < 0.99;
		bool passCHF = true;
		bool passNCH = true;
		bool passCEF = true;
		if (fabs(jet->eta()) < 2.4) {
			passCEF = jet->CEF() < 0.99;
			passCHF = jet->CHF() > 0;
			passNCH = jet->NCH() > 0;
		}
		passesJetID = passNOD && passCEF && passNHF && passNEF && passCHF && passNCH;

	} else { //assume CaloJet
		bool passesEMF = jet->emf() > 0.01;
		bool passesN90Hits = jet->n90Hits() > 1;
		bool passesFHPD = jet->fHPD() < 0.98;
		passesJetID = passesEMF && passesN90Hits && passesFHPD;
	}
	return passesPtAndEta && passesJetID;

}

bool TopPairMuPlusJetsReferenceSelection2011::isGoodJet(const JetPointer jet) const {
	bool passesPtAndEta = jet->pt() > 20 && fabs(jet->eta()) < 2.4;
	bool passesJetID(false);
	JetAlgorithm::value algo = jet->getUsedAlgorithm();
	if (algo == JetAlgorithm::CA08PF || algo == JetAlgorithm::PF2PAT) { //PFJet
		bool passNOD = jet->NOD() > 1;
		bool passNHF = jet->NHF() < 0.99;
		bool passNEF = jet->NEF() < 0.99;
		bool passCHF = jet->CEF() < 0.99;
		bool passNCH = true;
		bool passCEF = true;
		if (fabs(jet->eta()) < 2.4) {
			passCHF = jet->CHF() > 0;
			passNCH = jet->NCH() > 0;
		}
		passesJetID = passNOD && passCEF && passNHF && passNEF && passCHF && passNCH;
		;
	} else { //assume CaloJet
		bool passesEMF = jet->emf() > 0.01;
		bool passesN90Hits = jet->n90Hits() > 1;
		bool passesFHPD = jet->fHPD() < 0.98;
		passesJetID = passesEMF && passesN90Hits && passesFHPD;
	}
	return passesPtAndEta && passesJetID;

}

bool PseudoTopAnalyser::passesEventSelection( const MCParticlePointer pseudoLepton, const ParticlePointer pseudoNeutrino, const JetCollection pseudoJets, const MCParticleCollection pseudoBs, const ParticleCollection allPseudoLeptons, const ParticlePointer pseudoMET ) {

	// Event selection taken from here : https://twiki.cern.ch/twiki/bin/view/LHCPhysics/ParticleLevelTopDefinitions
	unsigned int numberGoodLeptons = 0;
	unsigned int numberVetoLeptons = 0;
	ParticlePointer leadingLepton;
	for ( unsigned int leptonIndex = 0; leptonIndex < allPseudoLeptons.size(); ++ leptonIndex ) {
		const ParticlePointer lepton = allPseudoLeptons.at(leptonIndex);

		// Check if this is a good signal type lepton
		if ( lepton->pt() > minLeptonPt_ && fabs(lepton->eta()) < maxLeptonAbsEta_ ) {
			++numberGoodLeptons;
			if ( leadingLepton == 0 ) leadingLepton = lepton;
		}
		
		// Check if this is a veto lepton
		if ( lepton->pt() > minVetoLeptonPt_ && fabs(lepton->eta()) < maxVetoLeptonAbsEta_ ) {
			++numberVetoLeptons;
		}
	}

	// Neutrino pt sum
	bool passesNeutrinoSumPt = false;
	if ( pseudoMET != 0 ) {
		if ( pseudoMET->pt() > minNeutrinoSumPt_ ) passesNeutrinoSumPt = true;
	}

	// W MT
	bool passesWMT = false;
	if ( leadingLepton != 0 && pseudoMET != 0 ) {
		double genMT = sqrt( 2 * leadingLepton->pt() * pseudoMET->pt() * ( 1 - cos(leadingLepton->phi() - pseudoMET->phi() ) ) );
		if (genMT > minWMt_) passesWMT = true;
	}


	// Jets
	unsigned int numberGoodJets = 0;
	unsigned int numberGoodBJets = 0;

	for ( unsigned int jetIndex = 0; jetIndex < pseudoJets.size(); ++ jetIndex ) {
		const JetPointer jet = pseudoJets.at(jetIndex);

		// Check if this is a good jet
		if ( jet->pt() > minJetPt_ && fabs(jet->eta()) < maxJetAbsEta_ ) {
			++numberGoodJets;

			// Check if this is also a good b jet
			if ( fabs( jet->partonFlavour() ) == 5 ) {
				++numberGoodBJets;
			}
		}
	}
	
	if ( numberGoodLeptons == 1 && numberVetoLeptons <= 1 && passesNeutrinoSumPt && passesWMT && numberGoodJets >= minNJets_ && numberGoodBJets >= minNBJets_ ) {
		return true;
	}
	else return false;

}

double BTagWeight::getGEfficiency(const JetPointer jet) const {

	const double jetPt = jet->pt();
	const double jetEta = jet->eta();
	// std::cout << "Jet Pt : " << jetPt << ", Jet Eta : " << jetEta << std::endl;
	int binNumber = Globals::gluonJet->FindBin( jetPt , jetEta );
	float BTagEff = Globals::gluonJet->GetBinContent( binNumber );	
	// std::cout << "Gluon Jet, B Tag Efficiency : " << BTagEff << std::endl;
	return BTagEff;
}

bool HLTriggerQCDAnalyser::passesTriggerAnalysisSelection(const EventPtr event) const {
	const ElectronCollection electrons(event->Electrons());
	const JetCollection jets(event->Jets());
	if (electrons.size() == 0 || jets.size() < 3)
		return false;

	unsigned int nElectrons(0);
	for (unsigned int index = 0; index < electrons.size(); ++index) {
		const ElectronPointer electron(electrons.at(index));
		if (fabs(electron->eta()) < 2.5 && electron->pt() > 20)
			++nElectrons;
		//if more than 2 electrons passing the selection of > 20GeV, reject event
	}
	const ElectronPointer mostEnergeticElectron(electrons.front());
	//clean jets against electron
	JetCollection cleanedJets;

	for (unsigned int index = 0; index < jets.size(); ++index) {
		const JetPointer jet(jets.at(index));
		if (!jet->isWithinDeltaR(0.3, mostEnergeticElectron))
			cleanedJets.push_back(jet);
	}

	unsigned int nCleanedJetsAbove30GeV(0), nCleanedJetsAbove45GeV(0);
	for (unsigned int index = 0; index < cleanedJets.size(); ++index) {
		const JetPointer jet(cleanedJets.at(index));
		if (jet->pt() > 45.)
			++nCleanedJetsAbove45GeV;
		if (jet->pt() > 30.)
			++nCleanedJetsAbove30GeV;
	}

	return nElectrons == 1
			&& (nCleanedJetsAbove45GeV >= 3
					|| (nCleanedJetsAbove45GeV >= 2 && nCleanedJetsAbove30GeV >= 3 && event->runnumber() >= 194270));
}

void BTagEff::analyse(const EventPtr event) {

	histMan_->setCurrentHistogramFolder(histogramFolder_);
	treeMan_->setCurrentFolder(histogramFolder_);
	int NJets = 0;
	const JetCollection allJets = event->Jets();

	for (unsigned int jetIndex = 0; jetIndex < allJets.size(); ++jetIndex) {
		const JetPointer jet(allJets.at(jetIndex));

		bool isLoose = false;
		bool isMedium = false;
		bool isTight = false;

		double jetPt = jet->pt();
		double jetEta = jet->eta();

		if (jetPt < 25 || abs(jetEta) > 2.4) continue;
		// double jetCSV = jet->getBTagDiscriminator(BtagAlgorithm::CombinedSecondaryVertexV2, BtagAlgorithm::MEDIUM);
		double jetCSV = jet->getBTagDiscriminator(BAT::BtagAlgorithm::value::CombinedSecondaryVertexV2);

		// https://twiki.cern.ch/twiki/bin/viewauth/CMS/BtagRecommendation74X50ns
		if (jetCSV > 0.605) isLoose = true;
		if (jetCSV > 0.890) isMedium = true;
		if (jetCSV > 0.970) isTight = true;

		unsigned int partonFlavour = abs(jet->partonFlavour());
		// const bool isBTagged = jet->isBJet(BtagAlgorithm::CombinedSecondaryVertexV2, BtagAlgorithm::MEDIUM);
		// cout << jet->isBJet(BtagAlgorithm::CombinedSecondaryVertexV2, BtagAlgorithm::MEDIUM) << endl;
		
		treeMan_->Fill("pt", jetPt);
		treeMan_->Fill("eta", jetEta);
		treeMan_->Fill("CSV", jetCSV);
		treeMan_->Fill("partonFlavour", partonFlavour);
		treeMan_->Fill("isLoose", isLoose);
		treeMan_->Fill("isMedium", isMedium);
		treeMan_->Fill("isTight", isTight);
		++NJets;
	}


	treeMan_->Fill("NJets", NJets);
}

double BTagWeight::getAverageUDSGEfficiency(const JetCollection jets) const {
	std::vector<double> efficiencies;

	for (unsigned int index = 0; index < jets.size(); ++index) {
		const JetPointer jet(jets.at(index));
		double efficiency(0);
		//these numbers are for CSVM only
		double pt = jet->pt();
		if (pt < 20) {
			continue;
		} else if (pt > 670) {
			efficiency = getMeanUDSGEfficiency(670.);
		} else {
			efficiency = getMeanUDSGEfficiency(pt);
		}
		efficiencies.push_back(efficiency);
	}
	double sumOfEfficiencies = std::accumulate(efficiencies.begin(), efficiencies.end(), 0.0);
	if (efficiencies.size() == 0)
		return 1.;
	else
		return sumOfEfficiencies / efficiencies.size();
}

void BTagEff::analyse(const EventPtr event) {

	treeMan_->setCurrentFolder(histogramFolder_);
	int NJets = 0;
	int NBJets = 0; // How many medium b jets
	const JetCollection cleanedJets = event->CleanedJets();

	int selectionCriteria = -1;
	if ( event->PassesElectronTriggerAndSelectionNoB() ) selectionCriteria = SelectionCriteria::ElectronPlusJetsReference;
	else if ( event->PassesMuonTriggerAndSelectionNoB() ) selectionCriteria = SelectionCriteria::MuonPlusJetsReference;
	const LeptonPointer signalLepton = event->getSignalLepton( selectionCriteria );

	// unsigned int nParton = 0;
	for (unsigned int jetIndex = 0; jetIndex < cleanedJets.size(); ++jetIndex) {
		const JetPointer jet(cleanedJets.at(jetIndex));

		bool isLoose = false;
		bool isMedium = false;
		bool isTight = false;

		double jetPt = jet->pt();
		double jetEta = jet->eta();

		if (jetPt < 25 || fabs(jetEta) > 2.4) continue;

		double jetCSV = jet->getBTagDiscriminator(BAT::BtagAlgorithm::value::CombinedSecondaryVertexV2);

		// https://twiki.cern.ch/twiki/bin/viewauth/CMS/BtagRecommendation76X
		if (jetCSV > 0.460) {
			isLoose = true;
		}
		if (jetCSV > 0.800) {
			isMedium = true;
			++NBJets;
		}
		if (jetCSV > 0.935) {
			isTight = true;
		}

		unsigned int partonFlavour = abs(jet->partonFlavour());
		unsigned int hadronFlavour = abs(jet->hadronFlavour());

		treeMan_->Fill("pt", jetPt);
		treeMan_->Fill("eta", jetEta);
		treeMan_->Fill("CSV", jetCSV);
		treeMan_->Fill("partonFlavour", partonFlavour);
		treeMan_->Fill("hadronFlavour", hadronFlavour);
		treeMan_->Fill("isLoose", isLoose);
		treeMan_->Fill("isMedium", isMedium);
		treeMan_->Fill("isTight", isTight);
		++NJets;
	}

	treeMan_->Fill("NJets", NJets);
	treeMan_->Fill("NBJets", NBJets);
	treeMan_->Fill("EventWeight", event->weight());
	treeMan_->Fill("PUWeight", event->PileUpWeight());

	if ( selectionCriteria == SelectionCriteria::ElectronPlusJetsReference ) {
		double electronEfficiencyCorrection = 1;
		if ( !event->isRealData() ) {
			const ElectronPointer signalElectron(boost::static_pointer_cast<Electron>(signalLepton));
			electronEfficiencyCorrection = signalElectron->getEfficiencyCorrection( 0 );
		}
		treeMan_->Fill("ElectronEfficiencyCorrection",electronEfficiencyCorrection);
	}
	else if ( selectionCriteria == SelectionCriteria::MuonPlusJetsReference ) {
		double muonEfficiencyCorrection = 1;
		if ( !event->isRealData() ) {
			const MuonPointer signalMuon(boost::static_pointer_cast<Muon>(signalLepton));
			muonEfficiencyCorrection = signalMuon->getEfficiencyCorrection( 0 );
		}
		treeMan_->Fill("MuonEfficiencyCorrection",muonEfficiencyCorrection);
	}
}

double BTagWeight::weight(const JetCollection jets, const int systematic) const {
	float bTaggedMCJet = 1.0;
	float nonBTaggedMCJet = 1.0;
	float bTaggedDataJet = 1.0;
	float nonBTaggedDataJet = 1.0;

	for (unsigned int index = 0; index < jets.size(); ++index) {
		// Info on this jet
		const JetPointer jet(jets.at(index));

		double jetPt = jet->pt();
		if ( jetPt < 25 ) continue;

		// If the pt of the jet is outside the pt range of the SFs,
		// use the pt at the upper/lower edge and double the uncertainty.
		bool ptOutOfRange = false;
		if ( jetPt <= 30 || jet->pt() >= 670 ) {
			ptOutOfRange = true;
		}

		const unsigned int partonFlavour = abs( jet->partonFlavour() );
		const bool isBTagged = jet->isBJet();

		// Get scale factor for this jet
		const double sf = jet->getBTagSF( 0 );
		double sf_up = jet->getBTagSF( 1 );
		double sf_down = jet->getBTagSF( -1 );

		if ( ptOutOfRange ) {
			sf_up = sf + 2 * ( sf_up - sf );
			sf_down = sf - 2 * ( sf - sf_down );

			if ( sf_up < 0 ) sf_up = 0;
			if ( sf_down < 0 ) sf_down = 0;
		}

		// Get efficiency for this jet
		const double eff = getEfficiency( partonFlavour, jet );
		double sfToUse = sf;
		if ( systematic == 1 ) {
			sfToUse = sf_up;
		}
		else if ( systematic == -1 ) {
			sfToUse = sf_down;
		}

		if ( isBTagged ) {
			bTaggedMCJet *= eff;
			if ( eff*sfToUse > 1 ) {
				bTaggedDataJet *= 1;
			}
			else if ( eff*sfToUse < 0 ) {
				bTaggedDataJet *= 0;
			}
			else {
				bTaggedDataJet *= eff*sfToUse;
			}
		}
		else {
			nonBTaggedMCJet *= ( 1 - eff );

			if ( eff*sfToUse > 1 ) {
				nonBTaggedDataJet *= 0;
			}
			else if ( eff*sfToUse < 0 ) {
				nonBTaggedDataJet *= 1;
			}
			else {
				bTaggedDataJet *= ( 1 - eff*sfToUse );
			}
		}
		// if ( nonBTaggedMCJet < 0 || nonBTaggedDataJet < 0 ) {
		// 	cout << nonBTaggedMCJet << " " << nonBTaggedDataJet << endl;
		// 	cout << eff << " " << sfToUse << endl;
		// }
	}

	double bTagWeight = (nonBTaggedDataJet * bTaggedDataJet) / (nonBTaggedMCJet * bTaggedMCJet);

	return bTagWeight;
}

void JetAnalyser::analyse(const EventPtr event) {

	histMan_->setCurrentHistogramFolder(histogramFolder_);
	treeMan_->setCurrentFolder(histogramFolder_);
	
	weight_ = event->weight() * prescale_ * scale_;
	// const JetCollection jets = event->Jets();
	const JetCollection jets(event->CleanedJets());
	const JetCollection bjets(event->CleanedBJets());

	unsigned int numberOfBJets = event->NJets(bjets);
	unsigned int numberOfJets = event->NJets(jets);

	treeMan_->Fill("NJets", numberOfJets);
	treeMan_->Fill("NBJets", numberOfBJets);

	for (unsigned int index = 0; index < jets.size(); ++index) {
		const JetPointer jet(jets.at(index));

		if (jet->pt() < 25 ) continue;

		if (index < 5) {
			stringstream temp;
			temp << "jet" << (index + 1);
			string nthJet = temp.str();

			treeMan_->Fill( nthJet + "_pt", jet->pt());
			treeMan_->Fill( nthJet + "_eta", jet->eta());
			treeMan_->Fill( nthJet + "_phi", jet->phi());
		}

		treeMan_->Fill("pt", jet->pt());
		treeMan_->Fill("eta", jet->eta());
		treeMan_->Fill("phi", jet->phi());
	}

	for (unsigned int index = 0; index < bjets.size(); ++index) {
		const JetPointer bJet(bjets.at(index));

		if ( bJet->pt() < 25 ) continue;

		treeMan_->Fill("bjet_pt", bJet->pt());
		treeMan_->Fill("bjet_eta", bJet->eta());
		treeMan_->Fill("bjet_phi", bJet->phi());
	}

	// if ( numberOfJets < 4 ) {
	// 	cout << "Fewer than 4 good cleaned jets with pt > 25 GeV : " << numberOfJets << endl;
	// 	cout << numberOfJets << " " << numberOfBJets << endl;
	// }


	treeMan_->Fill("EventWeight", weight_ );

	if ( event->PassesMuonTriggerAndSelection() ) {
		const LeptonPointer signalLepton = event->getSignalLepton( SelectionCriteria::MuonPlusJetsReference );
		const MuonPointer signalMuon(boost::static_pointer_cast<Muon>(signalLepton));
		double efficiencyCorrection = event->isRealData() ? 1. : signalMuon->getEfficiencyCorrection( 0 );
		treeMan_->Fill("MuonEfficiencyCorrection", efficiencyCorrection);
	}
	else if ( event->PassesElectronTriggerAndSelection() ) {
		const LeptonPointer signalLepton = event->getSignalLepton( SelectionCriteria::ElectronPlusJetsReference );
		const ElectronPointer signalElectron(boost::static_pointer_cast<Electron>(signalLepton));
		double efficiencyCorrection = event->isRealData() ? 1. : signalElectron->getEfficiencyCorrection( 0 );	
		treeMan_->Fill("ElectronEfficiencyCorrection", efficiencyCorrection);
	}
}