C++ Jet::Pt方法代码示例

Jet MuXboostedBTagging::Result(std::vector<Lepton> const &muons,  Jet const &jet, Jet &core) const
{
    INFO0;
    // Keep track of the smallest x for any muon
    auto min_x = std::numeric_limits<double>::max(); // Large enough, simple significand
    // Iterate through each muon, add it to the core, and calculate the resulting boost invariant
    auto p4_neutrino_correction = Jet {};
    for (auto const &muon : muons) {
        // Add back the muon and the neutrino
        core += 2 * muon;
        // Set a hard ceiling on subjet mass, for poor reconstruction
        auto dot = muon.Spatial().Dot(core.Spatial());
        auto min = std::min(core.Mass(), MaxSubJetMass());
        if (dot == 0_eV * eV || min == 0_eV) continue;
        auto cross = muon.Spatial().Cross(core.Spatial()).Mag();
        auto x_core = static_cast<double>(core.Energy() * cross / dot / min);
        CHECK(x_core >= 0, x_core);
        if (x_core < 0) return jet;
        // Add the neutrino to original jet IFF the muon passes the boosted test
        if (x_core <= XMax()) p4_neutrino_correction += muon;
        min_x = std::min(min_x, x_core);
    }

    auto result = jet + p4_neutrino_correction;
    CHECK(result.Pt() > 0_eV, result.Pt());
    CHECK(!(result.Pt() != result.Pt()), result.Pt())
    CHECK(!(core.Pt() != core.Pt()), core.Pt())
    CHECK(!(min_x != min_x), min_x)
    if (result.Pt() <= 0_eV) return jet;
    result.SetInfo(jet.Info());
    result.Info().SetMuBTag(min_x, core.Pt() / result.Pt());
    return result;
}

double BTagEffService::GetEfficiency(BTagger const &bTagger, Jet const &jet)  const
{
    return GetEfficiency(bTagger, jet.Pt(), jet.Eta(), jet.Flavour(Jet::FlavourType::Hadron));
}

