CMS 3D CMS Logo
List of all members | Public Member Functions | Private Member Functions | Private Attributes
heppy::Hemisphere Class Reference

#include <Hemisphere.h>

Public Member Functions

void ClearAllNoLists ()
 
std::vector< float > getAxis1 ()
 
std::vector< float > getAxis2 ()
 
std::vector< int > getGrouping ()
 
int GetNumLoop ()
 
 Hemisphere ()
 
 Hemisphere (std::vector< float > Px_vector, std::vector< float > Py_vector, std::vector< float > Pz_vector, std::vector< float > E_vector, int seed_method, int hemisphere_association_method)
 
 Hemisphere (std::vector< float > Px_vector, std::vector< float > Py_vector, std::vector< float > Pz_vector, std::vector< float > E_vector)
 
void RejectISRDRmax (float drmax)
 
void RejectISRPtmax (float ptmax)
 
void SetDebug (int debug)
 
void SetDRminSeed1 (float rmin)
 
void SetMethod (int seed_method, int hemisphere_association_method)
 
void SetnItermax (int niter)
 
void SetNoAssoc (int object_number)
 
void SetNoSeed (int object_number)
 
 ~Hemisphere ()
 

Private Member Functions

int Reconstruct ()
 
int RejectISR ()
 

Private Attributes

std::vector< float > Axis1
 
std::vector< float > Axis2
 
int dbg
 
float dRminSeed1
 
int hemi_meth
 
int nItermax
 
int numLoop
 
std::vector< float > Object_E
 
std::vector< float > Object_Eta
 
std::vector< int > Object_Group
 
std::vector< int > Object_Noassoc
 
std::vector< int > Object_Noseed
 
std::vector< float > Object_P
 
std::vector< float > Object_Phi
 
std::vector< float > Object_Pt
 
std::vector< float > Object_Px
 
std::vector< float > Object_Py
 
std::vector< float > Object_Pz
 
int rejectISR
 
int rejectISRDR
 
float rejectISRDRmax
 
int rejectISRPt
 
float rejectISRPtmax
 
int seed_meth
 
int status
 

Detailed Description

Definition at line 28 of file Hemisphere.h.

Constructor & Destructor Documentation

heppy::Hemisphere::Hemisphere ( )
inline

Definition at line 75 of file Hemisphere.h.

75 {};
heppy::Hemisphere::Hemisphere ( std::vector< float >  Px_vector,
std::vector< float >  Py_vector,
std::vector< float >  Pz_vector,
std::vector< float >  E_vector,
int  seed_method,
int  hemisphere_association_method 
)

Definition at line 14 of file Hemisphere.cc.

References mps_fire::i, createfilelist::int, numLoop, Object_Noassoc, Object_Noseed, and Object_Px.

15  : Object_Px(Px_vector),
16  Object_Py(Py_vector), Object_Pz(Pz_vector), Object_E(E_vector), seed_meth(seed_method),
17  hemi_meth(hemisphere_association_method), status(0),
18  dRminSeed1(0.5),nItermax(100),
19  rejectISR(0), rejectISRPt(0), rejectISRPtmax(10000.),
20 rejectISRDR(0), rejectISRDRmax(100.), dbg(0) {
21  for(int i = 0; i < (int) Object_Px.size(); i++){
22  Object_Noseed.push_back(0);
23  Object_Noassoc.push_back(0);
24  }
25  numLoop =0;
26 }
mps_fire.i
i
Definition: mps_fire.py:156
heppy::Hemisphere::Object_Py
std::vector< float > Object_Py
Definition: Hemisphere.h:172
heppy::Hemisphere::Object_Noseed
std::vector< int > Object_Noseed
Definition: Hemisphere.h:180
heppy::Hemisphere::Object_E
std::vector< float > Object_E
Definition: Hemisphere.h:176
heppy::Hemisphere::Object_Noassoc
std::vector< int > Object_Noassoc
Definition: Hemisphere.h:181
heppy::Hemisphere::Object_Pz
std::vector< float > Object_Pz
Definition: Hemisphere.h:173
heppy::Hemisphere::Object_Px
std::vector< float > Object_Px
Definition: Hemisphere.h:171
heppy::Hemisphere::Hemisphere ( std::vector< float >  Px_vector,
std::vector< float >  Py_vector,
std::vector< float >  Pz_vector,
std::vector< float >  E_vector 
)

Definition at line 30 of file Hemisphere.cc.

References mps_fire::i, createfilelist::int, numLoop, Object_Noassoc, Object_Noseed, and Object_Px.

31  : Object_Px(Px_vector),
32  Object_Py(Py_vector), Object_Pz(Pz_vector), Object_E(E_vector), seed_meth(0),
33  hemi_meth(0), status(0),
34  dRminSeed1(0.5),nItermax(100),
35  rejectISR(0), rejectISRPt(0), rejectISRPtmax(10000.),
36 rejectISRDR(0), rejectISRDRmax(100.), dbg(0) {
37  for(int i = 0; i < (int) Object_Px.size(); i++){
38  Object_Noseed.push_back(0);
39  Object_Noassoc.push_back(0);
40  }
41  numLoop =0;
42 }
mps_fire.i
i
Definition: mps_fire.py:156
heppy::Hemisphere::Object_Py
std::vector< float > Object_Py
Definition: Hemisphere.h:172
heppy::Hemisphere::Object_Noseed
std::vector< int > Object_Noseed
Definition: Hemisphere.h:180
heppy::Hemisphere::Object_E
std::vector< float > Object_E
Definition: Hemisphere.h:176
heppy::Hemisphere::Object_Noassoc
std::vector< int > Object_Noassoc
Definition: Hemisphere.h:181
heppy::Hemisphere::Object_Pz
std::vector< float > Object_Pz
Definition: Hemisphere.h:173
heppy::Hemisphere::Object_Px
std::vector< float > Object_Px
Definition: Hemisphere.h:171
heppy::Hemisphere::~Hemisphere ( )
inline

Definition at line 82 of file Hemisphere.h.

References getAxis1(), getAxis2(), and getGrouping().

82 {};

Member Function Documentation

void heppy::Hemisphere::ClearAllNoLists ( )
inline

Definition at line 121 of file Hemisphere.h.

References mps_fire::i, createfilelist::int, Object_Noassoc, Object_Noseed, and status.

