00001 // emacs this is -*-C++-*- 00002 #ifndef JETUTILS_JETKINEMATICS_H 00003 #define JETUTILS_JETKINEMATICS_H 00004 00005 // #include "EventKernel/I4Momentum.h" 00006 00007 // //#include "FourMom/P4Help.h" 00008 // #include "FourMomUtils/P4Helpers.h" 00009 00010 // #include "JetEvent/Jet.h" 00011 // #include "JetUtils/JetDistances.h" 00012 00013 // #include <cmath> 00014 00015 // struct JetKinematics 00016 // { 00017 00018 // typedef CLHEP::Hep3Vector vector_t; 00019 00020 // static double dotProduct(const I4Momentum* pJet1, const I4Momentum* pJet2) 00021 // { 00022 // return 00023 // vector_t(pJet1->px(),pJet1->py(),pJet1->pz()).dot(vector_t(pJet2->px(), 00024 // pJet2->py(), 00025 // pJet2->pz())); 00026 // }; 00027 00028 // static double dotProductTrans(const I4Momentum* pJet1, 00029 // const I4Momentum* pJet2) 00030 // { 00031 // return 00032 // vector_t(pJet1->px(),pJet1->py(),double(0.)).dot(vector_t(pJet2->px(), 00033 // pJet2->py(), 00034 // double(0.))); 00035 // }; 00036 00037 // static vector_t crossProduct(const I4Momentum* pJet1, 00038 // const I4Momentum* pJet2) 00039 // { 00040 // return 00041 // vector_t(pJet1->px(),pJet1->py(),pJet1->pz()).cross(vector_t(pJet2->px(), 00042 // pJet2->py(), 00043 // pJet2-> 00044 // pz())); 00045 // }; 00046 00047 // static CLHEP::Hep3Vector crossProductTrans(const I4Momentum* pJet1, 00048 // const I4Momentum* pJet2) 00049 // { 00050 // return 00051 // vector_t(pJet1->px(),pJet1->py(),double(0.)).cross(vector_t(pJet2->px(), 00052 // pJet2->py(), 00053 // double(0.))); 00054 00055 // }; 00056 00057 // static double projPtContrib(const I4Momentum* pJet1, 00058 // const I4Momentum* pJet2) 00059 // { 00060 // if ( pJet1->pt() == 0 || pJet2->pt() == 0. ) return 0.; 00061 // // calculate projection (absolute) 00062 // return dotProductTrans(pJet1,pJet2)/pJet1->pt(); 00063 // }; 00064 00065 // static double perpPtContrib(const I4Momentum* pJet1, 00066 // const I4Momentum* pJet2) 00067 // { 00068 // if ( pJet1->pt() == 0 || pJet2->pt() == 0. ) return 0.; 00069 // // calculate perpendicular component (absolute) 00070 // return (crossProductTrans(pJet1,pJet2).mag())/pJet1->pt(); 00071 // }; 00072 00073 // static double integratedProjPt(const Jet* pJet, 00074 // const double& radius) 00075 // { 00076 // Jet::constituent_iterator fCon(pJet->begin()); 00077 // Jet::constituent_iterator lCon(pJet->end()); 00078 // double intPtProj(double(0.0)); 00079 // for ( ; fCon != lCon; fCon++ ) 00080 // { 00081 // if ( P4Helpers::deltaR(pJet,(*fCon)) < radius ) 00082 // { 00083 // intPtProj += projPtContrib(pJet,(*fCon)); 00084 // } 00085 // } 00086 // return intPtProj; 00087 // }; 00088 00089 00090 00091 00092 // template<class CONT> 00093 // static double ptWeightedWidth(CONT & c, const I4Momentum* axis){ 00094 // typename CONT::const_iterator it = c.begin(); 00095 // typename CONT::const_iterator itE = c.end(); 00096 // return ptWeightedWidth(it, itE, axis); 00097 // } 00098 00099 // template<class ITERATOR> 00100 // static double ptWeightedWidth(ITERATOR & it, ITERATOR & itE, const I4Momentum* axis){ 00101 // double jetPtSumConst = 0; 00102 // double jetWidth = 0; 00103 00104 // for( ; it!=itE; ++it) 00105 // { 00106 // double dR = JetDistances::deltaR(axis,(*it)); 00107 // double pt=(*it)->pt(); 00108 // jetWidth += dR * pt; 00109 // jetPtSumConst += pt; 00110 // } 00111 00112 // if(jetPtSumConst <= 0) 00113 // return -1; 00114 00115 // jetWidth = jetWidth / jetPtSumConst; 00116 // return jetWidth; 00117 // } 00118 00119 00120 // /// Functions to compute various kinematic quantities related to a collection of 00121 // /// 4Momentum in a single call. 00122 // /// As a compromise for a stable yet flexible interface, we simply use a vector<double> as 00123 // /// the returned object with an enum helper to flag each value in this vector. 00124 // struct CollectionQuantities { 00125 // enum { 00126 // N, SumPt, WeightedWidth, 00127 // Px,Py,Pz,E, 00128 // MAX 00129 // }; 00130 // Jet::hlv_t toHepLotentzVect(std::vector<double> & v){ return Jet::hlv_t(v[Px], v[Py], v[Pz], v[E]); } 00131 // }; 00132 00133 // template<class CONT> 00134 // static std::vector<double> computeCollectionQuantities(CONT &c,const I4Momentum* parent){ 00135 // typename CONT::const_iterator it = c.begin(); 00136 // typename CONT::const_iterator itE = c.end(); 00137 // return computeCollectionQuantities(it, itE, parent); 00138 // } 00139 // template<class ITERATOR> 00140 // static std::vector<double> computeCollectionQuantities(ITERATOR & it, ITERATOR & itE, const I4Momentum* parent){ 00141 // std::vector<double> outV(CollectionQuantities::MAX,0.0); 00142 // typedef CollectionQuantities CQ; 00143 // for( ; it!=itE; ++it) 00144 // { 00145 // double dR = JetDistances::deltaR(parent,(*it)); 00146 // double pt=(*it)->pt(); 00147 // outV[CQ::WeightedWidth] += dR * pt; 00148 // outV[CQ::SumPt] += pt; 00149 // outV[CQ::N]+=1; 00150 // outV[CQ::Px] += (*it)->px();outV[CQ::Py] += (*it)->py();outV[CQ::Pz] += (*it)->pz();outV[CQ::E] += (*it)->e(); 00151 // } 00152 00153 // if(outV[CQ::SumPt] > 0){ 00154 // outV[CQ::WeightedWidth] = outV[CQ::WeightedWidth] / outV[CQ::SumPt] ; 00155 // } else outV[CQ::WeightedWidth] = -1; 00156 00157 // return outV; 00158 // } 00159 00160 00161 }; 00162 #endif