FastJet 3.0.2
ConeClusterAlgo.hpp
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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends