Fork me on GitHub

source: git/external/fastjet/plugins/D0RunICone/ConeClusterAlgo.hpp@ c18dca6

Last change on this file since c18dca6 was d7d2da3, checked in by pavel <pavel@…>, 12 years ago

move branches/ModularDelphes to trunk
  • Property mode set to 100644
File size: 26.4 KB
RevLine 
[d7d2da3]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
72FASTJET_BEGIN_NAMESPACE
73
74namespace d0runi{
75
76using namespace std;
77
78//some utility functions
79inline float R2(float eta1, float phi1, float eta2, float phi2) {
80 return (eta1-eta2)*(eta1-eta2)+(phi1-phi2)*(phi1-phi2); }
81
82inline 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
87inline 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
93inline 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
116inline 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 >
136template < class CalItem >
137class 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
148ConeClusterAlgo() {}
149
150//constructor for cone jet algorithm
151ConeClusterAlgo( 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)
168ConeClusterAlgo( 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
190void 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
199void 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
417public:
418
419 vector< TemporaryJet > TempColl;
420
421};
422 /////////////////////////////////////////////////////////
423
424//template < class CalItem,class CalItemAddress,class CalIClusterChunk >
425template < class CalItem >
426//void ConeClusterAlgo <CalItem,CalItemAddress,CalIClusterChunk >::
427void ConeClusterAlgo <CalItem >::
428getItemsInCone(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 >
521template < class CalItem >
522//void ConeClusterAlgo <CalItem,CalItemAddress,CalIClusterChunk >::
523void ConeClusterAlgo <CalItem>::
524getItemsInCone_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 >
642template < class CalItem >
643//inline void ConeClusterAlgo <CalItem,CalItemAddress,CalIClusterChunk >::
644inline void ConeClusterAlgo <CalItem >::
645print(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 >
654template < class CalItem >
655//void ConeClusterAlgo <CalItem,CalItemAddress,CalIClusterChunk >::
656void ConeClusterAlgo <CalItem >::
657makeClusters(//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
909FASTJET_END_NAMESPACE
910
911#endif // CONECLUSTERALGO_H
912
913
914
Note: See TracBrowser for help on using the repository browser.