C++ JetPointer类代码示例

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;

}

extern void printPFJetExtra(const JetPointer jet) {
	cout << setw(30) << "NOD" << setw(30) << "CEF" << setw(30) << "NHF" << setw(30) << "NEF" << endl;
	cout << setw(30) << jet->NOD() << setw(30) << jet->CEF() << setw(30) << jet->NHF() << setw(30) << jet->NEF() << endl
			<< endl;

	cout << setw(30) << "CHF" << setw(30) << "NCH" << setw(30) << "" << setw(30) << "" << endl;
	cout << setw(30) << jet->CHF() << setw(30) << jet->NCH() << setw(30) << "" << setw(30) << "" << endl << endl;
}

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;
}

JetCollection Event::GetBJetCollection(const JetCollection& jets, BtagAlgorithm::value btagAlgorithm,
		BtagAlgorithm::workingPoint WP) const {
	JetCollection bjets;
	for (unsigned int index = 0; index < jets.size(); ++index) {
		const JetPointer jet = jets.at(index);
		if (jet->isBJet(btagAlgorithm, WP))
			bjets.push_back(jet);
	}

	return bjets;
}

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_);
}

void Event::selectGoodJets() {
    goodJets.clear();
    for (unsigned int index = 0; index < allJets.size(); ++index) {
        const JetPointer jet = allJets.at(index);
        if (jet->isGood()) {
            goodJets.push_back(jet);
        }
    }
    cleanGoodJets();
    for (unsigned int index = 0; index < goodJets.size(); ++index) {
        const JetPointer jet = goodJets.at(index);
        if (jet->isBJet(BtagAlgorithm::SimpleSecondaryVertexHighEffBTag))
            goodBJets.push_back(jet);
    }
}

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;
}

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;
}

void Analysis::analyse() {
	cout << "detected samples:" << endl;
	for (unsigned int sample = 0; sample < DataType::NUMBER_OF_DATA_TYPES; ++sample) {
		if (eventReader->getSeenDatatypes()[sample])
			cout << DataType::names[sample] << endl;
	}

	createHistograms();


	while (eventReader->hasNextEvent()) {

		initiateEvent();
		printNumberOfProccessedEventsEvery(Globals::printEveryXEvents);
		inspectEvents();
		const JetCollection jets(currentEvent->Jets());
		unsigned int numberOfJets(jets.size());
		unsigned int numberOfBJets(0);
		for (unsigned int index = 0; index < numberOfJets; ++index) {
			const JetPointer jet(currentEvent->Jets().at(index));
			if (jet->isBJet(BtagAlgorithm::CombinedSecondaryVertexV2, BtagAlgorithm::MEDIUM))
				++numberOfBJets;
		}
		histMan->setCurrentBJetBin(numberOfBJets);
		histMan->setCurrentJetBin(numberOfJets);

		vector<double> bjetWeights;
		if (currentEvent->isRealData()) {
			for (unsigned int index = 0; index <= numberOfBJets; ++index) {
				if (index == numberOfBJets)
					bjetWeights.push_back(1.);
				else
					bjetWeights.push_back(0);
			}
		} else
			bjetWeights = BjetWeights(jets, numberOfBJets);

		ttbar_plus_X_analyser_->analyse(currentEvent);
		if ( ( currentEvent->getDataType() == DataType::TTJets_amcatnloFXFX || currentEvent->getDataType() == DataType::TTJets_madgraphMLM ) && Globals::treePrefix_ == "" ) {
			pseudoTopAnalyser_->analyse(currentEvent);
			unfoldingRecoAnalyser_->analyse(currentEvent);
			partonAnalyser_->analyse(currentEvent);
			// likelihoodInputAnalyser_->analyse(currentEvent);
		}
		treeMan->FillTrees();
	}
}

const JetCollection TopPairMuPlusJetsReferenceSelection2011::cleanedJets(const EventPtr event) const {
	const JetCollection jets(event->Jets());
	JetCollection cleanedJets;

	if (!hasExactlyOneIsolatedLepton(event)) //if no signal lepton is found, can't clean jets, return them all!
		return jets;

	const LeptonPointer lepton(signalLepton(event));

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

	return cleanedJets;
}

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);
}

int HitFitAnalyser::positionOfLastTTBarJet(const JetCollection jets) {
	// Loop over jets and find position of last jet that comes from ttbar decay
	bool foundHadB = false;
	bool foundLepB = false;
	bool foundQ = false;
	bool foundQBar = false;
	for ( unsigned int jetIndex=0; jetIndex < jets.size(); ++jetIndex ) {
		JetPointer jet = jets[jetIndex];

		if ( jet->ttbar_decay_parton() == 3 ) foundQ = true;
		else if ( jet->ttbar_decay_parton() == 4 ) foundQBar = true;
		else if ( jet->ttbar_decay_parton() == 5 ) foundLepB = true;
		else if ( jet->ttbar_decay_parton() == 6 ) foundHadB = true;

		if ( foundQ && foundQBar && foundLepB && foundHadB ) return jetIndex;
	}
	return -1;

}

void Analysis::analyse() {
	createHistograms();
	cout << "detected samples:" << endl;
	for (unsigned int sample = 0; sample < DataType::NUMBER_OF_DATA_TYPES; ++sample) {
		if (eventReader->getSeenDatatypes()[sample])
			cout << DataType::names[sample] << endl;
	}
	while (eventReader->hasNextEvent()) {
		initiateEvent();
		printNumberOfProccessedEventsEvery(Globals::printEveryXEvents);
		inspectEvents();
		const JetCollection jets(currentEvent->Jets());
		unsigned int numberOfJets(jets.size());
		unsigned int numberOfBJets(0);
		for (unsigned int index = 0; index < numberOfJets; ++index) {
			const JetPointer jet(currentEvent->Jets().at(index));
			if (jet->isBJet(BtagAlgorithm::CombinedSecondaryVertex, BtagAlgorithm::MEDIUM))
				++numberOfBJets;
		}
		histMan->setCurrentBJetBin(numberOfBJets);
		histMan->setCurrentJetBin(numberOfJets);

		vector<double> bjetWeights;
		if (currentEvent->isRealData()) {
			for (unsigned int index = 0; index <= numberOfBJets; ++index) {
				if (index == numberOfBJets)
					bjetWeights.push_back(1.);
				else
					bjetWeights.push_back(0);
			}
		} else
			bjetWeights = BjetWeights(jets, numberOfBJets);

		eventcountAnalyser->analyse(currentEvent);
//		mttbarAnalyser->analyse(currentEvent);
		ttbar_plus_X_analyser_->analyse(currentEvent);
		diffVariablesAnalyser->analyse(currentEvent);
//		binningAnalyser->analyse(currentEvent);
	}
}

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();
}

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;
}