121  {
122  for (int i = 0; i < (int)Object_Noseed.size(); ++i) {
123  Object_Noassoc[i] = 0;
124  Object_Noseed[i] = 0;
125  status = 0;
126  }
127  }
mps_fire.i
i
Definition: mps_fire.py:156
heppy::Hemisphere::Object_Noseed
std::vector< int > Object_Noseed
Definition: Hemisphere.h:180
heppy::Hemisphere::Object_Noassoc
std::vector< int > Object_Noassoc
Definition: Hemisphere.h:181
vector< float > heppy::Hemisphere::getAxis1 ( )

Definition at line 44 of file Hemisphere.cc.

References Axis1, Reconstruct(), RejectISR(), rejectISR, and status.

Referenced by ~Hemisphere().

44  {
45  if (status != 1) {
46  if (rejectISR == 0) {
47  this->Reconstruct();
48  } else {
49  this->RejectISR();
50  }
51  }
52  return Axis1;
53 }
heppy::Hemisphere::Axis1
std::vector< float > Axis1
Definition: Hemisphere.h:183
vector< float > heppy::Hemisphere::getAxis2 ( )

Definition at line 54 of file Hemisphere.cc.

References Axis2, Reconstruct(), RejectISR(), rejectISR, and status.

Referenced by ~Hemisphere().

54  {
55  if (status != 1) {
56  if (rejectISR == 0) {
57  this->Reconstruct();
58  } else {
59  this->RejectISR();
60  }
61  }
62  return Axis2;
63 }
heppy::Hemisphere::Axis2
std::vector< float > Axis2
Definition: Hemisphere.h:184
vector< int > heppy::Hemisphere::getGrouping ( )

Definition at line 65 of file Hemisphere.cc.

References Object_Group, Reconstruct(), RejectISR(), rejectISR, and status.

Referenced by ~Hemisphere().

65  {
66  if (status != 1) {
67  if (rejectISR == 0) {
68  this->Reconstruct();
69  } else {
70  this->RejectISR();
71  }
72  }
73  return Object_Group;
74 }
heppy::Hemisphere::Object_Group
std::vector< int > Object_Group
Definition: Hemisphere.h:179
int heppy::Hemisphere::GetNumLoop ( )
inline

Definition at line 161 of file Hemisphere.h.

References numLoop, Reconstruct(), and RejectISR().

161 {return numLoop;}
int heppy::Hemisphere::Reconstruct ( )
private

Definition at line 76 of file Hemisphere.cc.

References Axis1, Axis2, gather_cfg::cout, dbg, MuonCkfTrajectoryBuilder_cfi::deltaPhi, deltaR(), dRminSeed1, hemi_meth, mps_fire::i, createfilelist::int, reco::tau::disc::InvariantMass(), cmsBatch::log, nItermax, numLoop, Object_E, Object_Eta, Object_Group, Object_Noassoc, Object_Noseed, Object_P, Object_Phi, Object_Pt, Object_Px, Object_Py, Object_Pz, seed_meth, mathSSE::sqrt(), status, funct::tan(), and theta().

Referenced by getAxis1(), getAxis2(), getGrouping(), GetNumLoop(), and RejectISR().

