// Example from the Oxbridge Stransverse Mass Library -- oxbridgekinetics. // See http://www.hep.phy.cam.ac.uk/~lester/mt2/index.html // Authors: Christopher Lester and Alan Barr #include "ExampleEvent.h" #include "Mt2/Advanced_Mt2_332_Calculator.h" #include inline double sq(const double x) { return x*x; } int main(int argc, char * argv[]) { // For the example. we now need some momenta and masses from which // to calculate MT2. // // In "reality" we would get these momenta from an ntuple or // from a physics event. // // As this is just an example program, we will instead get // some "example" momenta from the class "ExampleEvent" // defined in "ExampleEvent.h"as follows: ExampleEvent exampleEvent; while(true){ exampleEvent.resetRandom332(); Mt2::LorentzTransverseVector ltv_Vis_A = exampleEvent. ltv_Vis_A(); Mt2::LorentzTransverseVector ltv_Vis_B = exampleEvent. ltv_Vis_B(); const double Q = 4.*(exampleEvent.flipshoot()); Mt2::TwoVector visSum = (ltv_Vis_A.transverse() + ltv_Vis_B.transverse()); Mt2::TwoVector pT_Miss = visSum * Q; double chi = exampleEvent.invis_mass(); std::cout << "Going to calculate MT2 with\n" << " ltv_Vis_A = " << ltv_Vis_A << "\n" << " ltv_Vis_B = " << ltv_Vis_B << "\n" << " pT_Miss = " << pT_Miss << "\n" << " invis_mass = " << chi << std::endl; // Now that we have some visiable momenta and some missing transverse // momentum we can calculate MT2. // First we create the object that is going to do the calculation // of MT2 for us. // // For this example we will use a modification of the "332" aglorithm // that was originally defined in SUSYPhys. Our modification (called // Basic_Mt2_332_Calculator) is basically the same as // SUSYPhys_Mt2_222_Calculator except that we remove the assumption // that visible particles are massless. Mt2::Advanced_Mt2_332_Calculator mt2Calculator(true); // Could tell the MT2 calculating object to be verbose, and print out // debug messages while it is thinking ... but we won't: // mt2Calculator.setDebug(true); // Now we can actually calculate MT2: const double mt2 = mt2Calculator.mt2_332( ltv_Vis_A, ltv_Vis_B, pT_Miss, chi); Mt2::SolutionType type = mt2Calculator.lastSolutionType(); const double ma = ltv_Vis_A.mass(); const double mb = ltv_Vis_B.mass(); const double AT = ltv_Vis_A.contralinearDot(ltv_Vis_B); /* The following is correct for Q=0 or Q=1 or ma=mb, and has "tuned" powers (the second (1.+Q) term and teh pow(fabs(1+Q),2) term) to make it do "well" elsewhere. */ const double other = chi*chi + (1.+Q)*(ma*ma + mb*mb)/2. + (1.+Q)*sq(mb*mb-ma*ma) / (2.*(2.*AT - ma*ma - mb*mb)) - AT*Q + sqrt( (AT*AT - ma*ma *mb*mb) * ( sq((mb*mb-ma*ma)/(2.*AT-ma*ma-mb*mb))*pow(fabs(1.+Q),2) + (4.* chi*chi)/(2.* AT - ma*ma - mb*mb) + Q*Q ) ); // Now we print out the result: std::cout << "ANSWER: mt2 = " << mt2 << " for " << mt2Calculator.algorithmName() << " algorithm" << std::endl; std::cout << "COMPARE " << type << " " << mt2 << " " << sqrt(other) << " " << Q << " " << 2.*AT - ma*ma - mb*mb <<" " << AT*AT - ma*ma *mb*mb << std::endl; } return 0; }