/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/TopQuarkAnalysis/TopPairBSM/src/METzCalculator.cc

Go to the documentation of this file.

#include "METzCalculator.h"
#include "TMath.h"

METzCalculator::METzCalculator() {
        isComplex_ = false;
}

METzCalculator::~METzCalculator() {
}

double
METzCalculator::Calculate(int type) {

        double M_W  = 80.4;
        double M_mu =  0.10566;
        double emu = lepton_.energy();
        double pxmu = lepton_.px();
        double pymu = lepton_.py();
        double pzmu = lepton_.pz();
        double pxnu = MET_.px();
        double pynu = MET_.py();
        double pznu = 0.;

                
        double a = M_W*M_W - M_mu*M_mu + 2.0*pxmu*pxnu + 2.0*pymu*pynu;
        double A = 4.0*(emu*emu - pzmu*pzmu);
        double B = -4.0*a*pzmu;
        double C = 4.0*emu*emu*(pxnu*pxnu + pynu*pynu) - a*a;

        double tmproot = B*B - 4.0*A*C;

        if (tmproot<0) {
                isComplex_= true;
                pznu = - B/(2*A); // take real part of complex roots

                //std::cout << " Neutrino Solutions: complex, real part " << pznu << std::endl;
                }
        else {
                        isComplex_ = false;
                        double tmpsol1 = (-B + TMath::Sqrt(tmproot))/(2.0*A);
                        double tmpsol2 = (-B - TMath::Sqrt(tmproot))/(2.0*A);

                        //std::cout << " Neutrino Solutions: " << tmpsol1 << ", " << tmpsol2 << std::endl;
                        
                        if (type == 0 ) {
                                // two real roots, pick the one closest to pz of muon
                                if (TMath::Abs(tmpsol2-pzmu) < TMath::Abs(tmpsol1-pzmu)) { pznu = tmpsol2;}
                                else pznu = tmpsol1;
                                // if pznu is > 300 pick the most central root
                                if ( pznu > 300. ) {
                                        if (TMath::Abs(tmpsol1)<TMath::Abs(tmpsol2) ) pznu = tmpsol1;
                                        else pznu = tmpsol2;
                                }
                        }
                        if (type == 1 ) {
                                // two real roots, pick the one closest to pz of muon
                                if (TMath::Abs(tmpsol2-pzmu) < TMath::Abs(tmpsol1-pzmu)) { pznu = tmpsol2;}
                                else pznu = tmpsol1;
                        }
                        if (type == 2 ) {
                                // pick the most central root.
                                if (TMath::Abs(tmpsol1)<TMath::Abs(tmpsol2) ) pznu = tmpsol1;
                                else pznu = tmpsol2;
                        }
                        if (type == 3 ) {
                                // pick the largest value of the cosine
                                TVector3 p3w, p3mu;
                                p3w.SetXYZ(pxmu+pxnu, pymu+pynu, pzmu+ tmpsol1);
                                p3mu.SetXYZ(pxmu, pymu, pzmu );
                                
                                double sinthcm1 = 2.*(p3mu.Perp(p3w))/M_W;
                                p3w.SetXYZ(pxmu+pxnu, pymu+pynu, pzmu+ tmpsol2);
                                double sinthcm2 = 2.*(p3mu.Perp(p3w))/M_W;

                                double costhcm1 = TMath::Sqrt(1. - sinthcm1*sinthcm1);
                                double costhcm2 = TMath::Sqrt(1. - sinthcm2*sinthcm2);

                                if ( costhcm1 > costhcm2 ) pznu = tmpsol1;
                                else pznu = tmpsol2;
                        }
                
        }

        //Particle neutrino;
        //neutrino.setP4( LorentzVector(pxnu, pynu, pznu, TMath::Sqrt(pxnu*pxnu + pynu*pynu + pznu*pznu ))) ;

        return pznu;
}