76  {
77 
78 // Performs the actual hemisphere reconstrucion
79 //
80 // definition of the vectors used internally:
81 // Object_xxx :
82 // xxx = Px, Py, Pz, E for input values
83 // xxx = P, Pt, Eta, Phi, Group for internal use
84 // Axis1 : final hemisphere axis 1
85 // Axis2 : final hemisphere axis 2
86 // Sum1_xxx : hemisphere 1 being updated during the association iterations
87 // Sum2_xxx : hemisphere 2 being updated during the association iterations
88 // NewAxis1_xxx, NewAxis1_xxx : temporary axes for calculation in association methods 2 and 3
89 
90 
91  numLoop=0; // initialize numLoop for Zero
92  int vsize = (int) Object_Px.size();
93  if((int) Object_Py.size() != vsize || (int) Object_Pz.size() != vsize){
94  cout << "WARNING!!!!! Input vectors have different size! Fix it!" << endl;
95  return 0;
96  }
97  if (dbg > 0) {
98 // cout << " Hemisphere method, vsn = " << hemivsn << endl;
99  cout << " Hemisphere method " << endl;
100  }
101 
102  // clear some vectors if method reconstruct() is called again
103  if (!Object_P.empty()){
104  Object_P.clear();
105  Object_Pt.clear();
106  Object_Eta.clear();
107  Object_Phi.clear();
108  Object_Group.clear();
109  Axis1.clear();
110  Axis2.clear();
111  }
112  // initialize the vectors
113  for(int j = 0; j < vsize; ++j){
114  Object_P.push_back(0);
115  Object_Pt.push_back(0);
116  Object_Eta.push_back(0);
117  Object_Phi.push_back(0);
118  Object_Group.push_back(0);
119  }
120  for(int j = 0; j < 5; ++j){
121  Axis1.push_back(0);
122  Axis2.push_back(0);
123  }
124 
125  // compute additional quantities for vectors Object_xxx
126  float theta;
127  for (int i = 0; i <vsize; ++i){
128  Object_P[i] = sqrt(Object_Px[i]*Object_Px[i]+Object_Py[i]*Object_Py[i]+Object_Pz[i]*Object_Pz[i]);
129  if (Object_P[i] > Object_E[i]+0.001) {
130  cout << "WARNING!!!!! Object " << i << " has E = " << Object_E[i]
131  << " less than P = " << Object_P[i] << " *** Fix it!" << endl;
132  return 0;
133  }
134  Object_Pt[i] = sqrt(Object_Px[i]*Object_Px[i]+Object_Py[i]*Object_Py[i]);
135  // protection for div by 0
136  if (fabs(Object_Pz[i]) > 0.001) {
137  theta = atan(sqrt(Object_Px[i]*Object_Px[i]+Object_Py[i]*Object_Py[i])/Object_Pz[i]);
138  }
139  else {
140  theta = 1.570796327;
141  }
142  if (theta < 0.) {theta = theta + 3.141592654;}
143  Object_Eta[i] = -log(tan(0.5*theta));
144  Object_Phi[i] = atan2(Object_Py[i], Object_Px[i]);
145  if (dbg > 0) {
146  cout << " Object " << i << " Eta = " << Object_Eta[i]
147  << " Phi = " << Object_Phi[i] << endl;
148  }
149  }
150 
151  if (dbg > 0) {
152  cout << endl;
153  cout << " Seeding method = " << seed_meth << endl;
154  }
155  // I_Max and J_Max are indices of the seeds in the vectors
156  int I_Max = -1;
157  int J_Max = -1;
158 
159  // determine the seeds for seed method 1
160  if (seed_meth == 1) {
161 
162  float P_Max = 0.;
163  float DeltaRP_Max = 0.;
164 
165  // take highest momentum object as first seed
166  for (int i = 0; i < vsize; ++i){
167  Object_Group[i] = 0;
168  //cout << "Object_Px[i] = " << Object_Px[i] << ", Object_Py[i] = " << Object_Py[i]
169  //<< ", Object_Pz[i] = " << Object_Pz[i] << " << endl;
170  if (Object_Noseed[i] == 0 && P_Max < Object_P[i]){
171  P_Max = Object_P[i];
172  I_Max = i;
173  }
174  }
175  // if 1st seed is found, save it as initial hemisphere 1 axis
176  if (I_Max >= 0){
177  Axis1[0] = Object_Px[I_Max] / Object_P[I_Max];
178  Axis1[1] = Object_Py[I_Max] / Object_P[I_Max];
179  Axis1[2] = Object_Pz[I_Max] / Object_P[I_Max];
180  Axis1[3] = Object_P[I_Max];
181  Axis1[4] = Object_E[I_Max];
182  } else {
183  // cout << " This is an empty event." << endl;
184  return 0;
185  }
186 
187  // take as second seed the object with largest DR*P w.r.t. the first seed
188  for (int i = 0; i < vsize; ++i){
189  // float DeltaR = sqrt((Object_Eta[i] - Object_Eta[I_Max])*(Object_Eta[i] - Object_Eta[I_Max])
190  // + (Util::DeltaPhi(Object_Phi[i], Object_Phi[I_Max]))*(Util::DeltaPhi(Object_Phi[i], Object_Phi[I_Max])) );
191  float DeltaR = sqrt((Object_Eta[i] - Object_Eta[I_Max])*(Object_Eta[i] - Object_Eta[I_Max])
192  + (deltaPhi(Object_Phi[i], Object_Phi[I_Max]))*(deltaPhi(Object_Phi[i], Object_Phi[I_Max])) );
193  if (Object_Noseed[i] == 0 && DeltaR > dRminSeed1) {
194  float DeltaRP = DeltaR * Object_P[i];
195  if (DeltaRP > DeltaRP_Max){
196  DeltaRP_Max = DeltaRP;
197  J_Max = i;
198  }
199  }
200  }
201  // if 2nd seed is found, save it as initial hemisphere 2 axis
202  if (J_Max >= 0){
203  Axis2[0] = Object_Px[J_Max] / Object_P[J_Max];
204  Axis2[1] = Object_Py[J_Max] / Object_P[J_Max];
205  Axis2[2] = Object_Pz[J_Max] / Object_P[J_Max];
206  Axis2[3] = Object_P[J_Max];
207  Axis2[4] = Object_E[J_Max];
208  } else {
209  // cout << " This is a MONOJET." << endl;
210  return 0;
211  }
212  if (dbg > 0) {
213  cout << " Axis 1 is Object = " << I_Max << endl;
214  cout << " Axis 2 is Object = " << J_Max << endl;
215  }
216 
217  // determine the seeds for seed methods 2 and 3
218  } else if (seed_meth == 2 || seed_meth == 3) {
219 
220  float Mass_Max = 0.;
221  float InvariantMass = 0.;
222 
223  // maximize the invariant mass of two objects
224  for (int i = 0; i < vsize; ++i){
225  Object_Group[i] = 0;
226  if (Object_Noseed[i] == 0){
227  for (int j = i+1; j < vsize; ++j){
228  if (Object_Noseed[j] == 0){
229  // either the invariant mass
230  if (seed_meth == 2){
231  InvariantMass = (Object_E[i] + Object_E[j])* (Object_E[i] + Object_E[j])
232  - (Object_Px[i] + Object_Px[j])* (Object_Px[i] + Object_Px[j])
233  - (Object_Py[i] + Object_Py[j])* (Object_Py[i] + Object_Py[j])
234  - (Object_Pz[i] + Object_Pz[j])* (Object_Pz[i] + Object_Pz[j]) ;
235  }
236  // or the transverse mass
237  else if (seed_meth == 3){
238  float pti = sqrt(Object_Px[i]*Object_Px[i] + Object_Py[i]*Object_Py[i]);
239  float ptj = sqrt(Object_Px[j]*Object_Px[j] + Object_Py[j]*Object_Py[j]);
240  InvariantMass = 2. * (pti*ptj - Object_Px[i]*Object_Px[j]
241  - Object_Py[i]*Object_Py[j] );
242  }
243  if (Mass_Max < InvariantMass){
244  Mass_Max = InvariantMass;
245  I_Max = i;
246  J_Max = j;
247  }
248  }
249  }
250  }
251  }
252 
253  // if both seeds are found, save them as initial hemisphere axes
254  if (J_Max > 0) {
255  Axis1[0] = Object_Px[I_Max] / Object_P[I_Max];
256  Axis1[1] = Object_Py[I_Max] / Object_P[I_Max];
257  Axis1[2] = Object_Pz[I_Max] / Object_P[I_Max];
258  Axis1[3] = Object_P[I_Max];
259  Axis1[4] = Object_E[I_Max];
260 
261  Axis2[0] = Object_Px[J_Max] / Object_P[J_Max];
262  Axis2[1] = Object_Py[J_Max] / Object_P[J_Max];
263  Axis2[2] = Object_Pz[J_Max] / Object_P[J_Max];
264  Axis2[3] = Object_P[J_Max];
265  Axis2[4] = Object_E[J_Max];
266  } else {
267  // cout << " This is a MONOJET." << endl;
268  return 0;
269  }
270  if (dbg > 0) {
271  cout << " Axis 1 is Object = " << I_Max << endl;
272  cout << " Axis 2 is Object = " << J_Max << endl;
273  }
274 
275  } else if (seed_meth == 4) {
276 
277  float P_Max1 = 0.;
278  float P_Max2 = 0.;
279 
280  // take largest Pt object as first seed
281  for (int i = 0; i < vsize; ++i){
282  Object_Group[i] = 0;
283  if (Object_Noseed[i] == 0 && P_Max1 < Object_Pt[i]){
284  P_Max1 = Object_Pt[i];
285  I_Max = i;
286  }
287  }
288  if(I_Max < 0) return 0;
289 
290  // take second largest Pt object as second seed, but require dR(seed1, seed2) > dRminSeed1
291  for (int i = 0; i < vsize; ++i){
292  if( i == I_Max) continue;
293  // float DeltaR = Util::GetDeltaR(Object_Eta[i], Object_Eta[I_Max], Object_Phi[i], Object_Phi[I_Max]);
294  float DeltaR = deltaR(Object_Eta[i], Object_Eta[I_Max], Object_Phi[i], Object_Phi[I_Max]);
295  if (Object_Noseed[i] == 0 && P_Max2 < Object_Pt[i] && DeltaR > dRminSeed1){
296  P_Max2 = Object_Pt[i];
297  J_Max = i;
298  }
299  }
300  if(J_Max < 0) return 0;
301 
302  // save first seed as initial hemisphere 1 axis
303  if (I_Max >= 0){
304  Axis1[0] = Object_Px[I_Max] / Object_P[I_Max];
305  Axis1[1] = Object_Py[I_Max] / Object_P[I_Max];
306  Axis1[2] = Object_Pz[I_Max] / Object_P[I_Max];
307  Axis1[3] = Object_P[I_Max];
308  Axis1[4] = Object_E[I_Max];
309  }
310 
311  // save second seed as initial hemisphere 2 axis
312  if (J_Max >= 0){
313  Axis2[0] = Object_Px[J_Max] / Object_P[J_Max];
314  Axis2[1] = Object_Py[J_Max] / Object_P[J_Max];
315  Axis2[2] = Object_Pz[J_Max] / Object_P[J_Max];
316  Axis2[3] = Object_P[J_Max];
317  Axis2[4] = Object_E[J_Max];
318  }
319 
320  if (dbg > 0) {
321  cout << " Axis 1 is Object = " << I_Max << " with Pt " << Object_Pt[I_Max]<< endl;
322  cout << " Axis 2 is Object = " << J_Max << " with Pt " << Object_Pt[J_Max]<< endl;
323  }
324 
325  } else if ( !(seed_meth == 0 && (hemi_meth == 8 || hemi_meth ==9) ) ) {
326  cout << "Please give a valid seeding method!" << endl;
327  return 0;
328  }
329 
330 
331  // seeding done
332  // now do the hemisphere association
333 
334  if (dbg > 0) {
335  cout << endl;
336  cout << " Association method = " << hemi_meth << endl;
337  }
338 
339  bool I_Move = true;
340 
341  // iterate to associate all objects to hemispheres (methods 1 to 3 only)
342  // until no objects are moved from one to the other hemisphere
343  // or the maximum number of iterations is reached
344  while (I_Move && (numLoop < nItermax) && hemi_meth != 8 && hemi_meth !=9){
345 
346  I_Move = false;
347  numLoop++;
348  if (dbg > 0) {
349  cout << " Iteration = " << numLoop << endl;
350  }
351  if(numLoop == nItermax-1){
352  cout << " Hemishpere: warning - reaching max number of iterations " << endl;
353  }
354 
355  // initialize the current sums of Px, Py, Pz, E for the two hemispheres
356  float Sum1_Px = 0.;
357  float Sum1_Py = 0.;
358  float Sum1_Pz = 0.;
359  float Sum1_E = 0.;
360  float Sum2_Px = 0.;
361  float Sum2_Py = 0.;
362  float Sum2_Pz = 0.;
363  float Sum2_E = 0.;
364 
365  // associate the objects for method 1
366  if (hemi_meth == 1) {
367 
368  for (int i = 0; i < vsize; ++i){
369  if (Object_Noassoc[i] == 0){
370  float P_Long1 = Object_Px[i]*Axis1[0] + Object_Py[i]*Axis1[1] + Object_Pz[i]*Axis1[2];
371  float P_Long2 = Object_Px[i]*Axis2[0] + Object_Py[i]*Axis2[1] + Object_Pz[i]*Axis2[2];
372  if (P_Long1 >= P_Long2){
373  if (Object_Group[i] != 1){
374  I_Move = true;
375  }
376  Object_Group[i] = 1;
377  Sum1_Px += Object_Px[i];
378  Sum1_Py += Object_Py[i];
379  Sum1_Pz += Object_Pz[i];
380  Sum1_E += Object_E[i];
381  } else {
382  if (Object_Group[i] != 2){
383  I_Move = true;
384  }
385  Object_Group[i] = 2;
386  Sum2_Px += Object_Px[i];
387  Sum2_Py += Object_Py[i];
388  Sum2_Pz += Object_Pz[i];
389  Sum2_E += Object_E[i];
390  }
391  }
392  }
393 
394  // associate the objects for methods 2 and 3
395  } else if (hemi_meth == 2 || hemi_meth == 3){
396 
397  for (int i = 0; i < vsize; ++i){
398  // add the seeds to the sums, as they remain fixed
399  if (i == I_Max) {
400  Object_Group[i] = 1;
401  Sum1_Px += Object_Px[i];
402  Sum1_Py += Object_Py[i];
403  Sum1_Pz += Object_Pz[i];
404  Sum1_E += Object_E[i];
405  } else if (i == J_Max) {
406  Object_Group[i] = 2;
407  Sum2_Px += Object_Px[i];
408  Sum2_Py += Object_Py[i];
409  Sum2_Pz += Object_Pz[i];
410  Sum2_E += Object_E[i];
411 
412  // for the other objects
413  } else {
414  if (Object_Noassoc[i] == 0){
415 
416  // only 1 object maximum is moved in a given iteration
417  if(!I_Move){
418  // initialize the new hemispheres as the current ones
419  float NewAxis1_Px = Axis1[0] * Axis1[3];
420  float NewAxis1_Py = Axis1[1] * Axis1[3];
421  float NewAxis1_Pz = Axis1[2] * Axis1[3];
422  float NewAxis1_E = Axis1[4];
423  float NewAxis2_Px = Axis2[0] * Axis2[3];
424  float NewAxis2_Py = Axis2[1] * Axis2[3];
425  float NewAxis2_Pz = Axis2[2] * Axis2[3];
426  float NewAxis2_E = Axis2[4];
427  // subtract the object from its hemisphere
428  if (Object_Group[i] == 1){
429  NewAxis1_Px = NewAxis1_Px - Object_Px[i];
430  NewAxis1_Py = NewAxis1_Py - Object_Py[i];
431  NewAxis1_Pz = NewAxis1_Pz - Object_Pz[i];
432  NewAxis1_E = NewAxis1_E - Object_E[i];
433  } else if (Object_Group[i] == 2) {
434  NewAxis2_Px = NewAxis2_Px - Object_Px[i];
435  NewAxis2_Py = NewAxis2_Py - Object_Py[i];
436  NewAxis2_Pz = NewAxis2_Pz - Object_Pz[i];
437  NewAxis2_E = NewAxis2_E - Object_E[i];
438  }
439 
440  // compute the invariant mass squared with each hemisphere (method 2)
441  float mass1 = NewAxis1_E
442  - ((Object_Px[i]*NewAxis1_Px + Object_Py[i]*NewAxis1_Py
443  + Object_Pz[i]*NewAxis1_Pz) / Object_P[i]);
444  float mass2 = NewAxis2_E
445  - ((Object_Px[i]*NewAxis2_Px + Object_Py[i]*NewAxis2_Py
446  + Object_Pz[i]*NewAxis2_Pz) / Object_P[i]);
447  // or the Lund distance (method 3)
448  if (hemi_meth == 3) {
449  mass1 *= NewAxis1_E/((NewAxis1_E+Object_E[i])*(NewAxis1_E+Object_E[i]));
450  mass2 *= NewAxis2_E/((NewAxis2_E+Object_E[i])*(NewAxis2_E+Object_E[i]));
451  }
452  // and associate the object to the best hemisphere and add it to the sum
453  if (mass1 < mass2) {
454  //if (Object_Group[i] != 1){
455  if (Object_Group[i] != 1 && Object_Group[i] != 0){
456  I_Move = true;
457  }
458  Object_Group[i] = 1;
459  Sum1_Px += Object_Px[i];
460  Sum1_Py += Object_Py[i];
461  Sum1_Pz += Object_Pz[i];
462  Sum1_E += Object_E[i];
463  } else {
464  //if (Object_Group[i] != 2){
465  if (Object_Group[i] != 2 && Object_Group[i] != 0){
466  I_Move = true;
467  }
468  Object_Group[i] = 2;
469  Sum2_Px += Object_Px[i];
470  Sum2_Py += Object_Py[i];
471  Sum2_Pz += Object_Pz[i];
472  Sum2_E += Object_E[i];
473  }
474 
475  // but if a previous object was moved, add all other associated objects to the sum
476  } else {
477  if (Object_Group[i] == 1){
478  Sum1_Px += Object_Px[i];
479  Sum1_Py += Object_Py[i];
480  Sum1_Pz += Object_Pz[i];
481  Sum1_E += Object_E[i];
482  } else if (Object_Group[i] == 2){
483  Sum2_Px += Object_Px[i];
484  Sum2_Py += Object_Py[i];
485  Sum2_Pz += Object_Pz[i];
486  Sum2_E += Object_E[i];
487  }
488  }
489  }
490  } // end loop over objects, Sum1_ and Sum2_ are now the updated hemispheres
491  }
492 
493  } else {
494  cout << "Please give a valid hemisphere association method!" << endl;
495  return 0;
496  }
497 
498  // recomputing the axes for next iteration
499 
500  Axis1[3] = sqrt(Sum1_Px*Sum1_Px + Sum1_Py*Sum1_Py + Sum1_Pz*Sum1_Pz);
501  if (Axis1[3] < 0.0001) {
502  cout << "ZERO objects in group 1! " << endl;
503  } else {
504  Axis1[0] = Sum1_Px / Axis1[3];
505  Axis1[1] = Sum1_Py / Axis1[3];
506  Axis1[2] = Sum1_Pz / Axis1[3];
507  Axis1[4] = Sum1_E;
508  }
509  Axis2[3] = sqrt(Sum2_Px*Sum2_Px + Sum2_Py*Sum2_Py + Sum2_Pz*Sum2_Pz);
510  if (Axis2[3] < 0.0001) {
511  cout << " ZERO objects in group 2! " << endl;
512  } else {
513  Axis2[0] = Sum2_Px / Axis2[3];
514  Axis2[1] = Sum2_Py / Axis2[3];
515  Axis2[2] = Sum2_Pz / Axis2[3];
516  Axis2[4] = Sum2_E;
517  }
518 
519  if (dbg > 0) {
520  cout << " Grouping = ";
521  for (int i=0;i<vsize;i++){
522  cout << " " << Object_Group[i];
523  }
524  cout << endl;
525  }
526 
527  if (numLoop <= 1) I_Move = true;
528 
529  } // end of iteration
530 
531  // associate all objects to hemispheres for method 8
532  if (hemi_meth == 8 || hemi_meth == 9) {
533  float sumtot = 0.;
534  for (int i = 0; i < vsize; ++i){
535  if (Object_Noassoc[i] != 0) continue;
536  if(hemi_meth==8){ sumtot += Object_E[i]; }
537  else sumtot += Object_Pt[i];
538  }
539  float sumdiff = fabs(sumtot - 2*Object_E[0]);
540  float sum1strt = 0.;
541  int ibest = 0, jbest = 0;
542 
543  // start by choosing object 0 in hemi 1, all others in hemi 2
544  // then add each of the other objects one by one to hemi 1
545  // then add object 1 to hemi one and add each of the others in turn, etc
546  if(vsize > 2){
547  for (int i = 0; i < vsize-1; ++i){
548  if (Object_Noassoc[i] != 0) continue;
549  if (hemi_meth ==8 ){sum1strt += Object_E[i];}
550  else {sum1strt += Object_Pt[i];}
551  for (int j = i+1; j < vsize; ++j){
552  if (Object_Noassoc[j] != 0) continue;
553  float sum1_E =0;
554  if (hemi_meth ==8 ){
555  sum1_E = sum1strt + Object_E[j];
556  }else {
557  sum1_E = sum1strt + Object_Pt[j];
558  }
559  float sum2_E = sumtot - sum1_E;
560  if(sumdiff >= fabs(sum1_E - sum2_E)) {
561  sumdiff = fabs(sum1_E - sum2_E);
562  ibest = i;
563  jbest = j;
564  }
565  }
566  }
567  }
568 
569  // then store the best combination into the hemisphere axes
570  float Sum1_Px=0, Sum1_Py=0, Sum1_Pz=0, Sum1_E=0;
571  float Sum2_Px=0, Sum2_Py=0, Sum2_Pz=0, Sum2_E=0;
572  for (int i = 0; i < vsize; ++i){
573  if (Object_Noassoc[i] != 0) Object_Group[i] = 0;
574  else if (i <= ibest || i == jbest) {
575  Sum1_Px += Object_Px[i];
576  Sum1_Py += Object_Py[i];
577  Sum1_Pz += Object_Pz[i];
578  Sum1_E += Object_E[i];
579  Object_Group[i] = 1;
580  } else {
581  Sum2_Px += Object_Px[i];
582  Sum2_Py += Object_Py[i];
583  Sum2_Pz += Object_Pz[i];
584  Sum2_E += Object_E[i];
585  Object_Group[i] = 2;
586  }
587  }
588  Axis1[3] = sqrt(Sum1_Px*Sum1_Px + Sum1_Py*Sum1_Py + Sum1_Pz*Sum1_Pz);
589  Axis1[0] = Sum1_Px / Axis1[3];
590  Axis1[1] = Sum1_Py / Axis1[3];
591  Axis1[2] = Sum1_Pz / Axis1[3];
592  Axis1[4] = Sum1_E;
593  Axis2[3] = sqrt(Sum2_Px*Sum2_Px + Sum2_Py*Sum2_Py + Sum2_Pz*Sum2_Pz);
594  Axis2[0] = Sum2_Px / Axis2[3];
595  Axis2[1] = Sum2_Py / Axis2[3];
596  Axis2[2] = Sum2_Pz / Axis2[3];
597  Axis2[4] = Sum2_E;
598  }
599 
600  status = 1;
601  return 1;
602 }
DeltaR
Definition: DeltaR.py:1
mps_fire.i
i
Definition: mps_fire.py:156
heppy::Hemisphere::Object_Phi
std::vector< float > Object_Phi
Definition: Hemisphere.h:177
heppy::Hemisphere::Object_P
std::vector< float > Object_P
Definition: Hemisphere.h:174
theta
Geom::Theta< T > theta() const
Definition: Basic3DVectorLD.h:170
heppy::Hemisphere::Axis2
std::vector< float > Axis2
Definition: Hemisphere.h:184
heppy::Hemisphere::Object_Eta
std::vector< float > Object_Eta
Definition: Hemisphere.h:178
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
heppy::Hemisphere::Object_Py
std::vector< float > Object_Py
Definition: Hemisphere.h:172
heppy::Hemisphere::Object_Noseed
std::vector< int > Object_Noseed
Definition: Hemisphere.h:180
funct::tan
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
heppy::Hemisphere::Object_E
std::vector< float > Object_E
Definition: Hemisphere.h:176
deltaR
double deltaR(double eta1, double eta2, double phi1, double phi2)
Definition: TreeUtility.cc:17
heppy::Hemisphere::Object_Pt
std::vector< float > Object_Pt
Definition: Hemisphere.h:175
heppy::Hemisphere::Object_Noassoc
std::vector< int > Object_Noassoc
Definition: Hemisphere.h:181
heppy::Hemisphere::Object_Pz
std::vector< float > Object_Pz
Definition: Hemisphere.h:173
heppy::Hemisphere::Object_Px
std::vector< float > Object_Px
Definition: Hemisphere.h:171
MuonCkfTrajectoryBuilder_cfi.deltaPhi
deltaPhi
Definition: MuonCkfTrajectoryBuilder_cfi.py:27
gather_cfg.cout
cout
Definition: gather_cfg.py:145
heppy::Hemisphere::Axis1
std::vector< float > Axis1
Definition: Hemisphere.h:183
heppy::Hemisphere::Object_Group
std::vector< int > Object_Group
Definition: Hemisphere.h:179
int heppy::Hemisphere::RejectISR ( )
private

