44{
45 if (GenStableCharged.empty()) {
47 return nullptr;
48 }
49 std::vector< Trk::GenParticleJet >* myJets
50 = new std::vector< Trk::GenParticleJet >;
51 ATH_MSG_DEBUG (
"start jetMCFinder() with input GenParticle-vector "
52 << "of size " << GenStableCharged.size() );
53
54 myJets->clear();
55 double totalEnergyFromTracks=0.;
56 int partIterations=0;
57 double yijmin=100.;
58 bool combineParticles = false;
59
60 std::vector<Trk::GenParticleJet>::iterator iAtMin, jAtMin;
61
62 for( auto i = GenStableCharged.begin() ; i < GenStableCharged.end(); ++i){
63 totalEnergyFromTracks = totalEnergyFromTracks + (*i)->momentum().e();
64 Trk::GenParticleJet tempPJ;
65 tempPJ.
addParticle( *i,
int(i - GenStableCharged.begin()) );
66 myJets->push_back(tempPJ);
67 }
68
69 int numParts = myJets->size();
70 do {
71 combineParticles = false;
72 yijmin=100.;
73 for( std::vector<Trk::GenParticleJet>::iterator i = myJets->begin() ; i<myJets->
end(); ++
i) {
74
75 for( std::vector<Trk::GenParticleJet>::iterator j = myJets->begin() ; j<myJets->
end(); ++j) {
76
77 if (i!=j) {
78
79 double yij=0;
80 double absI= sqrt((
i->getMomentum()).dot(
i->getMomentum()));
81 double absJ= sqrt(((*j).getMomentum()).dot((*j).getMomentum()));
82 double IdotJ= (((*i).getMomentum()).
dot((*j).getMomentum()));
83 if (fabs(absI*absJ)<0.000001 || fabs(totalEnergyFromTracks)<0.00001) {
84
86 << "can't continue calculations." );
87 delete myJets;
88 return nullptr;
89
90 } else {
91 yij=2*(
i->getEnergy())*(j->getEnergy())/
pow(totalEnergyFromTracks,2)*(1-(IdotJ/(absI*absJ)));
92 if (yij<yijmin) {
93
94 combineParticles=true;
95 yijmin=yij;
97 jAtMin=j;
98
99 };
100 };
101 }
102 }
103 }
104 partIterations++;
105
106
107
108 if (combineParticles) {
109
111 << " ("<<(*iAtMin).getIndicesInEvent().size()<<")"
112 << " and "<< (*jAtMin).getNumParticles()
113 << " ("<<(*jAtMin).getIndicesInEvent().size()<<")" );
114 if((*iAtMin).getNumParticles()>(*jAtMin).getNumParticles()) {
115
116 auto partsTemp = (*jAtMin).getParticles();
117 std::vector<int> indexTemp = (*jAtMin).getIndicesInEvent();
118
119 if (!partsTemp.empty()) {
120
121 std::vector<int>::iterator ki=indexTemp.begin();
122 for (auto k =partsTemp.begin(); k!=partsTemp.end(); ++k, ++ki) {
123 (*iAtMin).addParticle(*k,*ki);
124 }
125 } else {
126
127 ATH_MSG_WARNING (
"No particles in this jet, logic failed, stop finder." );
128 delete myJets;
129 return nullptr;
130 }
131 myJets->erase(jAtMin);
132 }
133 else {
134
135 auto partsTemp = (*iAtMin).getParticles();
136 std::vector<int> indexTemp = (*iAtMin).getIndicesInEvent();
137 if (!partsTemp.empty()) {
138
139 std::vector<int>::iterator ki=indexTemp.begin();
140 for (auto k =partsTemp.begin(); k!=partsTemp.end(); ++k, ++ki) {
141
142 (*jAtMin).addParticle(*k,*ki);
143 }
144 }
145 else {
146 ATH_MSG_WARNING (
"No particles in this jet, logic failed, stop finder." );
147 delete myJets;
148 return nullptr;
149 }
150 myJets->erase(iAtMin);
151 }
152 }
153 int currentNumParts = 0;
154 for(std::vector<Trk::GenParticleJet>::iterator k = myJets->begin(); k<myJets->
end(); ++
k)
155 currentNumParts = currentNumParts +
k->getNumParticles();
156
157 if(numParts != currentNumParts){
159 delete myJets;
160 return nullptr;
161 }
162 if(partIterations==999)
163 ATH_MSG_INFO (
"JetMCFinder cut off because of too many (1000) iterations");
164
165 }
while (yijmin<
m_yijCut && partIterations<1000);
166 return myJets;
167}
#define ATH_MSG_VERBOSE(x)
#define ATH_MSG_WARNING(x)
constexpr int pow(int base, int exp) noexcept
void addParticle(HepMC::ConstGenParticlePtr part, int indexInEvent=-1)
dot(G, fn, nodesToHighlight=[])
1.14.0