CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HemisphereAlgo.cc
Go to the documentation of this file.
1 #include "PhysicsTools/PatAlgos/interface/HemisphereAlgo.h"
2 
3 #include "FWCore/MessageLogger/interface/MessageLogger.h"
4 
5 #include "DataFormats/Math/interface/deltaPhi.h"
6 #include "DataFormats/Math/interface/deltaR.h"
7 
8 using namespace std;
9 using std::vector;
10 
11 HemisphereAlgo::HemisphereAlgo(const std::vector<reco::CandidatePtr>& componentPtr_,const int seed_method, const int hemisphere_association_method )
12  : Object(componentPtr_), Object_Group() , Axis1(), Axis2(), seed_meth(seed_method), hemi_meth(hemisphere_association_method), status(0) {
13 
14  for(int i = 0; i < (int) Object.size(); i++){
15  Object_Noseed.push_back(0);
16  }
17 
18 }
19 
20 
21 
22 vector<float> HemisphereAlgo::getAxis1(){
23  if (status != 1){this->reconstruct();}
24  return Axis1;
25 }
26 vector<float> HemisphereAlgo::getAxis2(){
27  if (status != 1){this->reconstruct();}
28  return Axis2;
29 }
30 
31 vector<int> HemisphereAlgo::getGrouping(){
32  if (status != 1){this->reconstruct();}
33  return Object_Group;
34 }
35 
36 int HemisphereAlgo::reconstruct(){
37 
38  int vsize = (int) Object.size();
39 
40  LogDebug("HemisphereAlgo") << " HemisphereAlgo method ";
41 
42  Object_Group.clear();
43  Axis1.clear();
44  Axis2.clear();
45 
46  for(int j = 0; j < vsize; j++){
47  Object_Group.push_back(0);
48  }
49 
50  for(int j = 0; j < 5; j++){
51  Axis1.push_back(0);
52  Axis2.push_back(0);
53  }
54 
55 
56  for (int i = 0; i <vsize; i++){
57 
58  if ( (*(Object)[i]).p() > (*Object[i]).energy() + 0.001) {
59 
60  edm::LogWarning("HemisphereAlgo") << "Object " << i << " has E = " << (*Object[i]).energy()
61  << " less than P = " << (*Object[i]).p() ;
62 
63  }
64  }
65 
66 
67 
68  LogDebug("HemisphereAlgo") << " Seeding method = " << seed_meth;
69  int I_Max = -1;
70  int J_Max = -1;
71 
72  if (seed_meth == 1) {
73 
74  float P_Max = 0.;
75  float DeltaRP_Max = 0.;
76 
77  // take highest momentum object as first axis
78  for (int i=0;i<vsize;i++){
79  Object_Group[i] = 0;
80  if (Object_Noseed[i] == 0 && P_Max < (*Object[i]).p()){
81  P_Max = (*Object[i]).p();
82  I_Max = i;
83  }
84  }
85 
86  Axis1[0] = (*Object[I_Max]).px() / (*Object[I_Max]).p();
87  Axis1[1] = (*Object[I_Max]).py() / (*Object[I_Max]).p();
88  Axis1[2] = (*Object[I_Max]).pz() / (*Object[I_Max]).p();
89  Axis1[3] = (*Object[I_Max]).p();
90  Axis1[4] = (*Object[I_Max]).energy();
91 
92  // take as second axis
93  for (int i=0;i<vsize;i++){
94 
95  float DeltaR = deltaR((*Object[i]).eta(),(*Object[i]).phi(),(*Object[I_Max]).eta(),(*Object[I_Max]).phi()) ;
96 
97  if (Object_Noseed[i] == 0 && DeltaR > 0.5) {
98 
99  float DeltaRP = DeltaR * (*Object[i]).p();
100  if (DeltaRP > DeltaRP_Max){
101  DeltaRP_Max = DeltaRP;
102  J_Max = i;
103  }
104  }
105  }
106 
107  if (J_Max >=0){
108  Axis2[0] = (*Object[J_Max]).px() / (*Object[J_Max]).p();
109  Axis2[1] =(*Object[J_Max]).py() / (*Object[J_Max]).p();
110  Axis2[2] = (*Object[J_Max]).pz() / (*Object[J_Max]).p();
111  Axis2[3] = (*Object[J_Max]).p();
112  Axis2[4] = (*Object[J_Max]).energy();
113 
114  } else {
115  return 0;
116  }
117  LogDebug("HemisphereAlgo") << " Axis 1 is Object = " << I_Max;
118  LogDebug("HemisphereAlgo") << " Axis 2 is Object = " << J_Max;
119 
120 
121  } else if ( (seed_meth == 2) | (seed_meth == 3) ) {
122 
123  float Mass_Max = 0.;
124  float InvariantMass = 0.;
125 
126  // maximize the invariant mass of two objects
127  for (int i=0;i<vsize;i++){
128  Object_Group[i] = 0;
129  if (Object_Noseed[i] == 0){
130  for (int j=i+1;j<vsize;j++){
131  if (Object_Noseed[j] == 0){
132 
133  // either the invariant mass
134  if (seed_meth == 2){
135  InvariantMass = ((*Object[i]).energy() + (*Object[j]).energy())* ((*Object[i]).energy() + (*Object[j]).energy())
136  - ((*Object[i]).px() + (*Object[j]).px())* ((*Object[i]).px() + (*Object[j]).px())
137  - ((*Object[i]).py() + (*Object[j]).py())* ((*Object[i]).py() + (*Object[j]).py())
138  - ((*Object[i]).pz() + (*Object[j]).pz())* ((*Object[i]).pz() + (*Object[j]).pz()) ;
139  }
140  // or the transverse mass
141  else if (seed_meth == 3){
142  float pti = sqrt((*Object[i]).px()*(*Object[i]).px() + (*Object[i]).py()*(*Object[i]).py());
143  float ptj = sqrt((*Object[j]).px()*(*Object[j]).px() + (*Object[j]).py()*(*Object[j]).py());
144  InvariantMass = 2. * (pti*ptj - (*Object[i]).px()*(*Object[j]).px()
145  - (*Object[i]).py()*(*Object[j]).py() );
146  }
147  if ( Mass_Max < InvariantMass){
148  Mass_Max = InvariantMass;
149  I_Max = i;
150  J_Max = j;
151  }
152  }
153  }
154  }
155  }
156 
157  if (J_Max>0) {
158 
159  Axis1[0] = (*Object[I_Max]).px() / (*Object[I_Max]).p();
160  Axis1[1] = (*Object[I_Max]).py() / (*Object[I_Max]).p();
161  Axis1[2] = (*Object[I_Max]).pz() / (*Object[I_Max]).p();
162 
163  Axis1[3] = (*Object[I_Max]).p();
164  Axis1[4] = (*Object[I_Max]).energy();
165 
166  Axis2[0] = (*Object[J_Max]).px() / (*Object[J_Max]).p();
167  Axis2[1] =(*Object[J_Max]).py() / (*Object[J_Max]).p();
168  Axis2[2] = (*Object[J_Max]).pz() / (*Object[J_Max]).p();
169 
170  Axis2[3] = (*Object[J_Max]).p();
171  Axis2[4] = (*Object[J_Max]).energy();
172 
173  } else {
174  return 0;
175  }
176  LogDebug("HemisphereAlgo") << " Axis 1 is Object = " << I_Max;
177  LogDebug("HemisphereAlgo") << " Axis 2 is Object = " << J_Max;
178 
179 
180  } else {
181  throw cms::Exception("Configuration") << "Please give a valid seeding method!";
182  }
183 
184  // seeding done
185  // now do the hemisphere association
186 
187  LogDebug("HemisphereAlgo") << " Association method = " << hemi_meth;
188 
189  int numLoop = 0;
190  bool I_Move = true;
191 
192 
193  while (I_Move && (numLoop < 100)){
194 
195  I_Move = false;
196  numLoop++;
197  LogDebug("HemisphereAlgo") << " Iteration = " << numLoop;
198 
199  float Sum1_Px = 0.;
200  float Sum1_Py = 0.;
201  float Sum1_Pz = 0.;
202  float Sum1_P = 0.;
203  float Sum1_E = 0.;
204  float Sum2_Px = 0.;
205  float Sum2_Py = 0.;
206  float Sum2_Pz = 0.;
207  float Sum2_P = 0.;
208  float Sum2_E = 0.;
209 
210 
211  if (hemi_meth == 1) {
212 
213  for (int i=0;i<vsize;i++){
214  float P_Long1 = (*Object[i]).px()*Axis1[0] + (*Object[i]).py()*Axis1[1] + (*Object[i]).pz()*Axis1[2];
215  float P_Long2 = (*Object[i]).px()*Axis2[0]+ (*Object[i]).py()*Axis2[1] + (*Object[i]).pz()*Axis2[2];
216  if (P_Long1 >= P_Long2){
217  if (Object_Group[i] != 1){
218  I_Move = true;
219  }
220  Object_Group[i] = 1;
221  Sum1_Px += (*Object[i]).px();
222  Sum1_Py += (*Object[i]).py();
223  Sum1_Pz += (*Object[i]).pz();
224  Sum1_P += (*Object[i]).p();
225  Sum1_E += (*Object[i]).energy();
226  } else {
227  if (Object_Group[i] != 2){
228  I_Move = true;
229  }
230  Object_Group[i] = 2;
231  Sum2_Px += (*Object[i]).px();
232  Sum2_Py += (*Object[i]).py();
233  Sum2_Pz += (*Object[i]).pz();
234  Sum2_P += (*Object[i]).p();
235  Sum2_E += (*Object[i]).energy();
236  }
237  }
238 
239  } else if (hemi_meth == 2 || hemi_meth == 3){
240 
241  for (int i=0;i<vsize;i++){
242  if (i == I_Max) {
243  Object_Group[i] = 1;
244  Sum1_Px += (*Object[i]).px();
245  Sum1_Py += (*Object[i]).py();
246  Sum1_Pz += (*Object[i]).pz();
247  Sum1_P += (*Object[i]).p();
248  Sum1_E += (*Object[i]).energy();
249  } else if (i == J_Max) {
250  Object_Group[i] = 2;
251  Sum2_Px += (*Object[i]).px();
252  Sum2_Py += (*Object[i]).py();
253  Sum2_Pz += (*Object[i]).pz();
254  Sum2_P += (*Object[i]).p();
255  Sum2_E += (*Object[i]).energy();
256  } else {
257 
258 
259  if(!I_Move){
260 
261  float NewAxis1_Px = Axis1[0] * Axis1[3];
262  float NewAxis1_Py = Axis1[1] * Axis1[3];
263  float NewAxis1_Pz = Axis1[2] * Axis1[3];
264  float NewAxis1_E = Axis1[4];
265 
266  float NewAxis2_Px = Axis2[0] * Axis2[3];
267  float NewAxis2_Py = Axis2[1] * Axis2[3];
268  float NewAxis2_Pz = Axis2[2] * Axis2[3];
269  float NewAxis2_E = Axis2[4];
270 
271  if (Object_Group[i] == 1){
272 
273  NewAxis1_Px = NewAxis1_Px - (*Object[i]).px();
274  NewAxis1_Py = NewAxis1_Py - (*Object[i]).py();
275  NewAxis1_Pz = NewAxis1_Pz - (*Object[i]).pz();
276  NewAxis1_E = NewAxis1_E - (*Object[i]).energy();
277 
278  } else if (Object_Group[i] == 2) {
279 
280  NewAxis2_Px = NewAxis2_Px - (*Object[i]).px();
281  NewAxis2_Py = NewAxis2_Py - (*Object[i]).py();
282  NewAxis2_Pz = NewAxis2_Pz - (*Object[i]).pz();
283  NewAxis2_E = NewAxis2_E - (*Object[i]).energy();
284  }
285 
286 
287  float mass1 = NewAxis1_E - (((*Object[i]).px()*NewAxis1_Px + (*Object[i]).py()*NewAxis1_Py +
288  (*Object[i]).pz()*NewAxis1_Pz)/(*Object[i]).p());
289 
290  float mass2 = NewAxis2_E - (((*Object[i]).px()*NewAxis2_Px + (*Object[i]).py()*NewAxis2_Py +
291  (*Object[i]).pz()*NewAxis2_Pz)/(*Object[i]).p());
292 
293  if (hemi_meth == 3) {
294 
295  mass1 *= NewAxis1_E/((NewAxis1_E+(*Object[i]).energy())*(NewAxis1_E+(*Object[i]).energy()));
296 
297  mass2 *= NewAxis2_E/((NewAxis2_E+(*Object[i]).energy())*(NewAxis2_E+(*Object[i]).energy()));
298 
299  }
300 
301  if(mass1 < mass2) {
302  if (Object_Group[i] != 1){
303  I_Move = true;
304  }
305  Object_Group[i] = 1;
306 
307  Sum1_Px += (*Object[i]).px();
308  Sum1_Py += (*Object[i]).py();
309  Sum1_Pz += (*Object[i]).pz();
310  Sum1_P += (*Object[i]).p();
311  Sum1_E += (*Object[i]).energy();
312  } else {
313  if (Object_Group[i] != 2){
314  I_Move = true;
315  }
316  Object_Group[i] = 2;
317  Sum2_Px += (*Object[i]).px();
318  Sum2_Py += (*Object[i]).py();
319  Sum2_Pz += (*Object[i]).pz();
320  Sum2_P += (*Object[i]).p();
321  Sum2_E += (*Object[i]).energy();
322 
323  }
324 
325 
326  } else {
327 
328  if (Object_Group[i] == 1){
329  Sum1_Px += (*Object[i]).px();
330  Sum1_Py += (*Object[i]).py();
331  Sum1_Pz += (*Object[i]).pz();
332  Sum1_P += (*Object[i]).p();
333  Sum1_E += (*Object[i]).energy();
334  }
335  if (Object_Group[i] == 2){
336  Sum2_Px += (*Object[i]).px();
337  Sum2_Py += (*Object[i]).py();
338  Sum2_Pz += (*Object[i]).pz();
339  Sum2_P += (*Object[i]).p();
340  Sum2_E += (*Object[i]).energy();
341  }
342 
343 
344 
345  }
346 
347 
348  }
349  }
350  } else {
351  throw cms::Exception("Configuration")
352  << "Please give a valid hemisphere association method!";
353  }
354 
355  // recomputing the axes
356 
357  Axis1[3] = sqrt(Sum1_Px*Sum1_Px + Sum1_Py*Sum1_Py + Sum1_Pz*Sum1_Pz);
358  if (Axis1[3]<0.0001) {
359  edm::LogWarning("HemisphereAlgo") << "ZERO objects in group 1! ";
360  } else {
361  Axis1[0] = Sum1_Px / Axis1[3];
362  Axis1[1] = Sum1_Py / Axis1[3];
363  Axis1[2] = Sum1_Pz / Axis1[3];
364  Axis1[4] = Sum1_E;
365  }
366 
367 
368 
369  Axis2[3] = sqrt(Sum2_Px*Sum2_Px + Sum2_Py*Sum2_Py + Sum2_Pz*Sum2_Pz);
370  if (Axis2[3]<0.0001) {
371  edm::LogWarning("HemisphereAlgo") << "ZERO objects in group 2!";
372  } else {
373  Axis2[0] = Sum2_Px / Axis2[3];
374  Axis2[1] = Sum2_Py / Axis2[3];
375  Axis2[2] = Sum2_Pz / Axis2[3];
376  Axis2[4] = Sum2_E;
377  }
378 
379  LogDebug("HemisphereAlgo") << " Grouping = ";
380  for (int i=0;i<vsize;i++){
381  LogTrace("HemisphereAlgo") << " " << Object_Group[i];
382  }
383  LogTrace("HemisphereAlgo") << std::endl;
384 
385  }
386 
387 
388  status = 1;
389  return 1;
390 }
391 
392 
393 
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:501
i
int i
Definition: DBlmapReader.cc:9
HemisphereAlgo::Object_Noseed
std::vector< int > Object_Noseed
Definition: HemisphereAlgo.h:86
HemisphereAlgo::Axis2
std::vector< float > Axis2
Definition: HemisphereAlgo.h:89
DeltaR
Definition: deltaR.h:79
HemisphereAlgo::Axis1
std::vector< float > Axis1
Definition: HemisphereAlgo.h:88
HemisphereAlgo::Object_Group
std::vector< int > Object_Group
Definition: HemisphereAlgo.h:85
eta
T eta() const
Definition: Basic3DVectorLD.h:179
HemisphereAlgo::getAxis1
std::vector< float > getAxis1()
Definition: HemisphereAlgo.cc:22
HemisphereAlgo::HemisphereAlgo
HemisphereAlgo(const std::vector< reco::CandidatePtr > &componentRefs_, const int seed_method=0, const int hemisphere_association_method=0)
Definition: HemisphereAlgo.cc:11
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:48
j
int j
Definition: DBlmapReader.cc:9
LogTrace
#define LogTrace(id)
Definition: MessageLogger.h:502
HemisphereAlgo::getGrouping
std::vector< int > getGrouping()
Definition: HemisphereAlgo.cc:31
deltaR
double deltaR(double eta1, double eta2, double phi1, double phi2)
Definition: TreeUtility.cc:17
HemisphereAlgo::Object
std::vector< reco::CandidatePtr > Object
Definition: HemisphereAlgo.h:83
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
relval_parameters_module.energy
string energy
Definition: relval_parameters_module.py:29
ntuplemaker.status
tuple status
Definition: ntuplemaker.py:245
HemisphereAlgo::getAxis2
std::vector< float > getAxis2()
Definition: HemisphereAlgo.cc:26
phi
Definition: DDAxes.h:10
Generated for CMSSW Reference Manual by doxygen 1.8.5