Definition at line 604 of file Hemisphere.cc.

References Axis1, Axis2, MuonCkfTrajectoryBuilder_cfi::deltaPhi, HiRegitMuonDetachedTripletStep_cff::DeltaR, mps_fire::i, createfilelist::int, cmsBatch::log, Object_Eta, Object_Group, Object_Phi, Object_Px, Object_Py, Object_Pz, Reconstruct(), rejectISRDR, rejectISRDRmax, rejectISRPt, rejectISRPtmax, SetNoAssoc(), mathSSE::sqrt(), status, funct::tan(), and theta().

Referenced by getAxis1(), getAxis2(), getGrouping(), and GetNumLoop().

604  {
605 // tries to avoid including ISR jets into hemispheres
606 //
607 
608  // iterate to remove all ISR objects from the hemispheres
609  // until no ISR objects are found
610 // cout << " entered RejectISR() with rejectISRDR = " << rejectISRDR << endl;
611  bool I_Move = true;
612  while (I_Move) {
613  I_Move = false;
614  float valmax = 0.;
615  int imax = -1;
616 
617  // start by making a hemisphere reconstruction
618  int hemiOK = Reconstruct();
619  if (hemiOK == 0) {return 0;}
620 
621  // convert the axes into Px, Py, Pz, E vectors
622  float newAxis1_Px = Axis1[0] * Axis1[3];
623  float newAxis1_Py = Axis1[1] * Axis1[3];
624  float newAxis1_Pz = Axis1[2] * Axis1[3];
625  //float newAxis1_E = Axis1[4];
626  float newAxis2_Px = Axis2[0] * Axis2[3];
627  float newAxis2_Py = Axis2[1] * Axis2[3];
628  float newAxis2_Pz = Axis2[2] * Axis2[3];
629  //float newAxis2_E = Axis2[4];
630 
631  // loop over all objects associated to a hemisphere
632  int vsize = (int) Object_Px.size();
633  for (int i = 0; i < vsize; ++i){
634  if (Object_Group[i] == 1 || Object_Group[i] == 2){
635 // cout << " Object = " << i << ", Object_Group = " << Object_Group[i] << endl;
636 
637  // collect the hemisphere data
638  float newPx = 0.;
639  float newPy = 0.;
640  float newPz = 0.;
641  //float newE = 0.;
642  if (Object_Group[i] == 1){
643  newPx = newAxis1_Px;
644  newPy = newAxis1_Py;
645  newPz = newAxis1_Pz;
646  //newE = newAxis1_E;
647  } else if (Object_Group[i] == 2) {
648  newPx = newAxis2_Px;
649  newPy = newAxis2_Py;
650  newPz = newAxis2_Pz;
651  //newE = newAxis2_E;
652  }
653 
654  // compute the quantities to test whether the object is ISR
655  float ptHemi = 0.;
656  float hemiEta = 0.;
657  float hemiPhi = 0.;
658  if (rejectISRPt == 1) {
659  float plHemi = (Object_Px[i]*newPx + Object_Py[i]*newPy
660  + Object_Pz[i]*newPz)
661  / sqrt(newPx*newPx+newPy*newPy+newPz*newPz);
662  float pObjsq = Object_Px[i]*Object_Px[i] + Object_Py[i]*Object_Py[i]
663  + Object_Pz[i]*Object_Pz[i];
664  ptHemi = sqrt(pObjsq - plHemi*plHemi);
665  if (ptHemi > valmax) {
666  valmax = ptHemi;
667  imax = i;
668  }
669  } else if (rejectISRDR == 1) {
670  float theta = 1.570796327;
671  // compute the new hemisphere eta, phi and DeltaR
672  float pdiff = fabs(newPx-Object_Px[i]) + fabs(newPy-Object_Py[i])
673  + fabs(newPz-Object_Pz[i]);
674  if (pdiff > 0.001) {
675  if (fabs(newPz) > 0.001) {
676  theta = atan(sqrt(newPx*newPx+newPy*newPy)/newPz);
677  }
678  if (theta < 0.) {theta = theta + 3.141592654;}
679  hemiEta = -log(tan(0.5*theta));
680  hemiPhi = atan2(newPy, newPx);
681  // float DeltaR = sqrt((Object_Eta[i] - hemiEta)*(Object_Eta[i] - hemiEta)
682  // + (Util::DeltaPhi(Object_Phi[i], hemiPhi))*(Util::DeltaPhi(Object_Phi[i], hemiPhi)) );
683  float DeltaR = sqrt((Object_Eta[i] - hemiEta)*(Object_Eta[i] - hemiEta)
684  + (deltaPhi(Object_Phi[i], hemiPhi))*(deltaPhi(Object_Phi[i], hemiPhi)) );
685 // cout << " Object = " << i << ", DeltaR = " << DeltaR << endl;
686  if (DeltaR > valmax) {
687  valmax = DeltaR;
688  imax = i;
689  }
690  }
691  }
692  }
693  } // end loop over objects
694 
695  // verify whether the ISR tests are fulfilled
696  if (imax < 0) {
697  I_Move = false;
698  } else if (rejectISRPt == 1 && valmax > rejectISRPtmax) {
699  SetNoAssoc(imax);
700  I_Move = true;
701  } else if (rejectISRDR == 1 && valmax > rejectISRDRmax) {
702  SetNoAssoc(imax);
703  I_Move = true;
704  }
705 
706  } // end iteration over ISR objects
707 
708  status = 1;
709  return 1;
710 }
DeltaR
Definition: DeltaR.py:1
mps_fire.i
i
Definition: mps_fire.py:156
heppy::Hemisphere::Object_Phi
std::vector< float > Object_Phi
Definition: Hemisphere.h:177
theta
Geom::Theta< T > theta() const
Definition: Basic3DVectorLD.h:170
heppy::Hemisphere::Axis2
std::vector< float > Axis2
Definition: Hemisphere.h:184
heppy::Hemisphere::Object_Eta
std::vector< float > Object_Eta
Definition: Hemisphere.h:178
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
heppy::Hemisphere::Object_Py
std::vector< float > Object_Py
Definition: Hemisphere.h:172
HiRegitMuonDetachedTripletStep_cff.DeltaR
DeltaR
Definition: HiRegitMuonDetachedTripletStep_cff.py:22
funct::tan
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
heppy::Hemisphere::Object_Pz
std::vector< float > Object_Pz
Definition: Hemisphere.h:173
heppy::Hemisphere::Object_Px
std::vector< float > Object_Px
Definition: Hemisphere.h:171
MuonCkfTrajectoryBuilder_cfi.deltaPhi
deltaPhi
Definition: MuonCkfTrajectoryBuilder_cfi.py:27
heppy::Hemisphere::Axis1
std::vector< float > Axis1
Definition: Hemisphere.h:183
heppy::Hemisphere::Object_Group
std::vector< int > Object_Group
Definition: Hemisphere.h:179
heppy::Hemisphere::SetNoAssoc
void SetNoAssoc(int object_number)
Definition: Hemisphere.h:113
void heppy::Hemisphere::RejectISRDRmax ( float  drmax)
inline

