FastJet 3.0.2
|
00001 ////////////////////////////////////////////////////////////// 00002 // File: ConeClusterAlgo.hpp 00003 // 00004 // Author: G. Le Meur & F. Touze 00005 // 00006 // Created: 15-JUNE-1998 00007 // 00008 // Purpose: make jet clusters using fixed cone like algorithm 00009 // implemented in RUNI. 00010 // 00011 // Modified: 00012 // 28-OCT-1998 to use KinemUtil (S. Protopopescu) 00013 // 8-JAN-1999: Laurent Duflot 00014 // . correct bugs in getItemsInCone and updateEtaPhiEt for jets 00015 // overlapping the phi=0 line 00016 // . change abs(float) to fabs(float) 00017 // 1-NOV-1999: Laurent Duflot 00018 // . correct bug in makeCluster: when the temporary jet was emptied the eta 00019 // and phi were not set again. The main effect was a nearly zero 00020 // efficiency for jets at phi=pi (as seen by Volker Buescher) 00021 // 25-JAN-2000: Francois Touze 00022 // . change in updateEtaPhiEt : the method E() which returns energy doesn't 00023 // exist in MCparticle classe,... so use the 4-momentum components 00024 // . declare const the EnergyClusterCollection of seeds in makeClusters 00025 // 01-FEB-2000: Laurent Duflot 00026 // . add a missing break statement in the removal of shared items. Caused 00027 // an infinite loop on some events. 00028 // . correct typo in variable name. Change a variable name that was in 00029 // French. 00030 // . leave some debug printout (commented) 00031 // 15-Sep-2009 Lars Sonnenschein 00032 // extracted from D0 software framework and modified to remove subsequent dependencies 00033 // 00034 // 00035 // This file is distributed with FastJet under the terms of the GNU 00036 // General Public License (v2). Permission to do so has been granted 00037 // by Lars Sonnenschein and the D0 collaboration (see COPYING for 00038 // details) 00039 // 00040 // History of changes in FastJet compared to the original version of 00041 // ConeClusterAlgo.hpp 00042 // 00043 // 2011-12-13 Gregory Soyez <soyez@fastjet.fr> 00044 // 00045 // * added license information 00046 // 00047 // 2011-10-06 Gregory Soyez <soyez@fastjet.fr> 00048 // 00049 // * put the code in the fastjet::d0runi namespace 00050 // 00051 ////////////////////////////////////////////////////////////// 00052 00053 //#ifndef CONECLUSTERALGO_H 00054 //#define CONECLUSTERALGO_H 00055 00056 #ifndef D0RunIconeJets_CONECLUSTERALGO_H 00057 #define D0RunIconeJets_CONECLUSTERALGO_H 00058 00059 00060 //#include "EnergyClusterReco.hpp" 00061 #include <vector> 00062 #include <list> 00063 #include <utility> 00064 //#include "kinem_util/AnglesUtil.hpp" 00065 00066 #include <algorithm> 00067 #include <iostream> 00068 00069 #include "inline_maths.h" 00070 #include <fastjet/internal/base.hh> 00071 00072 FASTJET_BEGIN_NAMESPACE 00073 00074 namespace d0runi{ 00075 00076 using namespace std; 00077 00078 //some utility functions 00079 inline float R2(float eta1, float phi1, float eta2, float phi2) { 00080 return (eta1-eta2)*(eta1-eta2)+(phi1-phi2)*(phi1-phi2); } 00081 00082 inline float R2_bis(float eta1, float phi1, float eta2, float phi2) { 00083 //float dphi = kinem::delta_phi(phi1,phi2); 00084 float dphi = inline_maths::delta_phi(phi1,phi2); 00085 return (eta1-eta2)*(eta1-eta2)+dphi*dphi; } 00086 00087 inline float DELTA_r(float eta1,float eta2,float phi1,float phi2) { 00088 //float dphi = kinem::delta_phi(phi1,phi2); 00089 float dphi = inline_maths::delta_phi(phi1,phi2); 00090 return sqrt((eta1-eta2)*(eta1-eta2)+dphi*dphi); 00091 } 00092 00093 inline float E2eta(float* p) { 00094 float small= 1.E-05; 00095 float E[3]; 00096 if(p[3] < 0.0) { 00097 E[0]= -p[0]; 00098 E[1]= -p[1]; 00099 E[2]= -p[2]; 00100 } 00101 else { 00102 E[0]= p[0]; 00103 E[1]= p[1]; 00104 E[2]= p[2]; 00105 } 00106 float pperp= sqrt(E[0]*E[0]+E[1]*E[1])+small; 00107 float ptotal= sqrt(E[0]*E[0]+E[1]*E[1]+E[2]*E[2])+small; 00108 //float theta= atan2(pperp,E[2]); 00109 00110 float eta= 0.0; 00111 if(E[2] > 0.0) eta= log((ptotal+E[2])/pperp); 00112 else eta= log(pperp/(ptotal-E[2])); 00113 return eta; 00114 } 00115 00116 inline float E2phi(float* p) { 00117 float small= 1.E-05; 00118 float E[3]; 00119 if(p[3] < 0.0) { 00120 E[0]= -p[0]; 00121 E[1]= -p[1]; 00122 E[2]= -p[2]; 00123 } 00124 else { 00125 E[0]= p[0]; 00126 E[1]= p[1]; 00127 E[2]= p[2]; 00128 } 00129 float phi= atan2(E[1],E[0]+small); 00130 //if(phi < 0.0) phi+=kinem::TWOPI; 00131 if (phi < 0.0) phi += inline_maths::TWOPI; 00132 return phi; 00133 } 00134 00135 //template < class CalItem,class CalItemAddress,class CalIClusterChunk > 00136 template < class CalItem > 00137 class ConeClusterAlgo { 00138 // 00139 // Purpose: make calorimeter clusters using a cone algorithm from 00140 // preclusters created previously by the class ConePreClusterAlgo. 00141 // Items must have addresses and 4-momenta. 00142 // The algorithm is implemented with a template function makeClusters. 00143 // 00144 public : 00145 00146 00147 //default constructor 00148 ConeClusterAlgo() {} 00149 00150 //constructor for cone jet algorithm 00151 ConeClusterAlgo( float CONErad,float JETmne,float SPLifr): 00152 _CONErad(fabs(CONErad)), 00153 _JETmne(JETmne), 00154 _SPLifr(SPLifr), 00155 _TWOrad(0.), 00156 _D0_Angle(false), 00157 _Increase_Delta_R(true), 00158 _Kill_Far_Clusters(true), 00159 _Jet_Et_Min_On_Iter(true), 00160 _Far_Ratio(0.5), 00161 _Eitem_Negdrop(-1.0), 00162 _Et_Min_Ratio(0.5), 00163 _Thresh_Diff_Et(0.01) 00164 {} 00165 00166 //changing default thresholds & parameters 00167 // (declared by PARAMETER in RUNI code) 00168 ConeClusterAlgo( float CONErad,float JETmne,float SPLifr,float TWOrad, 00169 float Tresh_Diff_Et,bool D0_Angle,bool Increase_Delta_R, 00170 bool Kill_Far_Clusters,bool Jet_Et_Min_On_Iter, 00171 float Far_Ratio,float Eitem_Negdrop,float Et_Min_Ratio ): 00172 _CONErad(fabs(CONErad)), 00173 _JETmne(JETmne), 00174 _SPLifr(SPLifr), 00175 _TWOrad(TWOrad), 00176 _D0_Angle(D0_Angle), 00177 _Increase_Delta_R(Increase_Delta_R), 00178 _Kill_Far_Clusters(Kill_Far_Clusters), 00179 _Jet_Et_Min_On_Iter(Jet_Et_Min_On_Iter), 00180 _Far_Ratio(Far_Ratio), 00181 _Eitem_Negdrop(Eitem_Negdrop), 00182 _Et_Min_Ratio(Et_Min_Ratio), 00183 _Thresh_Diff_Et(Tresh_Diff_Et) 00184 {} 00185 00186 //destructor 00187 ~ConeClusterAlgo() {} 00188 00189 //to make jet clusters using cone algorithm 00190 void makeClusters(//const EnergyClusterReco* r, 00191 std::list<CalItem> &jets, 00192 list<const CalItem*> &itemlist, float Zvertex 00193 //, const EnergyClusterCollection<CalItemAddress> &preclu, 00194 //CalIClusterChunk* chunkptr 00195 //) const; 00196 ); 00197 00198 //print parameters of the algorithm 00199 void print(ostream &os)const; 00200 00201 //vector< TemporaryJet > TempColl; 00202 00203 00204 private : 00205 00206 float _CONErad; 00207 float _JETmne; 00208 float _SPLifr; 00209 00210 float _TWOrad; 00211 bool _D0_Angle; 00212 bool _Increase_Delta_R; 00213 bool _Kill_Far_Clusters; 00214 bool _Jet_Et_Min_On_Iter; 00215 float _Far_Ratio; 00216 float _Eitem_Negdrop; 00217 float _Et_Min_Ratio; 00218 float _Thresh_Diff_Et; 00219 00220 class TemporaryJet { 00221 00222 public: 00223 00224 00225 00226 TemporaryJet() {} 00227 00228 TemporaryJet(float eta,float phi) { 00229 _Eta=eta; 00230 _Phi=phi; 00231 } 00232 00233 ~TemporaryJet() {} 00234 00235 void addItem(const CalItem* tw) { 00236 _LItems.push_back(tw); 00237 } 00238 00239 void setEtaPhiEt(float eta,float phi,float pT) { 00240 _Eta= eta; 00241 _Phi= phi; 00242 _Et = pT; 00243 } 00244 00245 void erase() { 00246 _LItems.erase(_LItems.begin(),_LItems.end()); 00247 _Eta= 0.0; 00248 _Phi= 0.0; 00249 _Et = 0.0; 00250 } 00251 00252 bool share_jets(TemporaryJet &NewJet,float SharedFr,float SPLifr) { 00253 // 00254 // combined 00255 // 00256 if(SharedFr >= SPLifr) { 00257 typename list<const CalItem*>::iterator it; 00258 typename list<const CalItem*>::iterator end_of_old=_LItems.end(); 00259 for(it=NewJet._LItems.begin(); it!=NewJet._LItems.end(); it++) { 00260 typename list<const CalItem*>::iterator 00261 where = find(_LItems.begin(),end_of_old,*it); 00262 // if the item is not shared, add to this jet 00263 if(where == end_of_old) { 00264 _LItems.push_back(*it); 00265 } 00266 } 00267 NewJet.erase(); 00268 return false; 00269 } else { 00270 // 00271 // split 00272 // 00273 typename list<const CalItem*>::iterator it; 00274 for(it=NewJet._LItems.begin(); it!=NewJet._LItems.end(); ) { 00275 typename list<const CalItem*>::iterator 00276 where = find(_LItems.begin(),_LItems.end(),*it); 00277 if(where != _LItems.end()) { 00278 //float EtaItem=(*it)->eta(); 00279 //float PhiItem=(*it)->phi(); 00280 // stay closer to the RUNI conventions for negative E cells 00281 float pz[4]; 00282 (*it)->p4vec(pz); 00283 float EtaItem= E2eta(pz); 00284 float PhiItem= E2phi(pz); 00285 00286 float RadOld=R2_bis(_Eta,_Phi,EtaItem,PhiItem); 00287 float RadNew=R2_bis(NewJet.Eta(),NewJet.Phi(),EtaItem,PhiItem); 00288 if (RadNew > RadOld) { 00289 it = NewJet._LItems.erase(it); 00290 } 00291 else { 00292 _LItems.erase(where); 00293 ++it; 00294 } 00295 } 00296 else ++it; 00297 } 00298 return true; 00299 } 00300 } 00301 00302 00303 float dist_R2(TemporaryJet &jet) const { 00304 float deta= _Eta-jet.Eta(); 00305 //float dphi= kinem::delta_phi(_Phi,jet.Phi()); 00306 float dphi= inline_maths::delta_phi(_Phi,jet.Phi()); 00307 return (deta*deta+dphi*dphi); 00308 } 00309 00310 bool ItemInJet(const CalItem* tw) const { 00311 typename list<const CalItem*>::const_iterator 00312 where= find(_LItems.begin(),_LItems.end(),tw); 00313 if(where != _LItems.end()) return true; 00314 else return false; 00315 } 00316 00317 bool updateEtaPhiEt() { 00318 float ETsum = 0.0; 00319 float ETAsum= 0.0; 00320 float PHIsum= 0.0; 00321 float Esum= 0.0; 00322 typename list<const CalItem*>::iterator it; 00323 for(it=_LItems.begin(); it!=_LItems.end(); it++) { 00324 float ETk = (*it)->pT(); 00325 // now done in CalCell/CalTower if((*it)->E() < 0.0) ETk= -ETk; 00326 00327 //float ETAk= (*it)->eta(); 00328 //float PHIk= (*it)->phi(); 00329 float pz[4]; 00330 (*it)->p4vec(pz); 00331 float ETAk= E2eta(pz); 00332 // take care of the phi=0=2pi problem 00333 float PHIk= E2phi(pz); 00334 //if(fabs(PHIk-_Phi) > kinem::TWOPI-fabs(PHIk-_Phi)) 00335 if(fabs(PHIk-_Phi) > inline_maths::TWOPI-fabs(PHIk-_Phi)) 00336 { 00337 if(_Phi < PHIk) 00338 { 00339 //PHIk -= kinem::TWOPI; 00340 PHIk -= inline_maths::TWOPI; 00341 } 00342 else 00343 { 00344 //PHIk += kinem::TWOPI; 00345 PHIk += inline_maths::TWOPI; 00346 } 00347 } 00348 ETAsum+= ETAk*ETk; 00349 PHIsum+= PHIk*ETk; 00350 ETsum += ETk; 00351 // Esum+=(*it)->E(); use 4-momentum components 00352 Esum+= pz[3]; 00353 } 00354 if(ETsum <= 0.0) { 00355 _Eta= 0.0; 00356 _Phi= 0.0; 00357 _Et = 0.0; 00358 _E=0.; 00359 return false; 00360 } 00361 else { 00362 _Eta= ETAsum/ETsum; 00363 _Phi= PHIsum/ETsum; 00364 //if ( _Phi<0 ) _Phi+=kinem::TWOPI; 00365 if ( _Phi<0 ) _Phi+=inline_maths::TWOPI; 00366 _Et = ETsum; 00367 _E = Esum; 00368 return true; 00369 } 00370 } 00371 00372 void D0_Angle_updateEtaPhi() { 00373 float EXsum = 0.0; 00374 float EYsum = 0.0; 00375 float EZsum = 0.0; 00376 typename list<const CalItem*>::iterator it; 00377 for(it=_LItems.begin(); it!=_LItems.end(); it++) { 00378 float p[4]; 00379 (*it)->p4vec(p); 00380 EXsum += p[0]; 00381 EYsum += p[1]; 00382 EZsum += p[2]; 00383 } 00384 //_Phi=kinem::phi(EYsum,EXsum); 00385 _Phi=inline_maths::phi(EYsum,EXsum); 00386 //_Eta=kinem::eta(EXsum,EYsum,EZsum); 00387 _Eta=inline_maths::eta(EXsum,EYsum,EZsum); 00388 } 00389 00390 void getItems(list<const CalItem*> &ecv) const { 00391 ecv.clear(); //ls 27/Feb/2010 00392 typename list<const CalItem*>::const_iterator it; 00393 for(it=_LItems.begin(); it!=_LItems.end(); it++) { 00394 ecv.push_back(*it); 00395 } 00396 } 00397 00398 float Eta() {return _Eta;} 00399 float Phi() {return _Phi;} 00400 float Et() {return _Et;} 00401 float E() {return _E;} 00402 list<const CalItem*> &LItems() {return _LItems;} 00403 00404 private: 00405 list<const CalItem*> _LItems; 00406 float _Eta; 00407 float _Phi; 00408 float _Et; 00409 float _E; 00410 }; //class TemporaryJet 00411 00412 void getItemsInCone(list<const CalItem*> &tlist, float etaJet, float phiJet, 00413 float cone_radius, float zvertex_in) const; 00414 void getItemsInCone_bis(list<const CalItem*> &tlist, float etaJet, 00415 float phiJet,float cone_radius, float zvertex_in) const; 00416 00417 public: 00418 00419 vector< TemporaryJet > TempColl; 00420 00421 }; 00422 ///////////////////////////////////////////////////////// 00423 00424 //template < class CalItem,class CalItemAddress,class CalIClusterChunk > 00425 template < class CalItem > 00426 //void ConeClusterAlgo <CalItem,CalItemAddress,CalIClusterChunk >:: 00427 void ConeClusterAlgo <CalItem >:: 00428 getItemsInCone(list<const CalItem*> &tlist, float etaJet, float phiJet, 00429 float cone_radius, float zvertex_in) const { 00430 // 00431 // provide the list of Items (towers, Cells...) containing the energy from a 00432 // jet of a given cone size 00433 // 00434 float ZVERTEX_MAX=200.; 00435 float DMIN=80.; 00436 float DMAX=360.; 00437 float THETA_margin=0.022; 00438 00439 float zvertex=zvertex_in; 00440 float d1,d2; 00441 float phi_d1, phi_d2; 00442 float theta_E1, r1, r2, z1, z2; 00443 float theta_d1, theta_d2, eta_d1, eta_d2; 00444 00445 // Ignore very large vertex positions 00446 if (fabs(zvertex) > ZVERTEX_MAX ) zvertex=0.0; 00447 00448 if (zvertex >=0. ) { 00449 d1=fabs(DMIN-zvertex); 00450 d2=fabs(DMAX+zvertex); 00451 } else { 00452 d1=fabs(DMAX-zvertex); 00453 d2=fabs(DMIN+zvertex); 00454 } 00455 00456 // calculate theta of physics cone and find which eta's this intercepts 00457 // a the maximum points 00458 phi_d1 = phiJet+cone_radius; 00459 //theta_E1 = kinem::theta(etaJet+cone_radius); 00460 theta_E1 = inline_maths::theta(etaJet+cone_radius); 00461 z1 = zvertex+d1*cos(theta_E1); 00462 r1 = d1*sin(theta_E1); 00463 00464 phi_d2 = phiJet-cone_radius; 00465 //theta_E1 = kinem::theta(etaJet-cone_radius); 00466 theta_E1 = inline_maths::theta(etaJet-cone_radius); 00467 z2 = zvertex+d2*cos(theta_E1); 00468 r2 = d2*sin(theta_E1); 00469 00470 // maximum spread in detector theta 00471 theta_d1 = atan2(r1, z1); 00472 theta_d2 = atan2(r2, z2); 00473 00474 // make sure they stay in the calorimeter 00475 theta_d1=max(theta_d1, THETA_margin); 00476 theta_d2=max(theta_d2, THETA_margin); 00477 //theta_d1=min(kinem::PI-(double)THETA_margin, (double)theta_d1); 00478 theta_d1=min(inline_maths::PI-(double)THETA_margin, (double)theta_d1); 00479 //theta_d2=min(kinem::PI-(double)THETA_margin, (double)theta_d2); 00480 theta_d2=min(inline_maths::PI-(double)THETA_margin, (double)theta_d2); 00481 00482 //eta_d1 = kinem::eta(theta_d1); 00483 eta_d1 = inline_maths::eta(theta_d1); 00484 //eta_d2 = kinem::eta(theta_d2); 00485 eta_d2 = inline_maths::eta(theta_d2); 00486 00487 00488 typename list<const CalItem*>::iterator it; 00489 for (it=tlist.begin() ; it != tlist.end() ; ) { 00490 //float eta_cur= (*it)->eta(); 00491 //float phi_cur= (*it)->phi(); 00492 float pz[4]; 00493 (*it)->p4vec(pz); 00494 float eta_cur= E2eta(pz); 00495 float phi_cur= E2phi(pz); 00496 00497 bool accepted = eta_cur < eta_d1 && eta_cur > eta_d2; 00498 //if ( phi_d2>0 && phi_d1<kinem::TWOPI ) { 00499 if ( phi_d2>0 && phi_d1<inline_maths::TWOPI ) { 00500 accepted = accepted && phi_cur<phi_d1 && phi_cur>phi_d2; 00501 } 00502 else{ // case the cone overlap the phi=0=2pi line 00503 if ( phi_d2>0 ){ 00504 accepted = accepted && 00505 //((phi_cur>phi_d2 && phi_cur<kinem::TWOPI) || phi_cur<phi_d1-kinem::TWOPI); 00506 ((phi_cur>phi_d2 && phi_cur<inline_maths::TWOPI) || phi_cur<phi_d1-inline_maths::TWOPI); 00507 } 00508 else{ 00509 accepted = accepted && 00510 //((phi_cur<phi_d1 && phi_cur>0) || phi_cur>phi_d2+kinem::TWOPI); 00511 ((phi_cur<phi_d1 && phi_cur>0) || phi_cur>phi_d2+inline_maths::TWOPI); 00512 } 00513 } 00514 if ( ! accepted ) it = tlist.erase(it); 00515 else ++it; 00516 00517 } 00518 } 00519 ///////////////////////////////////////////////////////// 00520 //template < class CalItem,class CalItemAddress,class CalIClusterChunk > 00521 template < class CalItem > 00522 //void ConeClusterAlgo <CalItem,CalItemAddress,CalIClusterChunk >:: 00523 void ConeClusterAlgo <CalItem>:: 00524 getItemsInCone_bis(list<const CalItem*> &tlist, float etaJet, float phiJet, 00525 float cone_radius, float zvertex_in) const { 00526 // 00527 // provide the list of Items (towers, Cells...) containing the energy from a 00528 // jet of a given cone size 00529 // 00530 // WARNING: this is only to be used to compare to RUN I cone jets 00531 // 00532 float ZVERTEX_MAX=200.; 00533 float DMIN=80.; 00534 float DMAX=360.; 00535 float THETA_margin=0.022; 00536 00537 float zvertex=zvertex_in; 00538 float d1,d2; 00539 float phi_d1, phi_d2; 00540 float theta_E1, r1, r2, z1, z2; 00541 float theta_d1, theta_d2, eta_d1, eta_d2; 00542 00543 // Ignore very large vertex positions 00544 if (fabs(zvertex) > ZVERTEX_MAX ) zvertex=0.0; 00545 00546 if (zvertex >=0. ) { 00547 d1=fabs(DMIN-zvertex); 00548 d2=fabs(DMAX+zvertex); 00549 } else { 00550 d1=fabs(DMAX-zvertex); 00551 d2=fabs(DMIN+zvertex); 00552 } 00553 00554 // calculate theta of physics cone and find which eta's this intercepts 00555 // a the maximum points 00556 00557 phi_d1 = phiJet+cone_radius; 00558 //theta_E1 = kinem::theta(etaJet+cone_radius); 00559 theta_E1 = inline_maths::theta(etaJet+cone_radius); 00560 z1 = zvertex+d1*cos(theta_E1); 00561 r1 = d1*sin(theta_E1); 00562 00563 phi_d2 = phiJet-cone_radius; 00564 //theta_E1 = kinem::theta(etaJet-cone_radius); 00565 theta_E1 = inline_maths::theta(etaJet-cone_radius); 00566 z2 = zvertex+d2*cos(theta_E1); 00567 r2 = d2*sin(theta_E1); 00568 00569 // maximum spread in detector theta 00570 00571 theta_d1 = atan2(r1, z1); 00572 theta_d2 = atan2(r2, z2); 00573 00574 // make sure they stay in the calorimeter 00575 00576 theta_d1=max(theta_d1, THETA_margin); 00577 theta_d2=max(theta_d2, THETA_margin); 00578 //theta_d1=min(kinem::PI-(double)THETA_margin, (double)theta_d1); 00579 theta_d1=min(inline_maths::PI-(double)THETA_margin, (double)theta_d1); 00580 //theta_d2=min(kinem::PI-(double)THETA_margin, (double)theta_d2); 00581 theta_d2=min(inline_maths::PI-(double)THETA_margin, (double)theta_d2); 00582 00583 00584 //eta_d1 = kinem::eta(theta_d1); 00585 eta_d1 = inline_maths::eta(theta_d1); 00586 //eta_d2 = kinem::eta(theta_d2); 00587 eta_d2 = inline_maths::eta(theta_d2); 00588 00589 float signe; 00590 00591 if( eta_d1>=0.0 ) signe= 1.0; 00592 else signe= -1.0; 00593 int ietaMAX= eta_d1/0.1+signe; 00594 if(fabs(eta_d1)>=4.45) ietaMAX= 37*signe; 00595 else if(fabs(eta_d1)>=4.1) ietaMAX= 36*signe; 00596 else if(fabs(eta_d1)>=3.7) ietaMAX= 35*signe; 00597 else if(fabs(eta_d1)>=3.42) ietaMAX= 34*signe; 00598 else if(fabs(eta_d1)>=3.2) ietaMAX= 33*signe; 00599 00600 if( eta_d2>=0.0 ) signe= 1.0; 00601 else signe= -1.0; 00602 int ietaMIN= eta_d2/0.1+signe; 00603 if(fabs(eta_d2)>=4.45) ietaMIN= 37*signe; 00604 else if(fabs(eta_d2)>=4.1) ietaMIN= 36*signe; 00605 else if(fabs(eta_d2)>=3.7) ietaMIN= 35*signe; 00606 else if(fabs(eta_d2)>=3.42) ietaMIN= 34*signe; 00607 else if(fabs(eta_d2)>=3.2) ietaMIN= 33*signe; 00608 00609 //int iphiMAX= 64*phi_d1/(2.*kinem::PI)+1; 00610 int iphiMAX= 64*phi_d1/(2.*inline_maths::PI)+1; 00611 //int iphiMIN= 64*phi_d2/(2.*kinem::PI)+1; 00612 int iphiMIN= 64*phi_d2/(2.*inline_maths::PI)+1; 00613 00614 typename list<const CalItem*>::iterator it; 00615 for (it=tlist.begin() ; it != tlist.end() ; ) { 00616 //float eta_cur= (*it)->eta(); 00617 //float phi_cur= (*it)->phi(); 00618 int ieta= (*it)->address().ieta(); 00619 int iphi= (*it)->address().iphi(); 00620 00621 bool accepted = ieta<ietaMAX && ieta>ietaMIN; 00622 if ( iphiMIN>0 && iphiMAX<=64 ) { 00623 accepted = accepted && iphi<iphiMAX && iphi>iphiMIN; 00624 } 00625 else{ // case the cone overlap the phi=0=2pi line 00626 if ( iphiMIN>0 ){ 00627 accepted = accepted && 00628 ((iphi>iphiMIN && iphi<=64) || iphi<iphiMAX-64); 00629 } 00630 else{ 00631 accepted = accepted && 00632 ((iphi<iphiMAX && iphi>0) || iphi>iphiMIN+64); 00633 } 00634 } 00635 if ( ! accepted ) it = tlist.erase(it); 00636 else ++it; 00637 00638 } 00639 } 00640 ///////////////////////////////////////////////////////// 00641 //template < class CalItem,class CalItemAddress,class CalIClusterChunk > 00642 template < class CalItem > 00643 //inline void ConeClusterAlgo <CalItem,CalItemAddress,CalIClusterChunk >:: 00644 inline void ConeClusterAlgo <CalItem >:: 00645 print(ostream &os) const { 00646 os<<endl<<" CONE ALGORITHM, cone radius= "<<_CONErad<<endl<< 00647 " min E_T fraction= "<<_JETmne<<endl<< 00648 " minimum Delta_R separation between cones= "<<_TWOrad<<endl<< 00649 " shared E_T fraction threshold for combining jets= "<<_SPLifr<<endl; 00650 } 00651 ///////////////////////////////////////////////////////// 00652 00653 //template < class CalItem,class CalItemAddress,class CalIClusterChunk > 00654 template < class CalItem > 00655 //void ConeClusterAlgo <CalItem,CalItemAddress,CalIClusterChunk >:: 00656 void ConeClusterAlgo <CalItem >:: 00657 makeClusters(//const EnergyClusterReco* r, 00658 std::list<CalItem> &jets, 00659 list<const CalItem*> &itemlist, float Zvertex 00660 //, const EnergyClusterCollection<CalItemAddress> &preclu, 00661 //CalIClusterChunk* chunkptr 00662 //) const { 00663 ) { 00664 00665 // create an energy cluster collection for jets 00666 //EnergyClusterCollection<CalItemAddress>* ptrcol; 00667 //r->createClusterCollection(chunkptr, ptrcol); 00668 std::vector<const CalItem*> ecv; 00669 for ( typename std::list<const CalItem*>::iterator it = itemlist.begin(); 00670 it != itemlist.end(); it++) { 00671 ecv.push_back(*it); 00672 } 00673 00674 00675 // Initialize 00676 float Rcut= 1.E-06; 00677 if(_Increase_Delta_R) Rcut= 1.E-04; 00678 bool nojets; 00679 00680 //vector< TemporaryJet > TempColl; 00681 list< pair<float,float> > LTrack; 00682 00683 // get a vector with pointers to EnergyCluster in the collection 00684 //vector<const EnergyCluster<CalItemAddress>*> ecv; 00685 //preclu.getClusters(ecv); 00686 00687 // loop over all preclusters 00688 //typename vector<const EnergyCluster<CalItemAddress>*>::iterator jclu; 00689 typename std::vector<const CalItem*>::iterator jclu; 00690 for( jclu=ecv.begin(); jclu!=ecv.end(); jclu++ ) { 00691 ////const EnergyCluster<CalItemAddress>* ptr= *jclu; 00692 const CalItem* ptr= *jclu; 00693 //float PHIst= ptr->phi(); 00694 //float ETAst= ptr->eta(); 00695 float pz[4]; 00696 ptr->p4vec(pz); 00697 float ETAst= E2eta(pz); 00698 float PHIst= E2phi(pz); 00699 00700 //cout << "seed 4-vec "; 00701 //for ( int i = 0; i < 4; i++) cout << pz[i] << " "; 00702 //cout << endl; 00703 00704 nojets= false; 00705 // check to see if precluster is too close to a found jet 00706 if(_Kill_Far_Clusters) { 00707 list< pair<float,float> >::iterator kj; 00708 for(kj=LTrack.begin(); kj!=LTrack.end(); kj++) { 00709 if(DELTA_r((*kj).first,ETAst,(*kj).second,PHIst)<_Far_Ratio*_CONErad) { 00710 nojets= true; 00711 //cout << "seed too close ! skip " << endl; 00712 break; 00713 } 00714 } 00715 } 00716 if( nojets==false ) { 00717 TemporaryJet TJet(ETAst,PHIst); 00718 list< pair<int,float> > JETshare; 00719 00720 // start of cone building loop 00721 int trial= 0; 00722 do{ 00723 trial++; 00724 //cout << " trial " << trial << endl; 00725 00726 ETAst= TJet.Eta(); 00727 PHIst= TJet.Phi(); 00728 TJet.erase(); 00729 00730 //if(PHIst > kinem::TWOPI) PHIst= PHIst-kinem::TWOPI; 00731 if(PHIst > inline_maths::TWOPI) PHIst= PHIst-inline_maths::TWOPI; 00732 //if(PHIst < 0.0 ) PHIst= PHIst+kinem::TWOPI; 00733 if(PHIst < 0.0 ) PHIst= PHIst+inline_maths::TWOPI; 00734 //if( PHIst>kinem::TWOPI || PHIst<0.0 ) { 00735 if( PHIst>inline_maths::TWOPI || PHIst<0.0 ) { 00736 TJet.setEtaPhiEt(0.0,0.0,0.0); 00737 break; // end loop do (illegal jet PHI) 00738 } 00739 TJet.setEtaPhiEt(ETAst,PHIst,0.0); 00740 00741 // calculate eta & phi limits for cone 00742 list<const CalItem*> Twlist(itemlist); 00743 00744 getItemsInCone(Twlist,ETAst,PHIst,_CONErad,Zvertex); 00745 // only to compare with RUN I cone jets ! getItemsInCone_bis(Twlist,ETAst,PHIst,_CONErad,Zvertex); 00746 00747 // loop over the possible items for this cone 00748 typename list<const CalItem*>::iterator tk; 00749 for( tk= Twlist.begin(); tk!=Twlist.end(); tk++ ) { 00750 float ETk =(*tk)->pT(); 00751 // now done in CalCell/CalTower if((*tk)->E() < 0.0) ETk= -ETk; 00752 00753 if( ETk > _Eitem_Negdrop ) { 00754 //float ETAk=(*tk)->eta(); 00755 //float PHIk=(*tk)->phi(); 00756 float pz[4]; 00757 (*tk)->p4vec(pz); 00758 float ETAk= E2eta(pz); 00759 float PHIk= E2phi(pz); 00760 00761 float dphi= fabs(PHIk-PHIst); 00762 //if(dphi > kinem::TWOPI-dphi) { 00763 if(dphi > inline_maths::TWOPI-dphi) { 00764 //if(PHIst < PHIk) PHIk= PHIk-kinem::TWOPI; 00765 if(PHIst < PHIk) PHIk= PHIk-inline_maths::TWOPI; 00766 //else PHIk= PHIk+kinem::TWOPI; 00767 else PHIk= PHIk+inline_maths::TWOPI; 00768 } 00769 00770 if( R2_bis(ETAk,PHIk,ETAst,PHIst) <= _CONErad*_CONErad ) { 00771 TJet.addItem(*tk); 00772 } 00773 } 00774 }// end loop tk 00775 00776 if(TJet.updateEtaPhiEt()==false) { 00777 //cout << " negative E jet ! drop " << endl; 00778 break; 00779 } 00780 00781 // require some minimum ET on every iteration 00782 if(_Jet_Et_Min_On_Iter) { 00783 if( TJet.Et() < _JETmne*_Et_Min_Ratio ) { 00784 //cout << " too low ET jet ! drop " << endl; 00785 break; // end loop trial 00786 } 00787 } 00788 00789 //cout << " R2 = " << R2_bis(TJet.Eta(),TJet.Phi(),ETAst,PHIst) << 00790 // " Rcut = " << Rcut << endl; 00791 }while(R2_bis(TJet.Eta(),TJet.Phi(),ETAst,PHIst)>=Rcut && trial<=50); 00792 00793 if( TJet.Et() >= _JETmne ) { 00794 //cout << " jet accepted will check for overlaps " << endl; 00795 if(_D0_Angle) TJet.D0_Angle_updateEtaPhi(); 00796 //cout << " after TJet.D0_Angle_updateEtaPhi() " << endl; 00797 00798 // item also in another jet 00799 list<const CalItem*> Lst; 00800 TJet.getItems(Lst); 00801 typename list<const CalItem*>::iterator tk; 00802 for(tk=Lst.begin(); tk!=Lst.end(); tk++) { 00803 float ETk=(*tk)->pT(); 00804 // now done in CalCell/CalTower if((*tk)->E() < 0.0) ETk= -ETk; 00805 for(unsigned int kj=0; kj<TempColl.size(); kj++) { 00806 if(TempColl[kj].ItemInJet(*tk)==true) { 00807 list< pair<int,float> >::iterator pit; 00808 bool jetok= false; 00809 for(pit=JETshare.begin(); pit!=JETshare.end();pit++) { 00810 if((*pit).first == (int) kj) { 00811 jetok= true; 00812 (*pit).second+= ETk; 00813 break; 00814 } 00815 } 00816 if(jetok==false) JETshare.push_back(make_pair(kj,ETk)); 00817 } 00818 } 00819 } 00820 00821 if(JETshare.size() >0) { 00822 list< pair<int,float> >::iterator pit; 00823 float Ssum= 0.0; 00824 list< pair<int,float> >::iterator pitMAX=JETshare.begin(); 00825 for(pit=JETshare.begin(); pit!=JETshare.end(); pit++) { 00826 Ssum+= (*pit).second; 00827 if((*pit).second > (*pitMAX).second) pitMAX= pit; 00828 } 00829 00830 //int IJET= (*pitMAX).first; 00831 bool splshr; 00832 float Eleft= fabs(TJet.Et()-Ssum); 00833 float Djets= TempColl[(*pitMAX).first].dist_R2(TJet); 00834 if(Djets <= _TWOrad || Eleft <= _Thresh_Diff_Et) { 00835 TJet.erase(); 00836 splshr= false; 00837 } 00838 else { 00839 float SharedFr=Ssum/min(TempColl[(*pitMAX).first].Et(),TJet.Et()); 00840 if(JETshare.size() >1) { 00841 typename list<const CalItem*>::iterator tk; 00842 for(tk=TJet.LItems().begin(); tk!=TJet.LItems().end(); ) { 00843 bool found = false; 00844 list< pair<int,float> >::iterator pit; 00845 for(pit=JETshare.begin(); pit!=JETshare.end();pit++) { 00846 if((*pit).first!=(*pitMAX).first) { 00847 if(TempColl[(*pit).first].ItemInJet(*tk)==true) { 00848 tk = TJet.LItems().erase(tk); 00849 found = true; 00850 break; 00851 } 00852 } 00853 } 00854 if ( !found ) ++tk; 00855 } 00856 } 00857 00858 splshr= TempColl[(*pitMAX).first].share_jets(TJet,SharedFr,_SPLifr); 00859 00860 if( splshr==true ) { 00861 //cout << " jet splitted due to overlaps " << endl; 00862 TempColl[(*pitMAX).first].updateEtaPhiEt(); 00863 TJet.updateEtaPhiEt(); 00864 if(_D0_Angle) TJet.D0_Angle_updateEtaPhi(); 00865 if(_D0_Angle) TempColl[(*pitMAX).first].D0_Angle_updateEtaPhi(); 00866 TempColl.push_back(TJet); 00867 LTrack.push_back(make_pair(TJet.Eta(),TJet.Phi())); 00868 } 00869 else { 00870 //cout << " jet merged due to overlaps " << endl; 00871 TempColl[(*pitMAX).first].updateEtaPhiEt(); 00872 if(_D0_Angle) TempColl[(*pitMAX).first].D0_Angle_updateEtaPhi(); 00873 } 00874 } 00875 } 00876 else { 00877 TJet.updateEtaPhiEt(); 00878 if(_D0_Angle) TJet.D0_Angle_updateEtaPhi(); 00879 TempColl.push_back(TJet); 00880 LTrack.push_back(make_pair(TJet.Eta(),TJet.Phi())); 00881 } 00882 } //JETmne 00883 } //nojets 00884 }// end loop jclu 00885 00886 for(unsigned int i=0; i<TempColl.size(); i++) { 00887 //EnergyCluster<CalItemAddress>* ptrclu; 00888 CalItem ptrclu; 00889 //r->createCluster(ptrcol,ptrclu); 00890 list<const CalItem*> Vi; 00891 TempColl[i].getItems(Vi); 00892 typename list<const CalItem*>::iterator it; 00893 for(it=Vi.begin(); it!=Vi.end(); it++) { 00894 const CalItem* ptr= *it; 00895 //CalItemAddress addr= ptr->address(); 00896 float p[4]; 00897 ptr->p4vec(p); 00898 //float emE= ptr->emE(); 00899 //r->addClusterItem(ptrclu,addr,p,emE); 00900 ptrclu.Add(*ptr); 00901 } 00902 jets.push_back(ptrclu); 00903 } 00904 00905 }// end 00906 00907 } //namespace d0runi 00908 00909 FASTJET_END_NAMESPACE 00910 00911 #endif // CONECLUSTERALGO_H 00912 00913 00914