C++ JetCollection::size方法代码示例

//this can be done shorter if you introduce a function which passes the btagWeight
std::vector<double> BjetWeights(const JetCollection& jets, unsigned int numberOfBtags, int btagSystematicFactor,
		int lightJetSystematicFactor) {
	boost::scoped_ptr < BTagWeight > btagWeight(new BTagWeight(btagSystematicFactor, lightJetSystematicFactor));
	//get b-jets
	const JetCollection bjets(btagWeight->getBJets(jets));
	//get c-jets
	const JetCollection cjets(btagWeight->getCJets(jets));
	//get udsg jets
	const JetCollection udsgjets(btagWeight->getUDSGJets(jets));

	//get mean scale factors
	double SF_b = btagWeight->getAverageBScaleFactor(bjets);
	double SF_c = btagWeight->getAverageCScaleFactor(cjets);
	double SF_udsg = btagWeight->getAverageUDSGScaleFactor(udsgjets);
	//get mean efficiencies
	double mean_bJetEfficiency = btagWeight->getAverageBEfficiency();
	double mean_cJetEfficiency = btagWeight->getAverageCEfficiency();
	double mean_udsgJetEfficiency = btagWeight->getAverageUDSGEfficiency(udsgjets);

	std::vector<double> event_weights;
	for (unsigned int nTag = 0; nTag <= numberOfBtags; ++nTag) { // >= 4 is our last b-tag bin!
		btagWeight->setNumberOfBtags(nTag, 20);
		double event_weight = btagWeight->weight(bjets.size(), cjets.size(), udsgjets.size(), mean_bJetEfficiency,
				mean_cJetEfficiency, mean_udsgJetEfficiency, SF_b, SF_c, SF_udsg, numberOfBtags);
		event_weights.push_back(event_weight);
	}
	//all weights are inclusive. To get the weight for exclusive N b-tags ones has to subtract:
	for (unsigned int nTag = 0; nTag < numberOfBtags; ++nTag) {
		// w(N b-tags) = w(>= N) - w(>= N+1)
		event_weights.at(nTag) = event_weights.at(nTag) - event_weights.at(nTag + 1);
		//last weight, >= numberOfBjets jets, stays inclusive
	}
	return event_weights;
}

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 TTbarDiLeptonAnalyser::eESignalAnalysis(const EventPtr event) {

	if (topEERefSelection_->passesSelectionUpToStep(event, TTbarEEReferenceSelection::MetCut)) {

			const JetCollection jets(topEERefSelection_->cleanedJets(event));
			const JetCollection bJets(topEERefSelection_->cleanedBJets(event));
			unsigned int numberOfBjets(bJets.size());
			vector<double> bjetWeights;
			if (event->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);
			histMan_->setCurrentJetBin(jets.size());
			histMan_->setCurrentBJetBin(numberOfBjets);

			//this is for HT in the MET analyser but needs changing
			const LeptonPointer signalLepton = topEERefSelection_->signalLepton(event);
			const ElectronPointer signalElectron(boost::static_pointer_cast<Electron>(signalLepton));
			const PhotonCollection photons = topEERefSelection_->signalPhotons(event);

			//get dilepton collection
			const ElectronCollection electrons = topEERefSelection_->signalElectrons(event);
			const MuonCollection muons = topEERefSelection_->signalMuons(event);

			for (unsigned int weightIndex = 0; weightIndex < bjetWeights.size(); ++weightIndex) {

				double bjetWeight = bjetWeights.at(weightIndex);
				histMan_->setCurrentBJetBin(weightIndex);
				histMan_->setCurrentHistogramFolder(histogramFolder_ + "/EE/Ref selection");
				//MET
				metAnalyserEERefSelection_->setScale(bjetWeight);
				metAnalyserEERefSelection_->analyse(event, signalLepton);

				//jets
				jetAnalyserEERefSelection_->setScale(bjetWeight);
				jetAnalyserEERefSelection_->analyse(event);

				//DiElectron
				diElectronAnalyserEERefSelection_->setScale(bjetWeight);
			    diElectronAnalyserEERefSelection_->analyse(event, electrons);

			    //photons
			    photonAnalyserEERefSelection_->setScale(bjetWeight);
			    photonAnalyserEERefSelection_->analyse(event, photons, jets, electrons, muons);

			}
	}
}

void MET::adjForUnc(const JetCollection &jets)
{
	if (Jet::correctDirection == JetCorrDirection::NONE ||
		jets.size() <= 0)
		return;
	FourVector new4vec = getFourVector();
	for (unsigned int index = 0; index < jets.size(); ++index)
		new4vec += jets.at(index)->DiffVec;
	new4vec.SetPz(0.0);  // Make sure no longitudinal component
	new4vec.SetE(new4vec.Pt()); // Make sure it stays massless
	setFourVector(new4vec);
}

extern void printJets(const JetCollection jets) {
	for (unsigned short index = 0; index < jets.size(); ++index) {
		const JetPointer jet = jets.at(index);
		cout << "Jet " << index + 1 << endl;
		printJet(jet);
	}
}

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;

}

JetCollection BTagWeight::getBJets(const JetCollection jets) const {
	JetCollection bjets;
	for (unsigned int index = 0; index < jets.size(); ++index) {
		if (abs(jets.at(index)->partonFlavour()) == 5) //b-quark
			bjets.push_back(jets.at(index));
	}
	return bjets;
}