Definition at line 151 of file Hemisphere.h.

References rejectISR, rejectISRDR, rejectISRDRmax, rejectISRPt, and status.

151  {
152  rejectISRDRmax = drmax;
153  rejectISRDR = 1;
154  rejectISRPt = 0;
155  rejectISR = 1;
156  status = 0;
157  }
void heppy::Hemisphere::RejectISRPtmax ( float  ptmax)
inline

Definition at line 142 of file Hemisphere.h.

References rejectISR, rejectISRDR, rejectISRPt, rejectISRPtmax, and status.

142  {
143  rejectISRPtmax = ptmax;
144  rejectISRPt = 1;
145  rejectISRDR = 0;
146  rejectISR = 1;
147  status = 0;
148  }
void heppy::Hemisphere::SetDebug ( int  debug)
inline

Definition at line 160 of file Hemisphere.h.

References dbg, and debug.

160 { dbg = debug; }
debug
#define debug
Definition: HDRShower.cc:19
void heppy::Hemisphere::SetDRminSeed1 ( float  rmin)
inline

Definition at line 131 of file Hemisphere.h.

References dRminSeed1.

131  {
132  dRminSeed1 = rmin;
133  }
void heppy::Hemisphere::SetMethod ( int  seed_method,
int  hemisphere_association_method 
)
inline

