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

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

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

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

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

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

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

}

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

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

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

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

double BTagWeight::getAverageUDSGScaleFactor(const JetCollection& jets) const {
	std::vector<double> scaleFactors;

	for (unsigned int index = 0; index < jets.size(); ++index) {
		const Jet jet(jets.at(index));
		scaleFactors.push_back(getUDSGScaleFactor(jet));
	}
	double sumOfScaleFactors = std::accumulate(scaleFactors.begin(), scaleFactors.end(), 0.0);
	if (scaleFactors.size() == 0)
		return 1.;
	else
		return sumOfScaleFactors / scaleFactors.size();
}

double BTagWeight::weight ( const JetCollection& jets, int BTagSystematic, std::string MCSampleTag) const {
	float bTaggedMCJet = 1.0;
	float nonBTaggedMCJet = 1.0;
	float bTaggedDataJet = 1.0;
	float nonBTaggedDataJet = 1.0;

	float err1 = 0.0;
	float err2 = 0.0;
	float err3 = 0.0;
	float err4 = 0.0;

	for (unsigned int index = 0; index < jets.size(); ++index) {
		// Info on this jet
		const pat::Jet& jet (jets.at(index));
		const unsigned int partonFlavour = abs( jet.partonFlavour() );
		const bool isBTagged = jet.bDiscriminator("combinedSecondaryVertexBJetTags") >= 0.679;

		// Get scale factor for this jet
		std::vector<double> sfAndError = getScaleFactor( partonFlavour, jet, MCSampleTag );
		const double sf = sfAndError.at(0);
		const double sfError = sfAndError.at(1);

		// Get efficiency for this jet
		const double eff = getEfficiency( partonFlavour, jet, MCSampleTag );
		if ( isBTagged ) {
			bTaggedMCJet *= eff;
			bTaggedDataJet *= eff*sf;
			if(partonFlavour==5 || partonFlavour ==4) err1 += sfError/sf; ///correlated for b/c
			else err3 += sfError/sf; //correlated for light
		} else {
			nonBTaggedMCJet *= ( 1 - eff );
			nonBTaggedDataJet *= ( 1 - eff*sf );
			if(partonFlavour==5 || partonFlavour ==4 ) err2 += (-eff*sfError)/(1-eff*sf); /// /correlated for b/c
			else err4 += (-eff*sfError)/(1-eff*sf); ////correlated for light
		}
	}

	double bTagWeight = (nonBTaggedDataJet * bTaggedDataJet) / (nonBTaggedMCJet * bTaggedMCJet);
	double error_BTagWeight = sqrt(pow(err1 + err2, 2) + pow(err3 + err4, 2)) * bTagWeight; ///un-correlated for b/c and light


	if (BTagSystematic == +1) {
		return bTagWeight + error_BTagWeight;
	} else if (BTagSystematic == -1) {
		return bTagWeight - error_BTagWeight;
	} else if (BTagSystematic == 0) {
		return bTagWeight;
	}
	return 0.;
}