JetCollection BTagWeight::getCJets(const JetCollection& jets) const {
	JetCollection cjets;
	for (unsigned int index = 0; index < jets.size(); ++index) {
		if (abs(jets.at(index).partonFlavour()) == 4) //c-quark
			cjets.push_back(jets.at(index));
	}
	return cjets;
}

JetCollection BTagWeight::getUDSGJets(const JetCollection& jets) const {
	JetCollection udsgjets;
	for (unsigned int index = 0; index < jets.size(); ++index) {
		if (abs(jets.at(index).partonFlavour()) != 4 && abs(jets.at(index).partonFlavour()) != 5) //not a c- or b-quark
			udsgjets.push_back(jets.at(index));
	}
	return udsgjets;
}

void TTbar_plus_X_analyser::fillCommonTreesNoBSelection(const EventPtr event,  const unsigned int selectionCriteria, std::string folder ) {
	SelectionCriteria::selection selection = SelectionCriteria::selection(selectionCriteria);

	// Jets
	const JetCollection jets(event->CleanedJets());
	// B Jets
	const JetCollection bJets(event->CleanedBJets());
	// Lepton
	const LeptonPointer signalLepton = event->getSignalLepton( selection );
	// MET
	const METPointer MET_original(event->MET((METAlgorithm::value) 0));


	treeMan_->setCurrentFolder(folder);
	treeMan_->Fill("EventWeight", event->weight());
	treeMan_->Fill("PUWeight", event->PileUpWeight());
	treeMan_->Fill("PUWeight_up", event->PileUpWeight(1));
	treeMan_->Fill("PUWeight_down", event->PileUpWeight(-1));

	treeMan_->Fill("M3",Event::M3(jets));

	if ( Event::NJets(bJets) > 0 ) {
		treeMan_->Fill("M_bl",Event::M_bl(bJets, signalLepton));
		treeMan_->Fill("angle_bl",Event::angle_bl(bJets, signalLepton));
	}

	treeMan_->Fill("HT",Event::HT(jets));
	treeMan_->Fill("MET",MET_original->et());
	treeMan_->Fill("MET_phi",MET_original->phi());
	treeMan_->Fill("ST",Event::ST(jets, signalLepton, MET_original));
	treeMan_->Fill("WPT",Event::WPT(signalLepton, MET_original));
	treeMan_->Fill("MT",Event::MT(signalLepton, MET_original));

	treeMan_->Fill("NJets",Event::NJets(jets));
	treeMan_->Fill("NBJets",Event::NJets(bJets));
	treeMan_->Fill("NVertices",	event->Vertices().size());

	treeMan_->Fill("BJetWeight",event->BJetWeight());
	treeMan_->Fill("BJetEfficiencyCorrectionWeight",event->BJetEfficiencyCorrectionWeight());
	treeMan_->Fill("BJetAlternativeWeight",event->BJetAlternativeWeight());

	treeMan_->Fill("BJetUpWeight",event->BJetUpWeight());
	treeMan_->Fill("BJetDownWeight",event->BJetDownWeight());
	treeMan_->Fill("LightJetUpWeight",event->LightJetUpWeight());
	treeMan_->Fill("LightJetDownWeight",event->LightJetDownWeight());

	if (selection == SelectionCriteria::selection(SelectionCriteria::ElectronPlusJetsReference)){
		treeMan_->Fill("lepton_isolation", signalLepton->PFRelIso03DeltaBeta());
	}
	else if (selection == SelectionCriteria::selection(SelectionCriteria::MuonPlusJetsReference)){
		treeMan_->Fill("lepton_isolation", signalLepton->PFRelIso04DeltaBeta());
	}

	for (unsigned int index = 0; index < jets.size(); ++index) {
		treeMan_->Fill("jet_csv", jets.at(index)->getBTagDiscriminator(BtagAlgorithm::CombinedSecondaryVertexV2) );
	}	
	fillLeptonEfficiencyCorrectionBranches( event, selectionCriteria, signalLepton );	
}

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 Event::HT(const JetCollection jets) {
	double ht(0);
	//Take ALL the jets!
	for (unsigned int index = 0; index < jets.size(); ++index) {
		if(jets.at(index)->pt() > 20)
			ht += jets.at(index)->pt();
	}
	return ht;
}

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

double TtbarHypothesis::M3(const JetCollection jets) {
	double m3(0), max_pt(0);
	if (jets.size() >= 3) {
		for (unsigned int index1 = 0; index1 < jets.size() - 2; ++index1) {
			for (unsigned int index2 = index1 + 1; index2 < jets.size() - 1; ++index2) {
				for (unsigned int index3 = index2 + 1; index3 < jets.size(); ++index3) {
					FourVector m3Vector(
							jets.at(index1)->getFourVector() + jets.at(index2)->getFourVector()
									+ jets.at(index3)->getFourVector());
					double currentPt = m3Vector.Pt();
					if (currentPt > max_pt) {
						max_pt = currentPt;
						m3 = m3Vector.M();
					}
				}
			}
		}
	}

	return m3;
}

unsigned short Electron::getClosestJetIndex(const JetCollection& jets) const {
	unsigned short idOfClosest = 999;
	double closestDR = 999.;
	for (unsigned short index = 0; index < jets.size(); ++index) {
		double DR = deltaR(jets.at(index));
		if (DR < closestDR) {
			closestDR = DR;
			idOfClosest = index;
		}
	}
	return idOfClosest;
}