Definition at line 98 of file Hemisphere.h.

References hemi_meth, seed_meth, and status.

98  {
99  seed_meth = seed_method;
100  hemi_meth = hemisphere_association_method;
101  status = 0;
102  }
void heppy::Hemisphere::SetnItermax ( int  niter)
inline

Definition at line 137 of file Hemisphere.h.

References nItermax.

137  {
138  nItermax = niter;
139  }
void heppy::Hemisphere::SetNoAssoc ( int  object_number)
inline

Definition at line 113 of file Hemisphere.h.

References Object_Noassoc, Object_Noseed, and status.

Referenced by RejectISR().

113  {
114  Object_Noassoc[object_number] = 1;
115  Object_Noseed[object_number] = 1;
116  status = 0;
117  }
heppy::Hemisphere::Object_Noseed
std::vector< int > Object_Noseed
Definition: Hemisphere.h:180
heppy::Hemisphere::Object_Noassoc
std::vector< int > Object_Noassoc
Definition: Hemisphere.h:181
void heppy::Hemisphere::SetNoSeed ( int  object_number)
inline

Definition at line 106 of file Hemisphere.h.

References Object_Noseed, and status.

106  {
107  Object_Noseed[object_number] = 1;
108  status = 0;
109  }
heppy::Hemisphere::Object_Noseed
std::vector< int > Object_Noseed
Definition: Hemisphere.h:180

