Go to the documentation of this file.00001 #ifndef FORIA_UTILS_HH
00002 #define FORIA_UTILS_HH
00003 #include "ForIA/IFourMomentum.hh"
00004 #include <cmath>
00005
00006 namespace ForIA{
00007
00008 const double TWOPI = 2.0 * M_PI;
00009
00011 inline double modMPi2Pi(double angle){
00012
00013 if(angle < M_PI && angle > -M_PI) return angle;
00014
00015 double dPhi = fmod(angle + M_PI, TWOPI);
00016 return (angle > 0)? dPhi - M_PI: dPhi + M_PI;
00017 }
00018
00020 inline double mod2Pi(double angle){
00021 if(angle < TWOPI && angle > 0) return angle;
00022
00023 double dPhi = fmod(angle, TWOPI);
00024 return (angle >0)? dPhi: TWOPI + dPhi;
00025 }
00026
00028 inline double modPi(double angle){
00029 return fabs(modMPi2Pi(angle));
00030 }
00031
00033 inline double PT2( IFourMomentumConstPtr p1, IFourMomentumConstPtr p2, double smear1, double smear2) {
00034
00035 double px = p1->px()*smear1 + p2->px()*smear2;
00036 double py = p1->py()*smear1 + p2->py()*smear2;
00037
00038 return px*px + py*py;
00039 }
00040
00042 inline double PT2( IFourMomentumConstPtr p1, IFourMomentumConstPtr p2) {
00043
00044 double px = p1->px() + p2->px();
00045 double py = p1->py() + p2->py();
00046
00047 return px*px + py*px;
00048 }
00049
00051 inline double PT(IFourMomentumConstPtr p1, IFourMomentumConstPtr p2) {
00052 return sqrt(PT2(p1, p2));
00053 }
00054
00056 inline double PT( IFourMomentumConstPtr p1, IFourMomentumConstPtr p2, double smear1, double smear2) {
00057 return sqrt(PT2(p1, p2, smear1, smear2));
00058 }
00059
00061 inline double P2( IFourMomentumConstPtr p1, IFourMomentumConstPtr p2) {
00062
00063 double pz = p1->pz() + p2->pz();
00064 return PT2(p1, p2) + pz*pz;
00065 }
00066
00068 inline double P2( IFourMomentumConstPtr p1, IFourMomentumConstPtr p2, double smear1, double smear2) {
00069 double pz = p1->pz()*smear1 + p2->pz()*smear2;
00070 return PT2(p1, p2, smear1, smear2) + pz*pz;
00071 }
00072
00074 inline double P( IFourMomentumConstPtr p1, IFourMomentumConstPtr p2) {
00075 return sqrt(P2(p1, p2));
00076 }
00077
00079 inline double invMass( IFourMomentumConstPtr p1, IFourMomentumConstPtr p2) {
00080
00081 double e = p1->E() + p2->E();
00082
00083 return sqrt(e*e - P2(p1, p2));
00084 }
00085
00087 inline double invMass( IFourMomentumConstPtr p1, IFourMomentumConstPtr p2, double smear1, double smear2) {
00088
00089 double e = p1->E()*smear1 + p2->E()*smear2;
00090 return sqrt(e*e - P2(p1, p2, smear1, smear2));
00091 }
00092
00093 inline double deltaPhi( IFourMomentumConstPtr p1, IFourMomentumConstPtr p2) {
00094 double dPhi = p1->phi() - p2->phi();
00095 return modPi(dPhi);
00096
00097
00098
00099
00100
00101
00102
00103 }
00104
00105 inline double deltaR( IFourMomentumConstPtr p1, IFourMomentumConstPtr p2) {
00106 double dPhi = deltaPhi(p1,p2);
00107 double dEta = p1->eta() - p2->eta();
00108 return sqrt (dPhi*dPhi + dEta*dEta);
00109
00110 }
00111
00112 template<class T>
00113 inline double trackWidth(IFourMomentumConstPtr jet, const vector<boost::shared_ptr<T const> > &tracks){
00114
00115 double sumPT = 0.;
00116 double sumPTR = 0.;
00117 for(typename vector<boost::shared_ptr<T const> >::const_iterator track = tracks.begin();
00118 track != tracks.end(); ++track){
00119 sumPT += (*track)->PT();
00120 sumPTR += (*track)->PT() * deltaR(*track, jet);
00121 }
00122 return sumPTR / sumPT;
00123 }
00124
00125 inline double UtilsScaleNominal (float eta) {
00126
00127 if ( eta >= -4.8 && eta < -4.2) return 1.04;
00128 else if ( eta >= -4.2 && eta < -3.5) return 0.983;
00129 else if ( eta >= -3.5 && eta < -3.2) return 1.01;
00130 else if ( eta >= -3.2 && eta < -2.8) return 0.973;
00131 else if ( eta >= -2.8 && eta < -2.37) return 0.911;
00132
00133
00134 else if ( eta >= -2.37 && eta < -1.52) return 0.978;
00135 else if ( eta >= -1.52 && eta < -1.37) return 0.927;
00136 else if ( eta >= -1.37 && eta < 0) return 0.983;
00137 else if ( eta >= 0 && eta < 1.37) return 0.987;
00138 else if ( eta >= 1.37 && eta < 1.52) return 0.987;
00139 else if ( eta >= 1.52 && eta < 2.37) return 0.969;
00140
00141
00142 else if ( eta >= 2.37 && eta < 2.8) return 0.893;
00143 else if ( eta >= 2.8 && eta < 3.2) return 0.946;
00144 else if ( eta >= 3.2 && eta < 3.5) return 1.04;
00145 else if ( eta >= 3.5 && eta < 4.2) return 0.958;
00146 else if ( eta >= 4.2 && eta < 4.8) return 1.01;
00147
00148 return 1.0;
00149
00150 }
00151
00152
00153 inline double UtilsScaleUp (float eta) {
00154
00155 if ( eta >= -4.8 && eta < -4.2) return 1.071;
00156 else if ( eta >= -4.2 && eta < -3.5) return 1.081;
00157 else if ( eta >= -3.5 && eta < -3.2) return 1.132;
00158 else if ( eta >= -3.2 && eta < -2.8) return 1.052;
00159 else if ( eta >= -2.8 && eta < -2.37) return 1.059;
00160 else if ( eta >= -2.37 && eta < -1.52) return 1.043;
00161 else if ( eta >= -1.52 && eta < -1.37) return 1.091;
00162 else if ( eta >= -1.37 && eta < -0.8) return 1.044;
00163 else if ( eta >= -0.8 && eta < 0) return 1.033;
00164 else if ( eta >= 0 && eta < 0.8) return 1.033;
00165 else if ( eta >= 0.8 && eta < 1.37) return 1.044;
00166 else if ( eta >= 1.37 && eta < 1.52) return 1.086;
00167 else if ( eta >= 1.52 && eta < 2.37) return 1.043;
00168 else if ( eta >= 2.37 && eta < 2.8) return 1.06;
00169 else if ( eta >= 2.8 && eta < 3.2) return 1.053;
00170 else if ( eta >= 3.2 && eta < 3.5) return 1.123;
00171 else if ( eta >= 3.5 && eta < 4.2) return 1.083;
00172 else if ( eta >= 4.2 && eta < 4.8) return 1.073;
00173
00174
00175 return 1.0;
00176
00177 }
00178
00179 inline double UtilsScaleDown (float eta) {
00180
00181 if ( eta >= -4.8 && eta < -4.2) return 0.9776;
00182 else if ( eta >= -4.2 && eta < -3.5) return 0.967;
00183 else if ( eta >= -3.5 && eta < -3.2) return 0.903;
00184 else if ( eta >= -3.2 && eta < -2.8) return 0.9761;
00185 else if ( eta >= -2.8 && eta < -2.37) return 0.9727;
00186 else if ( eta >= -2.37 && eta < -1.52) return 0.957;
00187 else if ( eta >= -1.52 && eta < -1.37) return 0.9097;
00188 else if ( eta >= -1.37 && eta < -0.8) return 0.9545;
00189 else if ( eta >= -0.8 && eta < 0) return 0.9656;
00190 else if ( eta >= 0 && eta < 0.8) return 0.9657;
00191 else if ( eta >= 0.8 && eta < 1.37) return 0.9546;
00192 else if ( eta >= 1.37 && eta < 1.52) return 0.9152;
00193 else if ( eta >= 1.52 && eta < 2.37) return 0.9565;
00194 else if ( eta >= 2.37 && eta < 2.8) return 0.9725;
00195 else if ( eta >= 2.8 && eta < 3.2) return 0.9758;
00196 else if ( eta >= 3.2 && eta < 3.5) return 0.9119;
00197 else if ( eta >= 3.5 && eta < 4.2) return 0.9661;
00198 else if ( eta >= 4.2 && eta < 4.8) return 0.977;
00199
00200 return 1.0;
00201
00202 }
00203
00204
00205 inline double UtilsPi0Up (float eta) {
00206
00207 if ( eta >= -4.8 && eta < -4.2) return 1.023;
00208 else if ( eta >= -4.2 && eta < -3.5) return 1.034;
00209 else if ( eta >= -3.5 && eta < -3.2) return 1.11;
00210 else if ( eta >= -3.2 && eta < -2.8) return 1.025;
00211 else if ( eta >= -2.8 && eta < -2.37) return 1.029;
00212 else if ( eta >= -2.37 && eta < -1.52) return 1.02;
00213 else if ( eta >= -1.52 && eta < -1.37) return 1.178;
00214 else if ( eta >= -1.37 && eta < 0) return 1.025;
00215 else if ( eta >= 0 && eta < 1.37) return 1.025;
00216 else if ( eta >= 1.37 && eta < 1.52) return 1.178;
00217 else if ( eta >= 1.52 && eta < 2.37) return 1.02;
00218 else if ( eta >= 2.37 && eta < 2.8) return 1.029;
00219 else if ( eta >= 2.8 && eta < 3.2) return 1.024;
00220 else if ( eta >= 3.2 && eta < 3.5) return 1.104;
00221 else if ( eta >= 3.5 && eta < 4.2) return 1.034;
00222 else if ( eta >= 4.2 && eta < 4.8) return 1.023;
00223 }
00224
00225
00226 inline double UtilsPi0Down (float eta) {
00227
00228
00229 if ( eta >= -4.8 && eta < -4.2) return 0.977;
00230 else if ( eta >= -4.2 && eta < -3.5) return 0.968;
00231 else if ( eta >= -3.5 && eta < -3.2) return 0.902;
00232 else if ( eta >= -3.2 && eta < -2.8) return 0.977;
00233 else if ( eta >= -2.8 && eta < -2.37) return 0.975;
00234 else if ( eta >= -2.37 && eta < -1.52) return 0.979;
00235 else if ( eta >= -1.52 && eta < -1.37) return 0.825;
00236 else if ( eta >= -1.37 && eta < 0) return 0.969;
00237 else if ( eta >= 0 && eta < 1.37) return 0.969;
00238 else if ( eta >= 1.37 && eta < 1.52) return 0.825;
00239 else if ( eta >= 1.52 && eta < 2.37) return 0.979;
00240 else if ( eta >= 2.37 && eta < 2.8) return 0.975;
00241 else if ( eta >= 2.8 && eta < 3.2) return 0.977;
00242 else if ( eta >= 3.2 && eta < 3.5) return 0.908;
00243 else if ( eta >= 3.5 && eta < 4.2) return 0.968;
00244 else if ( eta >= 4.2 && eta < 4.8) return 0.977;
00245
00246
00247
00248
00249 }
00250
00251 inline double UtilsScaleMC (float eta, std::string scaleType) {
00252
00253 if(scaleType=="NONE") return 1.0;
00254 if(scaleType=="NOMINAL") return UtilsScaleNominal(eta);
00255 if(scaleType=="UP") return UtilsScaleNominal(eta)*UtilsScaleUp(eta);
00256 if(scaleType=="DOWN") return UtilsScaleNominal(eta)*UtilsScaleDown(eta);
00257
00258 return 1.;
00259 }
00260
00261 inline double UtilsClusterCellSyst (float eta, float p) {
00262
00263
00264 float feta = fabs(eta);
00265
00266 if (feta < 0.6) {
00267 if (p < 2000) return 0.05;
00268 else if (p < 4000) return 0.02;
00269 else return 0.01;
00270 }
00271 else if (feta < 1.1) {
00272 if (p < 1200) return 0.00;
00273 else if (p < 4000) return 0.05;
00274 else if (p < 10000) return 0.02;
00275 else return 0.0;
00276 }
00277 else if (feta < 1.4) {
00278 if (p < 1200) return 0.00;
00279 else if (p < 6000) return 0.05;
00280 else return 0.01;
00281 }
00282 else if (feta < 1.5) {
00283 if (p < 1800) return -0.1;
00284 else if (p < 2200) return -0.05;
00285 else return 0.02;
00286 }
00287 else if (feta < 1.8) {
00288 if (p < 1800) return -0.2;
00289 else if (p < 5000) return -0.1;
00290 else return 0.05;
00291 }
00292 else if (feta < 1.9) {
00293 if (p < 4000) return -0.1;
00294 else if (p < 6000) return -0.02;
00295 else return 0.05;
00296 }
00297 else if (feta < 2.4) {
00298 if (p < 3000) return -0.05;
00299 else return 0.05;
00300 }
00301 return 0.0;
00302
00303 }
00304
00305 inline double UtilsEPSyst (float eta, float p) {
00306
00307
00308 float feta = fabs(eta);
00309 if (feta < 0.6) {
00310 if (p < 800) return 0.1;
00311 else return 0.02;
00312 }
00313 else if (feta < 1.1) {
00314
00315 return 0.02;
00316 }
00317 else if (feta < 1.4) {
00318 if (p < 800) return 0.6;
00319 else return 0.01;
00320 }
00321 else if (feta < 1.5) {
00322 return 0.04;
00323 }
00324 else if (feta < 1.8) {
00325 if (p < 1800) return 0.05;
00326 else return 0.02;
00327 }
00328 else if (feta < 1.9) {
00329 if (p < 1800) return 0.025;
00330 else return 0.05;
00331 }
00332 else if (feta < 2.4) {
00333 if (p < 1800) return 0.025;
00334 else if (p < 12000) return 0.075;
00335 else return 0.0;
00336 }
00337
00338 }
00339
00340
00341 }
00342
00343 #endif
