TopEvent.C

Go to the documentation of this file.

// Class TopEvent
//
// Author: Benno List, Jenny Boehme
// Last update: $Date: 2005/01/12 10:11:46 $
//          by: $Author: blist $
// 
// Description: class to generate and fit top pair events
//               
#include "jbltools/kinfit/TopEvent.h"

#include "jbltools/kinfit/JetFitObject.h"
#include "jbltools/kinfit/MassConstraint.h"
#include "jbltools/kinfit/cernlib.h"

#include <iostream>              // - cout
using std::cout;
using std::endl;

// constructor: 
TopEvent::TopEvent()
: pxc (1, 0),
  pyc (0, 1)
 {
  for (int i = 0; i < NFV; ++i) fv[i] = 0;
  for (int i = 0; i < NBFO; ++i) bfo[i] = bfosmear[i] = 0;
  w1.setMass(80.4);
  w2.setMass(80.4);
};

//destructor: 
TopEvent::~TopEvent() {
  for (int i = 0; i < NFV; ++i) delete fv[i];
  for (int i = 0; i < NBFO; ++i) {
    delete bfo[i];
    delete bfosmear[i];
  }  
};

// generate x value
double TopEvent::genX(double lambda, double xmin) {
  // Generate an x value between xmin and 1 distributed like x^-lambda,
  // lambda mustnot be -1!
  assert (lambda != -1);
  using std::pow;
  FReal randoms[1];
  ranmar (randoms, 1);
  
  double ymin = pow (xmin, 1-lambda);
  return pow (randoms[0]*(1-ymin)+ymin, 1./(1-lambda));
};



// generate four vectors
void TopEvent::genEvent(){
  // generate 4-vectors of top-decay:
  // 0: top-top-system -> top1 top2
  // 1: top 1 -> W1 b1
  // 2: top 2 -> W2 b2
  // 3: W1 -> j11 j12
  // 4: W2 -> j21 j22
  // 5: b1
  // 6: j11
  // 7: j12
  // 8: b2
  // 9: j21
  //10: j22
  
  double mtop = 174;
  double mW   = 80.4;
  double mb   = 5.0;
  double mj   = 1.0;
  
  double lambda = 1.5;
  double Ebeam = 900;
  double s = 4*Ebeam*Ebeam;
  double shatmin = 4*mtop*mtop;
  double xmin = shatmin/s;
  
  // x1, x2 are distributed according to x^-lambda
  double x1, x2, shat;
  do {
    x1 = genX (lambda, xmin);
    x2 = genX (lambda, xmin);
    shat = x1*x2*s;
//    std::cout << "x1=" << x1 << ", x2=" << x2 << ", shat=" << shat << endl;
  } while (shat < shatmin);
  
  FourVector *toppair = fv[0] = new FourVector ((x1+x2)*Ebeam, 0, 0, (x1-x2)*Ebeam);
  FourVector *top1 = fv[1] = new FourVector (mtop, 0, 0, 0);
  FourVector *top2 = fv[2] = new FourVector (mtop, 0, 0, 0);
  
  toppair->decayto (*top1, *top2);
  
  FourVector *W1 = fv[3] = new FourVector (mW, 0, 0, 0);
  FourVector *W2 = fv[4] = new FourVector (mW, 0, 0, 0);
  FourVector *b1 = fv[5] = new FourVector (mb, 0, 0, 0);
  FourVector *b2 = fv[8] = new FourVector (mb, 0, 0, 0);
  
  top1->decayto (*W1, *b1);
  top2->decayto (*W2, *b2);
  
  FourVector *j11 = fv[6]  = new FourVector (mj, 0, 0, 0);
  FourVector *j12 = fv[7]  = new FourVector (mj, 0, 0, 0);
  FourVector *j21 = fv[9]  = new FourVector (mj, 0, 0, 0);
  FourVector *j22 = fv[10] = new FourVector (mj, 0, 0, 0);
  
  W1->decayto (*j11, *j12);
  W2->decayto (*j21, *j22);
  
  double Eresol = 0.35;     // 35% / sqrt (E)
  double thetaResol = 0.1;  // rad
  double phiResol = 0.1;    // rad
  
  for (int j = 0; j < 6; ++j) {
    int i = j+5;
    double E = fv[i]->getE();
    double theta = fv[i]->getTheta();
    double phi = fv[i]->getPhi();
    double EError = Eresol*sqrt(E);
    
    // Create fit object with true quantities for later comparisons
    bfo[j] = new JetFitObject (E, theta, phi, EError, thetaResol, phiResol, 0.);
    
    FReal randoms[3];
    rnorml (randoms, 3);
    
    // Create fit object with smeared quantities as fit input
    double ESmear = E + EError*randoms[0];
    double thetaSmear = theta + thetaResol*randoms[1];
    double phiSmear = phi + phiResol*randoms[2];
    
    bfosmear[j] = new JetFitObject (ESmear, thetaSmear, phiSmear, EError, thetaResol, phiResol, 0.);
    pxc.addToFOList (*bfosmear[j]);
    pyc.addToFOList (*bfosmear[j]);
    w.addToFOList (*bfosmear[j], j<3?1:2);
      
  }
  
  w1.addToFOList (*bfosmear[1]);
  w1.addToFOList (*bfosmear[2]);
  w2.addToFOList (*bfosmear[4]);
  w2.addToFOList (*bfosmear[5]);
  
};

// fit it!
int TopEvent::fitEvent (BaseFitter& fitter){
  
//   for (int i = 0; i < 6; ++i) 
//     cout << "true four-vector of jet " << i << ": " << *bfo[i] << endl;
//   for (int i = 0; i < 6; ++i) 
//     cout << "initial four-vector of jet " << i << ": " << *bfosmear[i] << endl;
  
  
  // reset lists of constraints and fitobjects
  fitter.reset();
   
  for (int i = 0; i < 6; i++) {
    assert (bfosmear[i]);
    fitter.addFitObject (*bfosmear[i]);
//     pxc.addToFOList (*bfosmear[i]);
//     pyc.addToFOList (*bfosmear[i]);
//     w.addToFOList (*bfosmear[i], i<3?1:2);
  }
//   
//   w1.addToFOList (*bfosmear[1]);
//   w1.addToFOList (*bfosmear[2]);
//   w2.addToFOList (*bfosmear[4]);
//   w2.addToFOList (*bfosmear[5]);
  
  fitter.addConstraint (pxc);
  fitter.addConstraint (pyc);
  fitter.addConstraint (w);
  fitter.addConstraint (w1);
  fitter.addConstraint (w2);
  
  double prob = fitter.fit();
  
//   cout << "fit probability = " << prob << endl;
//   for (int i = 0; i < 6; ++i) 
//     cout << "final four-vector of jet " << i << ": " << *bfosmear[i] << endl;
//   cout << "Constraint pxc: " << pxc.getValue() << endl;
//   cout << "Constraint pxc: " << pxc.getValue() << endl;
//   cout << "Constraint w:   " << w.getValue() << ", top mass: " << w.getMass() << endl;
//   cout << "Constraint w1:  " << w1.getValue() << ", W mass: " << w1.getMass() << endl;
//   cout << "Constraint w2:  " << w2.getValue() << ", W mass: " << w2.getMass() << endl;
   cout << "fit probability = " << prob << ", top mass: " << w.getMass() << endl;

   
   return fitter.getError();


};

---