FastJet  3.1.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
ConeClusterAlgo.hpp
1 //////////////////////////////////////////////////////////////
2 // File: ConeClusterAlgo.hpp
3 //
4 // Author: G. Le Meur & F. Touze
5 //
6 // Created: 15-JUNE-1998
7 //
8 // Purpose: make jet clusters using fixed cone like algorithm
9 // implemented in RUNI.
10 //
11 // Modified:
12 // 28-OCT-1998 to use KinemUtil (S. Protopopescu)
13 // 8-JAN-1999: Laurent Duflot
14 // . correct bugs in getItemsInCone and updateEtaPhiEt for jets
15 // overlapping the phi=0 line
16 // . change abs(float) to fabs(float)
17 // 1-NOV-1999: Laurent Duflot
18 // . correct bug in makeCluster: when the temporary jet was emptied the eta
19 // and phi were not set again. The main effect was a nearly zero
20 // efficiency for jets at phi=pi (as seen by Volker Buescher)
21 // 25-JAN-2000: Francois Touze
22 // . change in updateEtaPhiEt : the method E() which returns energy doesn't
23 // exist in MCparticle classe,... so use the 4-momentum components
24 // . declare const the EnergyClusterCollection of seeds in makeClusters
25 // 01-FEB-2000: Laurent Duflot
26 // . add a missing break statement in the removal of shared items. Caused
27 // an infinite loop on some events.
28 // . correct typo in variable name. Change a variable name that was in
29 // French.
30 // . leave some debug printout (commented)
31 // 15-Sep-2009 Lars Sonnenschein
32 // extracted from D0 software framework and modified to remove subsequent dependencies
33 //
34 //
35 // This file is distributed with FastJet under the terms of the GNU
36 // General Public License (v2). Permission to do so has been granted
37 // by Lars Sonnenschein and the D0 collaboration (see COPYING for
38 // details)
39 //
40 // History of changes in FastJet compared to the original version of
41 // ConeClusterAlgo.hpp
42 //
43 // 2011-12-13 Gregory Soyez <soyez@fastjet.fr>
44 //
45 // * added license information
46 //
47 // 2011-10-06 Gregory Soyez <soyez@fastjet.fr>
48 //
49 // * put the code in the fastjet::d0runi namespace
50 //
51 //////////////////////////////////////////////////////////////
52 
53 //#ifndef CONECLUSTERALGO_H
54 //#define CONECLUSTERALGO_H
55 
56 #ifndef D0RunIconeJets_CONECLUSTERALGO_H
57 #define D0RunIconeJets_CONECLUSTERALGO_H
58 
59 
60 //#include "EnergyClusterReco.hpp"
61 #include <vector>
62 #include <list>
63 #include <utility>
64 //#include "kinem_util/AnglesUtil.hpp"
65 
66 #include <algorithm>
67 #include <iostream>
68 
69 #include "inline_maths.h"
70 #include <fastjet/internal/base.hh>
71 
72 FASTJET_BEGIN_NAMESPACE
73 
74 namespace d0runi{
75 
76 using namespace std;
77 
78 //some utility functions
79 inline float R2(float eta1, float phi1, float eta2, float phi2) {
80  return (eta1-eta2)*(eta1-eta2)+(phi1-phi2)*(phi1-phi2); }
81 
82 inline float R2_bis(float eta1, float phi1, float eta2, float phi2) {
83  //float dphi = kinem::delta_phi(phi1,phi2);
84  float dphi = inline_maths::delta_phi(phi1,phi2);
85  return (eta1-eta2)*(eta1-eta2)+dphi*dphi; }
86 
87 inline float DELTA_r(float eta1,float eta2,float phi1,float phi2) {
88  //float dphi = kinem::delta_phi(phi1,phi2);
89  float dphi = inline_maths::delta_phi(phi1,phi2);
90  return sqrt((eta1-eta2)*(eta1-eta2)+dphi*dphi);
91 }
92 
93 inline float E2eta(float* p) {
94  float small= 1.E-05;
95  float E[3];
96  if(p[3] < 0.0) {
97  E[0]= -p[0];
98  E[1]= -p[1];
99  E[2]= -p[2];
100  }
101  else {
102  E[0]= p[0];
103  E[1]= p[1];
104  E[2]= p[2];
105  }
106  float pperp= sqrt(E[0]*E[0]+E[1]*E[1])+small;
107  float ptotal= sqrt(E[0]*E[0]+E[1]*E[1]+E[2]*E[2])+small;
108  //float theta= atan2(pperp,E[2]);
109 
110  float eta= 0.0;
111  if(E[2] > 0.0) eta= log((ptotal+E[2])/pperp);
112  else eta= log(pperp/(ptotal-E[2]));
113  return eta;
114 }
115 
116 inline float E2phi(float* p) {
117  float small= 1.E-05;
118  float E[3];
119  if(p[3] < 0.0) {
120  E[0]= -p[0];
121  E[1]= -p[1];
122  E[2]= -p[2];
123  }
124  else {
125  E[0]= p[0];
126  E[1]= p[1];
127  E[2]= p[2];
128  }
129  float phi= atan2(E[1],E[0]+small);
130  //if(phi < 0.0) phi+=kinem::TWOPI;
131  if (phi < 0.0) phi += inline_maths::TWOPI;
132  return phi;
133 }
134 
135 //template < class CalItem,class CalItemAddress,class CalIClusterChunk >
136 template < class CalItem >
137 class ConeClusterAlgo {
138  //
139  // Purpose: make calorimeter clusters using a cone algorithm from
140  // preclusters created previously by the class ConePreClusterAlgo.
141  // Items must have addresses and 4-momenta.
142  // The algorithm is implemented with a template function makeClusters.
143  //
144  public :
145 
146 
147 //default constructor
148 ConeClusterAlgo() {}
149 
150 //constructor for cone jet algorithm
151 ConeClusterAlgo( float CONErad,float JETmne,float SPLifr):
152  _CONErad(fabs(CONErad)),
153  _JETmne(JETmne),
154  _SPLifr(SPLifr),
155  _TWOrad(0.),
156  _D0_Angle(false),
157  _Increase_Delta_R(true),
158  _Kill_Far_Clusters(true),
159  _Jet_Et_Min_On_Iter(true),
160  _Far_Ratio(0.5),
161  _Eitem_Negdrop(-1.0),
162  _Et_Min_Ratio(0.5),
163  _Thresh_Diff_Et(0.01)
164  {}
165 
166 //changing default thresholds & parameters
167 // (declared by PARAMETER in RUNI code)
168 ConeClusterAlgo( float CONErad,float JETmne,float SPLifr,float TWOrad,
169  float Tresh_Diff_Et,bool D0_Angle,bool Increase_Delta_R,
170  bool Kill_Far_Clusters,bool Jet_Et_Min_On_Iter,
171  float Far_Ratio,float Eitem_Negdrop,float Et_Min_Ratio ):
172  _CONErad(fabs(CONErad)),
173  _JETmne(JETmne),
174  _SPLifr(SPLifr),
175  _TWOrad(TWOrad),
176  _D0_Angle(D0_Angle),
177  _Increase_Delta_R(Increase_Delta_R),
178  _Kill_Far_Clusters(Kill_Far_Clusters),
179  _Jet_Et_Min_On_Iter(Jet_Et_Min_On_Iter),
180  _Far_Ratio(Far_Ratio),
181  _Eitem_Negdrop(Eitem_Negdrop),
182  _Et_Min_Ratio(Et_Min_Ratio),
183  _Thresh_Diff_Et(Tresh_Diff_Et)
184  {}
185 
186 //destructor
187 ~ConeClusterAlgo() {}
188 
189 //to make jet clusters using cone algorithm
190 void makeClusters(//const EnergyClusterReco* r,
191  std::list<CalItem> &jets,
192  list<const CalItem*> &itemlist, float Zvertex
193  //, const EnergyClusterCollection<CalItemAddress> &preclu,
194  //CalIClusterChunk* chunkptr
195  //) const;
196  );
197 
198 //print parameters of the algorithm
199 void print(ostream &os)const;
200 
201  //vector< TemporaryJet > TempColl;
202 
203 
204  private :
205 
206  float _CONErad;
207  float _JETmne;
208  float _SPLifr;
209 
210  float _TWOrad;
211  bool _D0_Angle;
212  bool _Increase_Delta_R;
213  bool _Kill_Far_Clusters;
214  bool _Jet_Et_Min_On_Iter;
215  float _Far_Ratio;
216  float _Eitem_Negdrop;
217  float _Et_Min_Ratio;
218  float _Thresh_Diff_Et;
219 
220  class TemporaryJet {
221 
222  public:
223 
224 
225 
226  TemporaryJet() {}
227 
228  TemporaryJet(float eta,float phi) {
229  _Eta=eta;
230  _Phi=phi;
231  }
232 
233  ~TemporaryJet() {}
234 
235  void addItem(const CalItem* tw) {
236  _LItems.push_back(tw);
237  }
238 
239  void setEtaPhiEt(float eta,float phi,float pT) {
240  _Eta= eta;
241  _Phi= phi;
242  _Et = pT;
243  }
244 
245  void erase() {
246  _LItems.erase(_LItems.begin(),_LItems.end());
247  _Eta= 0.0;
248  _Phi= 0.0;
249  _Et = 0.0;
250  }
251 
252  bool share_jets(TemporaryJet &NewJet,float SharedFr,float SPLifr) {
253  //
254  // combined
255  //
256  if(SharedFr >= SPLifr) {
257  typename list<const CalItem*>::iterator it;
258  typename list<const CalItem*>::iterator end_of_old=_LItems.end();
259  for(it=NewJet._LItems.begin(); it!=NewJet._LItems.end(); it++) {
260  typename list<const CalItem*>::iterator
261  where = find(_LItems.begin(),end_of_old,*it);
262  // if the item is not shared, add to this jet
263  if(where == end_of_old) {
264  _LItems.push_back(*it);
265  }
266  }
267  NewJet.erase();
268  return false;
269  } else {
270  //
271  // split
272  //
273  typename list<const CalItem*>::iterator it;
274  for(it=NewJet._LItems.begin(); it!=NewJet._LItems.end(); ) {
275  typename list<const CalItem*>::iterator
276  where = find(_LItems.begin(),_LItems.end(),*it);
277  if(where != _LItems.end()) {
278  //float EtaItem=(*it)->eta();
279  //float PhiItem=(*it)->phi();
280  // stay closer to the RUNI conventions for negative E cells
281  float pz[4];
282  (*it)->p4vec(pz);
283  float EtaItem= E2eta(pz);
284  float PhiItem= E2phi(pz);
285 
286  float RadOld=R2_bis(_Eta,_Phi,EtaItem,PhiItem);
287  float RadNew=R2_bis(NewJet.Eta(),NewJet.Phi(),EtaItem,PhiItem);
288  if (RadNew > RadOld) {
289  it = NewJet._LItems.erase(it);
290  }
291  else {
292  _LItems.erase(where);
293  ++it;
294  }
295  }
296  else ++it;
297  }
298  return true;
299  }
300  }
301 
302 
303  float dist_R2(TemporaryJet &jet) const {
304  float deta= _Eta-jet.Eta();
305  //float dphi= kinem::delta_phi(_Phi,jet.Phi());
306  float dphi= inline_maths::delta_phi(_Phi,jet.Phi());
307  return (deta*deta+dphi*dphi);
308  }
309 
310  bool ItemInJet(const CalItem* tw) const {
311  typename list<const CalItem*>::const_iterator
312  where= find(_LItems.begin(),_LItems.end(),tw);
313  if(where != _LItems.end()) return true;
314  else return false;
315  }
316 
317  bool updateEtaPhiEt() {
318  float ETsum = 0.0;
319  float ETAsum= 0.0;
320  float PHIsum= 0.0;
321  float Esum= 0.0;
322  typename list<const CalItem*>::iterator it;
323  for(it=_LItems.begin(); it!=_LItems.end(); it++) {
324  float ETk = (*it)->pT();
325  // now done in CalCell/CalTower if((*it)->E() < 0.0) ETk= -ETk;
326 
327  //float ETAk= (*it)->eta();
328  //float PHIk= (*it)->phi();
329  float pz[4];
330  (*it)->p4vec(pz);
331  float ETAk= E2eta(pz);
332  // take care of the phi=0=2pi problem
333  float PHIk= E2phi(pz);
334  //if(fabs(PHIk-_Phi) > kinem::TWOPI-fabs(PHIk-_Phi))
335  if(fabs(PHIk-_Phi) > inline_maths::TWOPI-fabs(PHIk-_Phi))
336  {
337  if(_Phi < PHIk)
338  {
339  //PHIk -= kinem::TWOPI;
340  PHIk -= inline_maths::TWOPI;
341  }
342  else
343  {
344  //PHIk += kinem::TWOPI;
345  PHIk += inline_maths::TWOPI;
346  }
347  }
348  ETAsum+= ETAk*ETk;
349  PHIsum+= PHIk*ETk;
350  ETsum += ETk;
351  // Esum+=(*it)->E(); use 4-momentum components
352  Esum+= pz[3];
353  }
354  if(ETsum <= 0.0) {
355  _Eta= 0.0;
356  _Phi= 0.0;
357  _Et = 0.0;
358  _E=0.;
359  return false;
360  }
361  else {
362  _Eta= ETAsum/ETsum;
363  _Phi= PHIsum/ETsum;
364  //if ( _Phi<0 ) _Phi+=kinem::TWOPI;
365  if ( _Phi<0 ) _Phi+=inline_maths::TWOPI;
366  _Et = ETsum;
367  _E = Esum;
368  return true;
369  }
370  }
371 
372  void D0_Angle_updateEtaPhi() {
373  float EXsum = 0.0;
374  float EYsum = 0.0;
375  float EZsum = 0.0;
376  typename list<const CalItem*>::iterator it;
377  for(it=_LItems.begin(); it!=_LItems.end(); it++) {
378  float p[4];
379  (*it)->p4vec(p);
380  EXsum += p[0];
381  EYsum += p[1];
382  EZsum += p[2];
383  }
384  //_Phi=kinem::phi(EYsum,EXsum);
385  _Phi=inline_maths::phi(EYsum,EXsum);
386  //_Eta=kinem::eta(EXsum,EYsum,EZsum);
387  _Eta=inline_maths::eta(EXsum,EYsum,EZsum);
388  }
389 
390  void getItems(list<const CalItem*> &ecv) const {
391  ecv.clear(); //ls 27/Feb/2010
392  typename list<const CalItem*>::const_iterator it;
393  for(it=_LItems.begin(); it!=_LItems.end(); it++) {
394  ecv.push_back(*it);
395  }
396  }
397 
398  float Eta() {return _Eta;}
399  float Phi() {return _Phi;}
400  float Et() {return _Et;}
401  float E() {return _E;}
402  list<const CalItem*> &LItems() {return _LItems;}
403 
404  private:
405  list<const CalItem*> _LItems;
406  float _Eta;
407  float _Phi;
408  float _Et;
409  float _E;
410  }; //class TemporaryJet
411 
412  void getItemsInCone(list<const CalItem*> &tlist, float etaJet, float phiJet,
413  float cone_radius, float zvertex_in) const;
414  void getItemsInCone_bis(list<const CalItem*> &tlist, float etaJet,
415  float phiJet,float cone_radius, float zvertex_in) const;
416 
417 public:
418 
419  vector< TemporaryJet > TempColl;
420 
421 };
422  /////////////////////////////////////////////////////////
423 
424 //template < class CalItem,class CalItemAddress,class CalIClusterChunk >
425 template < class CalItem >
426 //void ConeClusterAlgo <CalItem,CalItemAddress,CalIClusterChunk >::
427 void ConeClusterAlgo <CalItem >::
428 getItemsInCone(list<const CalItem*> &tlist, float etaJet, float phiJet,
429  float cone_radius, float zvertex_in) const {
430 //
431 // provide the list of Items (towers, Cells...) containing the energy from a
432 // jet of a given cone size
433 //
434  float ZVERTEX_MAX=200.;
435  float DMIN=80.;
436  float DMAX=360.;
437  float THETA_margin=0.022;
438 
439  float zvertex=zvertex_in;
440  float d1,d2;
441  float phi_d1, phi_d2;
442  float theta_E1, r1, r2, z1, z2;
443  float theta_d1, theta_d2, eta_d1, eta_d2;
444 
445  // Ignore very large vertex positions
446  if (fabs(zvertex) > ZVERTEX_MAX ) zvertex=0.0;
447 
448  if (zvertex >=0. ) {
449  d1=fabs(DMIN-zvertex);
450  d2=fabs(DMAX+zvertex);
451  } else {
452  d1=fabs(DMAX-zvertex);
453  d2=fabs(DMIN+zvertex);
454  }
455 
456  // calculate theta of physics cone and find which eta's this intercepts
457  // a the maximum points
458  phi_d1 = phiJet+cone_radius;
459  //theta_E1 = kinem::theta(etaJet+cone_radius);
460  theta_E1 = inline_maths::theta(etaJet+cone_radius);
461  z1 = zvertex+d1*cos(theta_E1);
462  r1 = d1*sin(theta_E1);
463 
464  phi_d2 = phiJet-cone_radius;
465  //theta_E1 = kinem::theta(etaJet-cone_radius);
466  theta_E1 = inline_maths::theta(etaJet-cone_radius);
467  z2 = zvertex+d2*cos(theta_E1);
468  r2 = d2*sin(theta_E1);
469 
470  // maximum spread in detector theta
471  theta_d1 = atan2(r1, z1);
472  theta_d2 = atan2(r2, z2);
473 
474  // make sure they stay in the calorimeter
475  theta_d1=max(theta_d1, THETA_margin);
476  theta_d2=max(theta_d2, THETA_margin);
477  //theta_d1=min(kinem::PI-(double)THETA_margin, (double)theta_d1);
478  theta_d1=min(inline_maths::PI-(double)THETA_margin, (double)theta_d1);
479  //theta_d2=min(kinem::PI-(double)THETA_margin, (double)theta_d2);
480  theta_d2=min(inline_maths::PI-(double)THETA_margin, (double)theta_d2);
481 
482  //eta_d1 = kinem::eta(theta_d1);
483  eta_d1 = inline_maths::eta(theta_d1);
484  //eta_d2 = kinem::eta(theta_d2);
485  eta_d2 = inline_maths::eta(theta_d2);
486 
487 
488  typename list<const CalItem*>::iterator it;
489  for (it=tlist.begin() ; it != tlist.end() ; ) {
490  //float eta_cur= (*it)->eta();
491  //float phi_cur= (*it)->phi();
492  float pz[4];
493  (*it)->p4vec(pz);
494  float eta_cur= E2eta(pz);
495  float phi_cur= E2phi(pz);
496 
497  bool accepted = eta_cur < eta_d1 && eta_cur > eta_d2;
498  //if ( phi_d2>0 && phi_d1<kinem::TWOPI ) {
499  if ( phi_d2>0 && phi_d1<inline_maths::TWOPI ) {
500  accepted = accepted && phi_cur<phi_d1 && phi_cur>phi_d2;
501  }
502  else{ // case the cone overlap the phi=0=2pi line
503  if ( phi_d2>0 ){
504  accepted = accepted &&
505  //((phi_cur>phi_d2 && phi_cur<kinem::TWOPI) || phi_cur<phi_d1-kinem::TWOPI);
506  ((phi_cur>phi_d2 && phi_cur<inline_maths::TWOPI) || phi_cur<phi_d1-inline_maths::TWOPI);
507  }
508  else{
509  accepted = accepted &&
510  //((phi_cur<phi_d1 && phi_cur>0) || phi_cur>phi_d2+kinem::TWOPI);
511  ((phi_cur<phi_d1 && phi_cur>0) || phi_cur>phi_d2+inline_maths::TWOPI);
512  }
513  }
514  if ( ! accepted ) it = tlist.erase(it);
515  else ++it;
516 
517  }
518 }
519  /////////////////////////////////////////////////////////
520 //template < class CalItem,class CalItemAddress,class CalIClusterChunk >
521 template < class CalItem >
522 //void ConeClusterAlgo <CalItem,CalItemAddress,CalIClusterChunk >::
523 void ConeClusterAlgo <CalItem>::
524 getItemsInCone_bis(list<const CalItem*> &tlist, float etaJet, float phiJet,
525  float cone_radius, float zvertex_in) const {
526 //
527 // provide the list of Items (towers, Cells...) containing the energy from a
528 // jet of a given cone size
529 //
530 // WARNING: this is only to be used to compare to RUN I cone jets
531 //
532  float ZVERTEX_MAX=200.;
533  float DMIN=80.;
534  float DMAX=360.;
535  float THETA_margin=0.022;
536 
537  float zvertex=zvertex_in;
538  float d1,d2;
539  float phi_d1, phi_d2;
540  float theta_E1, r1, r2, z1, z2;
541  float theta_d1, theta_d2, eta_d1, eta_d2;
542 
543  // Ignore very large vertex positions
544  if (fabs(zvertex) > ZVERTEX_MAX ) zvertex=0.0;
545 
546  if (zvertex >=0. ) {
547  d1=fabs(DMIN-zvertex);
548  d2=fabs(DMAX+zvertex);
549  } else {
550  d1=fabs(DMAX-zvertex);
551  d2=fabs(DMIN+zvertex);
552  }
553 
554  // calculate theta of physics cone and find which eta's this intercepts
555  // a the maximum points
556 
557  phi_d1 = phiJet+cone_radius;
558  //theta_E1 = kinem::theta(etaJet+cone_radius);
559  theta_E1 = inline_maths::theta(etaJet+cone_radius);
560  z1 = zvertex+d1*cos(theta_E1);
561  r1 = d1*sin(theta_E1);
562 
563  phi_d2 = phiJet-cone_radius;
564  //theta_E1 = kinem::theta(etaJet-cone_radius);
565  theta_E1 = inline_maths::theta(etaJet-cone_radius);
566  z2 = zvertex+d2*cos(theta_E1);
567  r2 = d2*sin(theta_E1);
568 
569  // maximum spread in detector theta
570 
571  theta_d1 = atan2(r1, z1);
572  theta_d2 = atan2(r2, z2);
573 
574  // make sure they stay in the calorimeter
575 
576  theta_d1=max(theta_d1, THETA_margin);
577  theta_d2=max(theta_d2, THETA_margin);
578  //theta_d1=min(kinem::PI-(double)THETA_margin, (double)theta_d1);
579  theta_d1=min(inline_maths::PI-(double)THETA_margin, (double)theta_d1);
580  //theta_d2=min(kinem::PI-(double)THETA_margin, (double)theta_d2);
581  theta_d2=min(inline_maths::PI-(double)THETA_margin, (double)theta_d2);
582 
583 
584  //eta_d1 = kinem::eta(theta_d1);
585  eta_d1 = inline_maths::eta(theta_d1);
586  //eta_d2 = kinem::eta(theta_d2);
587  eta_d2 = inline_maths::eta(theta_d2);
588 
589  float signe;
590 
591  if( eta_d1>=0.0 ) signe= 1.0;
592  else signe= -1.0;
593  int ietaMAX= eta_d1/0.1+signe;
594  if(fabs(eta_d1)>=4.45) ietaMAX= 37*signe;
595  else if(fabs(eta_d1)>=4.1) ietaMAX= 36*signe;
596  else if(fabs(eta_d1)>=3.7) ietaMAX= 35*signe;
597  else if(fabs(eta_d1)>=3.42) ietaMAX= 34*signe;
598  else if(fabs(eta_d1)>=3.2) ietaMAX= 33*signe;
599 
600  if( eta_d2>=0.0 ) signe= 1.0;
601  else signe= -1.0;
602  int ietaMIN= eta_d2/0.1+signe;
603  if(fabs(eta_d2)>=4.45) ietaMIN= 37*signe;
604  else if(fabs(eta_d2)>=4.1) ietaMIN= 36*signe;
605  else if(fabs(eta_d2)>=3.7) ietaMIN= 35*signe;
606  else if(fabs(eta_d2)>=3.42) ietaMIN= 34*signe;
607  else if(fabs(eta_d2)>=3.2) ietaMIN= 33*signe;
608 
609  //int iphiMAX= 64*phi_d1/(2.*kinem::PI)+1;
610  int iphiMAX= 64*phi_d1/(2.*inline_maths::PI)+1;
611  //int iphiMIN= 64*phi_d2/(2.*kinem::PI)+1;
612  int iphiMIN= 64*phi_d2/(2.*inline_maths::PI)+1;
613 
614  typename list<const CalItem*>::iterator it;
615  for (it=tlist.begin() ; it != tlist.end() ; ) {
616  //float eta_cur= (*it)->eta();
617  //float phi_cur= (*it)->phi();
618  int ieta= (*it)->address().ieta();
619  int iphi= (*it)->address().iphi();
620 
621  bool accepted = ieta<ietaMAX && ieta>ietaMIN;
622  if ( iphiMIN>0 && iphiMAX<=64 ) {
623  accepted = accepted && iphi<iphiMAX && iphi>iphiMIN;
624  }
625  else{ // case the cone overlap the phi=0=2pi line
626  if ( iphiMIN>0 ){
627  accepted = accepted &&
628  ((iphi>iphiMIN && iphi<=64) || iphi<iphiMAX-64);
629  }
630  else{
631  accepted = accepted &&
632  ((iphi<iphiMAX && iphi>0) || iphi>iphiMIN+64);
633  }
634  }
635  if ( ! accepted ) it = tlist.erase(it);
636  else ++it;
637 
638  }
639 }
640  /////////////////////////////////////////////////////////
641 //template < class CalItem,class CalItemAddress,class CalIClusterChunk >
642 template < class CalItem >
643 //inline void ConeClusterAlgo <CalItem,CalItemAddress,CalIClusterChunk >::
644 inline void ConeClusterAlgo <CalItem >::
645 print(ostream &os) const {
646  os<<endl<<" CONE ALGORITHM, cone radius= "<<_CONErad<<endl<<
647  " min E_T fraction= "<<_JETmne<<endl<<
648  " minimum Delta_R separation between cones= "<<_TWOrad<<endl<<
649  " shared E_T fraction threshold for combining jets= "<<_SPLifr<<endl;
650 }
651  /////////////////////////////////////////////////////////
652 
653 //template < class CalItem,class CalItemAddress,class CalIClusterChunk >
654 template < class CalItem >
655 //void ConeClusterAlgo <CalItem,CalItemAddress,CalIClusterChunk >::
656 void ConeClusterAlgo <CalItem >::
657 makeClusters(//const EnergyClusterReco* r,
658  std::list<CalItem> &jets,
659  list<const CalItem*> &itemlist, float Zvertex
660  //, const EnergyClusterCollection<CalItemAddress> &preclu,
661  //CalIClusterChunk* chunkptr
662  //) const {
663  ) {
664 
665  // create an energy cluster collection for jets
666  //EnergyClusterCollection<CalItemAddress>* ptrcol;
667  //r->createClusterCollection(chunkptr, ptrcol);
668  std::vector<const CalItem*> ecv;
669  for ( typename std::list<const CalItem*>::iterator it = itemlist.begin();
670  it != itemlist.end(); it++) {
671  ecv.push_back(*it);
672  }
673 
674 
675  // Initialize
676  float Rcut= 1.E-06;
677  if(_Increase_Delta_R) Rcut= 1.E-04;
678  bool nojets;
679 
680  //vector< TemporaryJet > TempColl;
681  list< pair<float,float> > LTrack;
682 
683  // get a vector with pointers to EnergyCluster in the collection
684  //vector<const EnergyCluster<CalItemAddress>*> ecv;
685  //preclu.getClusters(ecv);
686 
687  // loop over all preclusters
688  //typename vector<const EnergyCluster<CalItemAddress>*>::iterator jclu;
689  typename std::vector<const CalItem*>::iterator jclu;
690  for( jclu=ecv.begin(); jclu!=ecv.end(); jclu++ ) {
691  ////const EnergyCluster<CalItemAddress>* ptr= *jclu;
692  const CalItem* ptr= *jclu;
693  //float PHIst= ptr->phi();
694  //float ETAst= ptr->eta();
695  float pz[4];
696  ptr->p4vec(pz);
697  float ETAst= E2eta(pz);
698  float PHIst= E2phi(pz);
699 
700  //cout << "seed 4-vec ";
701  //for ( int i = 0; i < 4; i++) cout << pz[i] << " ";
702  //cout << endl;
703 
704  nojets= false;
705  // check to see if precluster is too close to a found jet
706  if(_Kill_Far_Clusters) {
707  list< pair<float,float> >::iterator kj;
708  for(kj=LTrack.begin(); kj!=LTrack.end(); kj++) {
709  if(DELTA_r((*kj).first,ETAst,(*kj).second,PHIst)<_Far_Ratio*_CONErad) {
710  nojets= true;
711  //cout << "seed too close ! skip " << endl;
712  break;
713  }
714  }
715  }
716  if( nojets==false ) {
717  TemporaryJet TJet(ETAst,PHIst);
718  list< pair<int,float> > JETshare;
719 
720  // start of cone building loop
721  int trial= 0;
722  do{
723  trial++;
724  //cout << " trial " << trial << endl;
725 
726  ETAst= TJet.Eta();
727  PHIst= TJet.Phi();
728  TJet.erase();
729 
730  //if(PHIst > kinem::TWOPI) PHIst= PHIst-kinem::TWOPI;
731  if(PHIst > inline_maths::TWOPI) PHIst= PHIst-inline_maths::TWOPI;
732  //if(PHIst < 0.0 ) PHIst= PHIst+kinem::TWOPI;
733  if(PHIst < 0.0 ) PHIst= PHIst+inline_maths::TWOPI;
734  //if( PHIst>kinem::TWOPI || PHIst<0.0 ) {
735  if( PHIst>inline_maths::TWOPI || PHIst<0.0 ) {
736  TJet.setEtaPhiEt(0.0,0.0,0.0);
737  break; // end loop do (illegal jet PHI)
738  }
739  TJet.setEtaPhiEt(ETAst,PHIst,0.0);
740 
741  // calculate eta & phi limits for cone
742  list<const CalItem*> Twlist(itemlist);
743 
744  getItemsInCone(Twlist,ETAst,PHIst,_CONErad,Zvertex);
745  // only to compare with RUN I cone jets ! getItemsInCone_bis(Twlist,ETAst,PHIst,_CONErad,Zvertex);
746 
747  // loop over the possible items for this cone
748  typename list<const CalItem*>::iterator tk;
749  for( tk= Twlist.begin(); tk!=Twlist.end(); tk++ ) {
750  float ETk =(*tk)->pT();
751  // now done in CalCell/CalTower if((*tk)->E() < 0.0) ETk= -ETk;
752 
753  if( ETk > _Eitem_Negdrop ) {
754  //float ETAk=(*tk)->eta();
755  //float PHIk=(*tk)->phi();
756  float pz[4];
757  (*tk)->p4vec(pz);
758  float ETAk= E2eta(pz);
759  float PHIk= E2phi(pz);
760 
761  float dphi= fabs(PHIk-PHIst);
762  //if(dphi > kinem::TWOPI-dphi) {
763  if(dphi > inline_maths::TWOPI-dphi) {
764  //if(PHIst < PHIk) PHIk= PHIk-kinem::TWOPI;
765  if(PHIst < PHIk) PHIk= PHIk-inline_maths::TWOPI;
766  //else PHIk= PHIk+kinem::TWOPI;
767  else PHIk= PHIk+inline_maths::TWOPI;
768  }
769 
770  if( R2_bis(ETAk,PHIk,ETAst,PHIst) <= _CONErad*_CONErad ) {
771  TJet.addItem(*tk);
772  }
773  }
774  }// end loop tk
775 
776  if(TJet.updateEtaPhiEt()==false) {
777  //cout << " negative E jet ! drop " << endl;
778  break;
779  }
780 
781  // require some minimum ET on every iteration
782  if(_Jet_Et_Min_On_Iter) {
783  if( TJet.Et() < _JETmne*_Et_Min_Ratio ) {
784  //cout << " too low ET jet ! drop " << endl;
785  break; // end loop trial
786  }
787  }
788 
789  //cout << " R2 = " << R2_bis(TJet.Eta(),TJet.Phi(),ETAst,PHIst) <<
790  // " Rcut = " << Rcut << endl;
791  }while(R2_bis(TJet.Eta(),TJet.Phi(),ETAst,PHIst)>=Rcut && trial<=50);
792 
793  if( TJet.Et() >= _JETmne ) {
794  //cout << " jet accepted will check for overlaps " << endl;
795  if(_D0_Angle) TJet.D0_Angle_updateEtaPhi();
796  //cout << " after TJet.D0_Angle_updateEtaPhi() " << endl;
797 
798  // item also in another jet
799  list<const CalItem*> Lst;
800  TJet.getItems(Lst);
801  typename list<const CalItem*>::iterator tk;
802  for(tk=Lst.begin(); tk!=Lst.end(); tk++) {
803  float ETk=(*tk)->pT();
804  // now done in CalCell/CalTower if((*tk)->E() < 0.0) ETk= -ETk;
805  for(unsigned int kj=0; kj<TempColl.size(); kj++) {
806  if(TempColl[kj].ItemInJet(*tk)==true) {
807  list< pair<int,float> >::iterator pit;
808  bool jetok= false;
809  for(pit=JETshare.begin(); pit!=JETshare.end();pit++) {
810  if((*pit).first == (int) kj) {
811  jetok= true;
812  (*pit).second+= ETk;
813  break;
814  }
815  }
816  if(jetok==false) JETshare.push_back(make_pair(kj,ETk));
817  }
818  }
819  }
820 
821  if(JETshare.size() >0) {
822  list< pair<int,float> >::iterator pit;
823  float Ssum= 0.0;
824  list< pair<int,float> >::iterator pitMAX=JETshare.begin();
825  for(pit=JETshare.begin(); pit!=JETshare.end(); pit++) {
826  Ssum+= (*pit).second;
827  if((*pit).second > (*pitMAX).second) pitMAX= pit;
828  }
829 
830  //int IJET= (*pitMAX).first;
831  bool splshr;
832  float Eleft= fabs(TJet.Et()-Ssum);
833  float Djets= TempColl[(*pitMAX).first].dist_R2(TJet);
834  if(Djets <= _TWOrad || Eleft <= _Thresh_Diff_Et) {
835  TJet.erase();
836  splshr= false;
837  }
838  else {
839  float SharedFr=Ssum/min(TempColl[(*pitMAX).first].Et(),TJet.Et());
840  if(JETshare.size() >1) {
841  typename list<const CalItem*>::iterator tk;
842  for(tk=TJet.LItems().begin(); tk!=TJet.LItems().end(); ) {
843  bool found = false;
844  list< pair<int,float> >::iterator pit;
845  for(pit=JETshare.begin(); pit!=JETshare.end();pit++) {
846  if((*pit).first!=(*pitMAX).first) {
847  if(TempColl[(*pit).first].ItemInJet(*tk)==true) {
848  tk = TJet.LItems().erase(tk);
849  found = true;
850  break;
851  }
852  }
853  }
854  if ( !found ) ++tk;
855  }
856  }
857 
858  splshr= TempColl[(*pitMAX).first].share_jets(TJet,SharedFr,_SPLifr);
859 
860  if( splshr==true ) {
861  //cout << " jet splitted due to overlaps " << endl;
862  TempColl[(*pitMAX).first].updateEtaPhiEt();
863  TJet.updateEtaPhiEt();
864  if(_D0_Angle) TJet.D0_Angle_updateEtaPhi();
865  if(_D0_Angle) TempColl[(*pitMAX).first].D0_Angle_updateEtaPhi();
866  TempColl.push_back(TJet);
867  LTrack.push_back(make_pair(TJet.Eta(),TJet.Phi()));
868  }
869  else {
870  //cout << " jet merged due to overlaps " << endl;
871  TempColl[(*pitMAX).first].updateEtaPhiEt();
872  if(_D0_Angle) TempColl[(*pitMAX).first].D0_Angle_updateEtaPhi();
873  }
874  }
875  }
876  else {
877  TJet.updateEtaPhiEt();
878  if(_D0_Angle) TJet.D0_Angle_updateEtaPhi();
879  TempColl.push_back(TJet);
880  LTrack.push_back(make_pair(TJet.Eta(),TJet.Phi()));
881  }
882  } //JETmne
883  } //nojets
884  }// end loop jclu
885 
886  for(unsigned int i=0; i<TempColl.size(); i++) {
887  //EnergyCluster<CalItemAddress>* ptrclu;
888  CalItem ptrclu;
889  //r->createCluster(ptrcol,ptrclu);
890  list<const CalItem*> Vi;
891  TempColl[i].getItems(Vi);
892  typename list<const CalItem*>::iterator it;
893  for(it=Vi.begin(); it!=Vi.end(); it++) {
894  const CalItem* ptr= *it;
895  //CalItemAddress addr= ptr->address();
896  float p[4];
897  ptr->p4vec(p);
898  //float emE= ptr->emE();
899  //r->addClusterItem(ptrclu,addr,p,emE);
900  ptrclu.Add(*ptr);
901  }
902  jets.push_back(ptrclu);
903  }
904 
905 }// end
906 
907 } //namespace d0runi
908 
909 FASTJET_END_NAMESPACE
910 
911 #endif // CONECLUSTERALGO_H
912 
913 
914 
Generated on Tue Jan 6 2015 16:39:23 for FastJet by   doxygen 1.8.8