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