k0fit.C

Go to the documentation of this file.

// k0fit
// Test program to test track fitter classes:
// Generate a K0S, let it decay, and fit it
// B. List 24.11.04


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

#include "jbltools/kinfit/WWFitterFast.h"
#include "jbltools/kinfit/WWFitter.h"
#include "jbltools/kinfit/K0Event.h"
#include "jbltools/kinfit/ftypes.h"
#include "jbltools/kinfit/hbook.h"
#include "jbltools/kinfit/TrackFitObject.h"



int main(int argc, char** argv) {
  int failed = 0;
  int total = 2;
  
  // ntuple variables
  FInteger IERR;
  FReal FITPROB, CHI2;
  FReal V0TRUE[3], V1TRUE[3];
  FReal V0FIT[3], V1FIT[3], V2FIT[3];
  FReal P0TRUE[4], P1TRUE[4], P2TRUE[4];
  FReal P1SM[4], P2SM[4];
  FReal P0FIT[4], P1FIT[4], P2FIT[4];
  FReal PAR0TRUE[6], PAR1TRUE[6], PAR2TRUE[6];
  FReal PAR1SM[6], PAR2SM[6];
  FReal PAR1SMERR[6], PAR2SMERR[6];
  FReal PAR0FIT[6], PAR1FIT[6], PAR2FIT[6];
  FReal PAR0ERR[6], PAR1ERR[6], PAR2ERR[6];
  FReal DL2TRUE, DL2FIT, DL3TRUE, DL3FIT;
  FReal TTRUE, TFIT;
  FReal MTRUE, MFIT;
  
  
  FInteger LREC=1024, ISTAT=0;
  hropen (21, "K0FIT", "k0fit.hbook", "N", LREC, ISTAT);
  hbnt (999,"k0fit","");
  hbname (999, "FIT", IERR, "IERR:I");
  hbname (999, "FIT", FITPROB, "FITPROB:R");
  hbname (999, "FIT", CHI2, "CHI2:R");
  hbname (999, "FIT", V0TRUE[0], "V0TRUE(3):R");
  hbname (999, "FIT", V1TRUE[0], "V1TRUE(3):R");
  hbname (999, "FIT", V0FIT[0], "V0FIT(3):R");
  hbname (999, "FIT", V1FIT[0], "V1FIT(3):R");
  hbname (999, "FIT", V2FIT[0], "V2FIT(3):R");
  hbname (999, "FIT", P0TRUE[0], "P0TRUE(4):R");
  hbname (999, "FIT", P1TRUE[0], "P1TRUE(4):R");
  hbname (999, "FIT", P2TRUE[0], "P2TRUE(4):R");
  hbname (999, "FIT", P1SM[0], "P1SM(4):R");
  hbname (999, "FIT", P2SM[0], "P2SM(4):R");
  hbname (999, "FIT", P0FIT[0], "P0FIT(4):R");
  hbname (999, "FIT", P1FIT[0], "P1FIT(4):R");
  hbname (999, "FIT", P2FIT[0], "P2FIT(4):R");
  hbname (999, "FIT", PAR0TRUE[0], "PAR0TRUE(6):R");
  hbname (999, "FIT", PAR1TRUE[0], "PAR1TRUE(6):R");
  hbname (999, "FIT", PAR2TRUE[0], "PAR2TRUE(6):R");
  hbname (999, "FIT", PAR1SM[0], "PAR1SM(6):R");
  hbname (999, "FIT", PAR2SM[0], "PAR2SM(6):R");
  hbname (999, "FIT", PAR1SMERR[0], "PAR1SMERR(6):R");
  hbname (999, "FIT", PAR2SMERR[0], "PAR2SMERR(6):R");
  hbname (999, "FIT", PAR0FIT[0], "PAR0FIT(6):R");
  hbname (999, "FIT", PAR1FIT[0], "PAR1FIT(6):R");
  hbname (999, "FIT", PAR2FIT[0], "PAR2FIT(6):R");
  hbname (999, "FIT", PAR0ERR[0], "PAR0ERR(6):R");
  hbname (999, "FIT", PAR1ERR[0], "PAR1ERR(6):R");
  hbname (999, "FIT", PAR2ERR[0], "PAR2ERR(6):R");
  
  hbname (999, "FIT", DL2TRUE, "DL2TRUE:R");
  hbname (999, "FIT", DL2FIT, "DL2FIT:R");
  hbname (999, "FIT", DL3TRUE, "DL3TRUE:R");
  hbname (999, "FIT", DL3FIT, "DL3FIT:R");
  hbname (999, "FIT", TTRUE, "TTRUE:R");
  hbname (999, "FIT", TFIT, "TFIT:R");
  hbname (999, "FIT", MTRUE, "MTRUE:R");
  hbname (999, "FIT", MFIT, "MFIT:R");
  
  
  
  for (int i = 0; i < total; i++) {
    WWFitter fitter;
    K0Event evt;
    evt.genEvent();    
    
    if ((IERR = evt.fitEvent(fitter)) != 0) {
      failed++;
      cout << "IERR = " << IERR << endl;
    }
    FITPROB = fitter.getProbability();
    CHI2    = fitter.getChi2();
    
    // vertices
    for (int i = 0; i < 3; ++i) {
      V0TRUE[i] = evt.genvert[0]->getComponent(i);
      V1TRUE[i] = evt.genvert[1]->getComponent(i);
      V0FIT[i] = evt.rectrack[0]->getTrajectoryPoint(0).getComponent(i);
      V1FIT[i] = evt.rectrack[1]->getVertex(0).getComponent(i);
      V2FIT[i] = evt.rectrack[2]->getVertex(0).getComponent(i);
    }
    
    // 4-momenta
    for (int i = 0; i < 4; ++i) {
      P0TRUE[i] = evt.gentrack[0]->getMomentum(0).getComponent(i);
      P1TRUE[i] = evt.gentrack[1]->getMomentum(0).getComponent(i);
      P2TRUE[i] = evt.gentrack[2]->getMomentum(0).getComponent(i);
      P1SM[i] = evt.smtrack[1]->getMomentum(0).getComponent(i);
      P2SM[i] = evt.smtrack[2]->getMomentum(0).getComponent(i);
      P0FIT[i] = evt.rectrack[0]->getMomentum(0).getComponent(i);
      P1FIT[i] = evt.rectrack[1]->getMomentum(0).getComponent(i);
      P2FIT[i] = evt.rectrack[2]->getMomentum(0).getComponent(i);
    }
    // parameters & errors
    for (int i = 0; i < 6; ++i) {
      PAR0TRUE[i] = evt.gentrack[0]->getParam(i);
      PAR1TRUE[i] = evt.gentrack[1]->getParam(i);
      PAR2TRUE[i] = evt.gentrack[2]->getParam(i);
      PAR1SM[i] = evt.smtrack[1]->getParam(i);
      PAR2SM[i] = evt.smtrack[2]->getParam(i);
      PAR1SMERR[i] = evt.smtrack[1]->getError(i);
      PAR2SMERR[i] = evt.smtrack[2]->getError(i);
      PAR0FIT[i] = evt.rectrack[0]->getParam(i);
      PAR1FIT[i] = evt.rectrack[1]->getParam(i);
      PAR2FIT[i] = evt.rectrack[2]->getParam(i);
      PAR0ERR[i] = evt.rectrack[0]->getError(i);
      PAR1ERR[i] = evt.rectrack[1]->getError(i);
      PAR2ERR[i] = evt.rectrack[2]->getError(i);
    }
    
    // decay length, 2d/3d
    
    DL2TRUE = (*evt.genvert[1]-*evt.genvert[0]).getPt();
    DL2FIT = (evt.rectrack[0]->getVertex(0)-evt.rectrack[0]->getTrajectoryPoint(0)).getPt();
    DL3TRUE = (*evt.genvert[1]-*evt.genvert[0]).getP();
    DL3FIT = (evt.rectrack[0]->getVertex(0)-evt.rectrack[0]->getTrajectoryPoint(0)).getP();
  
    TTRUE = DL3TRUE/(2.998E-2*(evt.gentrack[0]->getMomentum(0).getBetaGamma().getMag()));
    TFIT = DL3FIT/(2.998E-2*(evt.rectrack[0]->getMomentum(0).getBetaGamma().getMag()));
  
    MTRUE = (evt.gentrack[1]->getMomentum(0) + evt.gentrack[2]->getMomentum(0)).getM();
    MFIT = (evt.rectrack[1]->getMomentum(0) + evt.rectrack[2]->getMomentum(0)).getM();
    
    hfnt(999);
  }
  
  std::cerr << "Failed fits: " << failed 
            << "(" << (100.0*failed)/total << "%)" << endl;
            
  hrout(999);
  hrend("K0FIT");          

  return 0;
}

---