Member Data Documentation

std::vector<float> heppy::Hemisphere::Axis1
private

Definition at line 183 of file Hemisphere.h.

Referenced by getAxis1(), Reconstruct(), and RejectISR().

std::vector<float> heppy::Hemisphere::Axis2
private

Definition at line 184 of file Hemisphere.h.

Referenced by getAxis2(), Reconstruct(), and RejectISR().

int heppy::Hemisphere::dbg
private

Definition at line 197 of file Hemisphere.h.

Referenced by Reconstruct(), and SetDebug().

float heppy::Hemisphere::dRminSeed1
private

Definition at line 190 of file Hemisphere.h.

Referenced by Reconstruct(), and SetDRminSeed1().

int heppy::Hemisphere::hemi_meth
private

Definition at line 188 of file Hemisphere.h.

Referenced by Reconstruct(), and SetMethod().

int heppy::Hemisphere::nItermax
private

Definition at line 191 of file Hemisphere.h.

Referenced by Reconstruct(), and SetnItermax().

int heppy::Hemisphere::numLoop
private

Definition at line 198 of file Hemisphere.h.

Referenced by GetNumLoop(), Hemisphere(), and Reconstruct().

std::vector<float> heppy::Hemisphere::Object_E
private

Definition at line 176 of file Hemisphere.h.