void fillPlots ( CollectionPtr & jets, 
		 CollectionPtr & core_jets, 
		 CollectionPtr & partons,
		 CollectionPtr & hgceeClusters,
		 CollectionPtr & hgchefClusters,
		 CollectionPtr & hgchebClusters,
		 likelihoodGetter & l, 
		 int index,
		 std::vector<TH1F*> quark_th1, std::vector<TH1F*> gluon_th1,
		 std::vector<TH2F*> quark_th2, std::vector<TH2F*> gluon_th2 ){

  int n_jets = jets -> GetSize();

  std::vector<double> likelihood_variables;

  for (int i = 0; i < n_jets; ++i){

    Jet jet = jets -> GetConstituent<Jet>(i);
    CollectionPtr all_core_jets     = core_jets;
    CollectionPtr matched_core_jets = core_jets -> SkimByRequireDRMatch<PFCoreJet, PFPrunedJet>( jet, 0.3 ) ;
    
    rechit_type my_type = NO_TYPE;
    int my_index = -1;
    jet.getLeadClusterTypeAndIndex(my_type, my_index);    
    double clusterLength = -1.;
    double clusterVolume = -1.;
    if ( my_type == NO_TYPE ) continue;

    if ( my_type == HGCEE ){
      HGCEECluster c(*hgceeClusters, my_index);
      clusterLength = c.getLength();
      clusterVolume = c.getVolume();
    }

    else if ( my_type == HGCHEF ){
      HGCHEFCluster c(*hgchefClusters, my_index);
      clusterLength = c.getLength();
      clusterVolume = c.getVolume();
    }

    else {
      HGCHEBCluster c(*hgchebClusters, my_index);
      clusterLength = c.getLength();
      clusterVolume = c.getVolume();
    }

    double core_over_total_pt = -1.;
    double lead_core_jet_dr   = -1.;
    double n_core_jets        = matched_core_jets -> GetSize();
    if ( n_core_jets > 0 ){
      PFCoreJet lead_core_jet = matched_core_jets -> GetLeadPtObject<PFCoreJet> ();
      double core_pt = lead_core_jet.Pt();
      core_over_total_pt = core_pt / jet.Pt();
      lead_core_jet_dr   = jet.DeltaR( &lead_core_jet );
    }
    
    int flavor = getFlavor ( jet , partons );
    if ( flavor == -1 ) continue;
    bool isQuarkJet ( flavor == 1 );
    bool isGluonJet ( flavor == 0 );


    likelihood_variables.clear();

    double weighted_depth = jet.getWeightedDepth();
    double weighted_depth_noEE = jet.getWeightedDepthNoEE();

    double profile[n_radii];
    jet.getProfile(n_radii, radii, profile);

    double profile50 = jet.getProfileRadius(n_radii, radii, profile, 0.5);
    double profile90 = jet.getProfileRadius(n_radii, radii, profile, 0.9);
    double profile95 = jet.getProfileRadius(n_radii, radii, profile, 0.95);
    double profile99 = jet.getProfileRadius(n_radii, radii, profile, 0.99);

    likelihood_variables.push_back ( jet.TrimmedNSubjets(index) );
    likelihood_variables.push_back ( jet.getWidth()   );
    likelihood_variables.push_back ( jet.getNPFCandidates() );
    likelihood_variables.push_back ( jet.getPTD() );
    likelihood_variables.push_back ( jet.NSubJettiness() );
    likelihood_variables.push_back ( jet.TrimmedMass(index) * 1000. );
    likelihood_variables.push_back ( profile50 );
    likelihood_variables.push_back ( core_over_total_pt );
    likelihood_variables.push_back ( lead_core_jet_dr );

    double ee_energy     = jet.getEEEnergy();
    double ee015_energy  = jet.getEEEnergy(0, 15);
    double ee1631_energy = jet.getEEEnergy(16, 31);
    double hef_energy    = jet.getHEFEnergy();
    double heb_energy    = jet.getHEBEnergy();
    double all_energy    = ee_energy + hef_energy + heb_energy;

    double hef11_energy  = jet.getHEFEnergy(1, 1);
    double hef22_energy  = jet.getHEFEnergy(2, 2);
    double hef34_energy  = jet.getHEFEnergy(3, 4);
    double hef56_energy  = jet.getHEFEnergy(5, 6);
    double hef712_energy = jet.getHEFEnergy(7,12);

    double likelihood  = l.getLikelihood ( likelihood_variables );

//.........这里部分代码省略.........

void JESBDTVars::FillBranches(EventContainer * evtObj){

  //Evaluate each distribution once per JES shift
  for (int i = 0; i < evtObj->jets[0].GetNumberOfJESCorrections(); i++){
    //First let's clear the things we already have
    selectedJet.clear();
    Jet2040.clear();
    BJet.clear();
    UntaggedJet.clear();

    //And make a dummy variable here
    Jet tempJet;
    TLorentzVector tempjet(0,0,0,0);

    //Number of loose jets
    for (auto jet : evtObj->alljets){
      tempJet = jet;
      TLorentzVector tempmet(0,0,0,0);
      tempJet.ShiftPtWithJESCorr(i,&tempmet);
      if (tempJet.Pt() > 20 && tempJet.Pt() < 40) Jet2040.push_back(tempJet);
    }

    for (auto jet : evtObj->jesShiftedJets[i]){
      selectedJet.push_back(jet);
      if (jet.IsTagged()) BJet.push_back(jet);
      else UntaggedJet.push_back(jet);
    }

    //Now set up the lepton and met variables
    TLorentzVector Lepton(00,0,0,0);
    TLorentzVector Miss(00,0,0,0);
    TLorentzVector Wlv(0,0,0,0);
    
    Miss = evtObj->metVecsJESShifted[i];

    if (evtObj->electronsToUsePtr->size() > 0){ // if this number is >0 we're in the electron channel. Otherwise use muons
      Lepton.SetPtEtaPhiE(evtObj->electronsToUsePtr->at(0).Pt(),evtObj->electronsToUsePtr->at(0).Eta(),evtObj->electronsToUsePtr->at(0).Phi(),evtObj->electronsToUsePtr->at(0).E());
    }
    else{
      Lepton.SetPtEtaPhiE(evtObj->muonsToUsePtr->at(0).Pt(),evtObj->muonsToUsePtr->at(0).Eta(),evtObj->muonsToUsePtr->at(0).Phi(),evtObj->muonsToUsePtr->at(0).E());
    }
    Wlv = Lepton+Miss;
   
    TLorentzVector W(0,0,0,0), Top(0,0,0,0);
    for (auto jet : UntaggedJet){
      W = W + jet;
    }

    TLorentzVector totalJets(0,0,0,0);
    for (auto jet : selectedJet){
      totalJets += jet;
    }

    //That should be all of the things we need to make the BDT variables, so let's make the variables now.
    _floatVecVars["M_DeltaRBJetLepton_JESShifts"][i] = fabs(BJet[0].DeltaR(Lepton));
    if (UntaggedJet.size() > 1) _floatVecVars["M_DeltaRlightjets_JESShifts"][i] = UntaggedJet.at(0).DeltaR(UntaggedJet.at(1));
    
    if (BJet[0].DeltaR(Lepton) > BJet[0].DeltaR(W)){
      _floatVecVars["M_topMass2_lep_JESShifts"][i] = -1;
    }
    else{
      Top = W + BJet[0];
      _floatVecVars["M_topMass2_lep_JESShifts"][i] = Top.M();
    }
    _floatVecVars["M_Pt_Lepton_JESShifts"][i] = Lepton.Pt();
  
  _floatVecVars["M_Pt_AllJetsLeptonMET_JESShifts"][i] = (Lepton + Miss + totalJets).Pt();

  _floatVecVars["M_DeltaRLeptonJet1_JESShifts"][i] = fabs(selectedJet.at(0).DeltaR(Lepton));
  
  if (selectedJet.size() > 2) _floatVecVars["M_Mass_Jet1Jet2Jet3LeptonMET_JESShifts"][i] = (selectedJet[0] + selectedJet[1] + selectedJet[2]).M();

  _floatVecVars["M_hadronicWmass_JESShifts"][i] = W.M();

  _floatVecVars["lightJet1CSV_JESShifts"][i] = UntaggedJet[0].GetbDiscriminator();

  }

  


}

std::vector<Jet> MuXboostedBTagging::CoreCandidates(std::vector<Lepton> const &muons, Jet const &jet) const
{
    INFO0;
    auto recluster_input = muons;
    for (auto const &consituent : jet.Constituents()) {
        if (consituent.Info().ContainsDetectorPart(DetectorPart::tower) && consituent.Pt() < min_tower_pt_ratio_ * jet.Pt()) continue; // Don't use
        recluster_input.emplace_back(consituent);
    }
    // Recluster the jet, to find core candidates
    auto cluster_sequence = boca::ClusterSequence {recluster_input, fastjet::JetDefinition(fastjet::antikt_algorithm, core_jet_radius / rad, &Settings::Recombiner())};
    // Get the core candidates (NOT sorted by pT until we remove muons)
    auto core_candidates = cluster_sequence.InclusiveJets();
    if (core_candidates.empty()) return core_candidates;
    cluster_sequence.NoLongerNeeded();

    // If a taggable muon is inside a core candidate, remove the muon p4
    for (auto const &muon : muons) for (auto &core_candidate : core_candidates) if (muon.DeltaRTo(core_candidate) < core_jet_radius) {
                core_candidate -= muon;
                break;
            }

    // Sort the core candidates by pT (highest to lowest)
    core_candidates = SortedByPt(core_candidates);
    return core_candidates;
}