10 {
11
12 std::vector<Gep::Cluster> constituents =
clusters;
13
14 std::vector<Gep::Jet> v_jets;
15
16 const unsigned int n_constituents = constituents.size();
17 std::vector< std::vector<int> > v_constitutents_in_jet_indices(n_constituents, std::vector<int>() );
18
19 for(
unsigned int i = 0;
i < n_constituents;
i++){ v_constitutents_in_jet_indices.at(i).push_back(i); }
20
21 int niter = 0;
22
23
24 while(!constituents.empty() && niter <
m_nIter){
25
26 const unsigned int n_clusters = constituents.size();
27
28
29 std::vector< std::vector< float > > DeltaR2(n_clusters, std::vector< float >(n_clusters, 0.0));
30 std::vector< float > InvPt2 (n_clusters, 0.0);
31 for (
unsigned int i = 0;
i < n_clusters;
i++){
32 InvPt2[
i] = 1/constituents.at(i).vec.Perp2();
34 for(unsigned int j = i+1; j < n_clusters; j++){
35 float deltaEta = constituents.at(i).vec.Eta() - constituents.at(j).vec.Eta();
36 float deltaPhi = constituents.at(i).vec.DeltaPhi(constituents.at(j).vec);
39 }
40 }
41
42
43
44 std::vector< int > minDeltaR2_indices(n_clusters, 0.0);
45 for (
unsigned int i = 0;
i < n_clusters;
i++ ) {
46 minDeltaR2_indices[
i] = std::distance( DeltaR2[i].
begin(), std::min_element(std::begin(DeltaR2[i]), std::end(DeltaR2[i])) );
47 }
48
49
50 std::vector< float > InvPt2DeltaR2(n_clusters, 0.0);
51 for (
unsigned int i = 0;
i < n_clusters;
i++ ) InvPt2DeltaR2[i] = InvPt2[i]*DeltaR2[i][minDeltaR2_indices[i]];
52
53
54 int minInvPt2DeltaR2_index = std::distance(InvPt2DeltaR2.begin() , std::min_element(std::begin(InvPt2DeltaR2), std::end(InvPt2DeltaR2)) );
55
56
57 if (InvPt2DeltaR2[minInvPt2DeltaR2_index] < InvPt2[minInvPt2DeltaR2_index]){
58
59
60 TVector3 momentum_vector(constituents.at(minInvPt2DeltaR2_index).vec.Px() + constituents.at(minDeltaR2_indices[minInvPt2DeltaR2_index]).vec.Px(),
61 constituents.at(minInvPt2DeltaR2_index).vec.Py() + constituents.at(minDeltaR2_indices[minInvPt2DeltaR2_index]).vec.Py(),
62 constituents.at(minInvPt2DeltaR2_index).vec.Pz() + constituents.at(minDeltaR2_indices[minInvPt2DeltaR2_index]).vec.Pz());
63
64
65 constituents.at(minInvPt2DeltaR2_index).vec.SetPtEtaPhiE(momentum_vector.Pt(),
66 constituents.at(minInvPt2DeltaR2_index).vec.Eta(),
67 constituents.at(minInvPt2DeltaR2_index).vec.Phi(),
68 momentum_vector.Mag() );
69
70
71 v_constitutents_in_jet_indices.at(minInvPt2DeltaR2_index).insert( v_constitutents_in_jet_indices.at(minInvPt2DeltaR2_index).end(),
72 v_constitutents_in_jet_indices.at(minDeltaR2_indices[minInvPt2DeltaR2_index]).begin(),
73 v_constitutents_in_jet_indices.at(minDeltaR2_indices[minInvPt2DeltaR2_index]).end() );
74
75 constituents.erase(constituents.begin() + minDeltaR2_indices[minInvPt2DeltaR2_index]);
76 v_constitutents_in_jet_indices.erase(v_constitutents_in_jet_indices.begin() + minDeltaR2_indices[minInvPt2DeltaR2_index]);
77 } else {
78
79
80 Gep::Jet jet;
81 jet.
vec.SetXYZT(constituents.at(minInvPt2DeltaR2_index).vec.Px(), constituents.at(minInvPt2DeltaR2_index).vec.Py(),
82 constituents.at(minInvPt2DeltaR2_index).vec.Pz(), constituents.at(minInvPt2DeltaR2_index).vec.E());
84 v_constitutents_in_jet_indices.at( minInvPt2DeltaR2_index).begin(),
85 v_constitutents_in_jet_indices.at( minInvPt2DeltaR2_index).end());
86
87 v_jets.push_back(std::move(jet));
88
89 constituents.erase(constituents.begin() + minInvPt2DeltaR2_index);
90 v_constitutents_in_jet_indices.erase(v_constitutents_in_jet_indices.begin() + minDeltaR2_indices[minInvPt2DeltaR2_index]);
91
92 }
93
94 niter++;
95 }
96
97 return v_jets;
98}
Scalar deltaPhi(const MatrixBase< Derived > &vec) const
double deltaEta(const I4Momentum &p1, const I4Momentum &p2)
Computes efficiently .
std::vector< int > constituentsIndices
1.14.0