Referenced by Reconstruct().

std::vector<float> heppy::Hemisphere::Object_Eta
private

Definition at line 178 of file Hemisphere.h.

Referenced by Reconstruct(), and RejectISR().

std::vector<int> heppy::Hemisphere::Object_Group
private

Definition at line 179 of file Hemisphere.h.

Referenced by getGrouping(), Reconstruct(), and RejectISR().

std::vector<int> heppy::Hemisphere::Object_Noassoc
private

Definition at line 181 of file Hemisphere.h.

Referenced by ClearAllNoLists(), Hemisphere(), Reconstruct(), and SetNoAssoc().

std::vector<int> heppy::Hemisphere::Object_Noseed
private

Definition at line 180 of file Hemisphere.h.

Referenced by ClearAllNoLists(), Hemisphere(), Reconstruct(), SetNoAssoc(), and SetNoSeed().

std::vector<float> heppy::Hemisphere::Object_P
private

Definition at line 174 of file Hemisphere.h.

Referenced by Reconstruct().

std::vector<float> heppy::Hemisphere::Object_Phi
private

Definition at line 177 of file Hemisphere.h.

Referenced by Reconstruct(), and RejectISR().

std::vector<float> heppy::Hemisphere::Object_Pt
private

Definition at line 175 of file Hemisphere.h.

Referenced by Reconstruct().

std::vector<float> heppy::Hemisphere::Object_Px
private

Definition at line 171 of file Hemisphere.h.

Referenced by Hemisphere(), Reconstruct(), and RejectISR().

std::vector<float> heppy::Hemisphere::Object_Py
private

Definition at line 172 of file Hemisphere.h.

Referenced by Reconstruct(), and RejectISR().

std::vector<float> heppy::Hemisphere::Object_Pz
private

Definition at line 173 of file Hemisphere.h.

Referenced by Reconstruct(), and RejectISR().

int heppy::Hemisphere::rejectISR
private

Definition at line 192 of file Hemisphere.h.

Referenced by getAxis1(), getAxis2(), getGrouping(), RejectISRDRmax(), and RejectISRPtmax().

int heppy::Hemisphere::rejectISRDR
private

Definition at line 195 of file Hemisphere.h.

Referenced by RejectISR(), RejectISRDRmax(), and RejectISRPtmax().

float heppy::Hemisphere::rejectISRDRmax
private

Definition at line 196 of file Hemisphere.h.

Referenced by RejectISR(), and RejectISRDRmax().

int heppy::Hemisphere::rejectISRPt
private

Definition at line 193 of file Hemisphere.h.

Referenced by RejectISR(), RejectISRDRmax(), and RejectISRPtmax().

float heppy::Hemisphere::rejectISRPtmax
private

Definition at line 194 of file Hemisphere.h.

Referenced by RejectISR(), and RejectISRPtmax().

int heppy::Hemisphere::seed_meth
private

Definition at line 187 of file Hemisphere.h.

Referenced by Reconstruct(), and SetMethod().

int heppy::Hemisphere::status
private

Definition at line 189 of file Hemisphere.h.

Referenced by dirstructure.Comparison::__make_image(), dirstructure.Comparison::__repr__(), ClearAllNoLists(), getAxis1(), getAxis2(), getGrouping(), Reconstruct(), RejectISR(), RejectISRDRmax(), RejectISRPtmax(), SetMethod(), SetNoAssoc(), and SetNoSeed().

Generated for CMSSW Reference Manual by doxygen 1.8.11