GlobalCaloTrigger/test/metTest.cpp

0001 #include <iostream>
0002 #include <cmath>
0003
0004 using std::cout;
0005 using std::endl;
0006
0007 struct etmiss_vec {
0008       unsigned mag;
0009       unsigned phi;
0010 };
0011 etmiss_vec calculate_etmiss_vec (const int ex, const int ey) ;
0012 unsigned phFine(unsigned ph) {
0013   unsigned r = ph % 18;
0014   return (r<9 ? r : 17-r);
0015 }
0016
0017 int main(int argc, char **argv)
0018 {
0019   const unsigned maxbin=9;
0020   unsigned histTru[maxbin]={0,0,0,0,0,0,0,0,0};
0021   unsigned histGct[maxbin]={0,0,0,0,0,0,0,0,0};
0022   for (int ex=-50; ex<50; ex++) { for (int ey=-50; ey<50; ey++) {
0023     //if ((ex*ex+ey*ey)<100) {
0024     if ((ex*ex+ey*ey)<1000 && (ex*ex+ey*ey)>=500) {
0025       etmiss_vec met = calculate_etmiss_vec(ex,ey);
0026       unsigned phiF = phFine(met.phi);
0027       if (phiF<=maxbin) { histGct[phiF]++; } else { cout << " unexpected phi " << phiF << endl; }
0028       float phiVal = (atan2( (float) ey, (float) ex))/M_PI;
0029       if (phiVal<0) { phiVal += 2.0; }
0030       unsigned phiTrue = ((unsigned) (36.*phiVal));
0031       unsigned phiT = phFine(phiTrue);
0032       if (phiT<=maxbin) { histTru[phiT]++; } else { cout << " unexpected phi " << phiT << endl; }
0033       if (met.phi != phiTrue) {
0034         cout << "ex " << ex << " ey " << ey << " gct phi " << met.phi
0035                                             << " tru phi " << phiTrue << endl;
0036       }
0037     }
0038   } }
0039   cout << "Phi bins Gct"; for (int b=0; b<maxbin; b++) { cout << "  " << histGct[b]; } cout << endl;
0040   cout << "Phi bins Tru"; for (int b=0; b<maxbin; b++) { cout << "  " << histTru[b]; } cout << endl;
0041   return 0;
0042 }
0043
0044 etmiss_vec calculate_etmiss_vec (const int ex, const int ey)
0045 {
0046   //---------------------------------------------------------------------------------
0047   //
0048   // Calculates magnitude and direction of missing Et, given measured Ex and Ey.
0049   //
0050   // The algorithm used is suitable for implementation in hardware, using integer
0051   // multiplication, addition and comparison and bit shifting operations.
0052   //
0053   // Proceed in two stages. The first stage gives a result that lies between
0054   // 92% and 100% of the true Et, with the direction measured in 45 degree bins.
0055   // The final precision depends on the number of factors used in corrFact.
0056   // The present version with eleven factors gives a precision of 1% on Et, and
0057   // finds the direction to the nearest 5 degrees.
0058   //
0059   //---------------------------------------------------------------------------------
0060   etmiss_vec result;
0061
0062   unsigned eneCoarse, phiCoarse;
0063   unsigned eneCorect, phiCorect;
0064
0065   const unsigned root2fact = 181;
0066   const unsigned corrFact[11] = {24, 39, 51, 60, 69, 77, 83, 89, 95, 101, 106};
0067   const unsigned corrDphi[11] = { 0,  1,  1,  2,  2,  3,  3,  3,  3,   4,   4};
0068
0069   bool s[3];
0070   unsigned Mx, My, Mw;
0071
0072   unsigned Dx, Dy;
0073   unsigned Fx, Fy;
0074   unsigned eFact;
0075
0076   unsigned b,phibin;
0077   bool midphi=false;
0078
0079   // Here's the coarse calculation, with just one multiply operation
0080   //
0081   My = (ey<0 ? -ey : ey);
0082   Mx = (ex<0 ? -ex : ex);
0083   Mw = (((Mx+My)*root2fact)+0x80)>>8;
0084
0085   s[0] = (ey<0);
0086   s[1] = (ex<0);
0087   s[2] = (My>Mx);
0088
0089   phibin = 0; b = 0;
0090   for (int i=0; i<3; i++) {
0091     if (s[i]) { b=1-b;} phibin = 2*phibin + b;
0092   }
0093
0094   unsigned m=(My>Mx ? My : Mx);
0095   eneCoarse = (Mw>m ? Mw : m );
0096   phiCoarse = phibin*9;
0097
0098   // For the fine calculation we multiply both input components
0099   // by all the factors in the corrFact list in order to find
0100   // the required corrections to the energy and angle
0101   //
0102   for (eFact=0; eFact<10; eFact++) {
0103     Fx = Mx;
0104     Fy = My;
0105     Dx = (Mx*corrFact[eFact])>>8;
0106     Dy = (My*corrFact[eFact])>>8;
0107     if         ((Dx>=Fy) || (Dy>=Fx))         {midphi=false; break;}
0108     if ((Fx+Dx)>=(Fy-Dy) && (Fy+Dy)>=(Fx-Dx)) {midphi=true;  break;}
0109   }
0110   eneCorect = (eneCoarse*(128+eFact))>>7;
0111   if (midphi ^ (b==1)) {
0112     phiCorect = phiCoarse + 8 - corrDphi[eFact];
0113   } else {
0114     phiCorect = phiCoarse + corrDphi[eFact];
0115   }
0116
0117   // Store the result of the calculation
0118   //
0119   result.mag = eneCorect;
0120   result.phi = phiCorect;
0121
0122   return result;
0123 }