VertexFit.h

Go to the documentation of this file.

#ifndef _VERTEXFIT_H
#define _VERTEXFIT_H

/*
VertexFit.h

Definition of the Vertex fit of dimuon events, this package is aimed at:
1. Use closest distance method to find z0 for two muon tracks
2. Use vertex fit to fit the z0 --- will be implemented in next version

Reference: CBM-SOFT-note-2006-001, by S. Gorbunov and I. Kisel, with some minor
modifications

Author: Kun Liu, liuk@fnal.gov
Created: 2-8-2012
*/

// Fun4All includes
#include <fun4all/SubsysReco.h>
#include <GlobalConsts.h>

#include <iostream>
#include <vector>

#include <TFile.h>
#include <TTree.h>
#include <TMatrixD.h>

#include "KalmanUtil.h"
#include "SRecEvent.h"
#include "FastTracklet.h"
#include "GenFitExtrapolator.h"

class KalmanTrack;
class KalmanFilter;
class PHField;
class TGeoManager;
class JobOptsSvc;

class VtxPar
{
public:
    VtxPar()
    {
        _r.ResizeTo(3, 1);
        _cov.ResizeTo(3, 3);

        _r.Zero();
        _cov.Zero();
    }

    void print()
    {
        SMatrix::printMatrix(_r, "Vertex position:");
        SMatrix::printMatrix(_cov, "Vertex covariance:");
    }

    TMatrixD _r;
    TMatrixD _cov;
};

class VertexFit : public SubsysReco
{
public:
    VertexFit(const std::string& name = "VertexFit");
    ~VertexFit();

    int Init(PHCompositeNode *topNode);
    int InitRun(PHCompositeNode *topNode);
    int process_event(PHCompositeNode *topNode);
    int End(PHCompositeNode *topNode);

    void enableOptimization() { optimize = true; }

    void setControlParameter(int nMaxIteration, double tolerance)
    {
        _max_iteration = nMaxIteration;
        _tolerance = tolerance;
    }

    int setRecEvent(SRecEvent* recEvent, int sign1 = 1, int sign2 = -1);

    void init();
    void addHypothesis(double z, double sigz = 50.) { z_start.push_back(z); sig_z_start.push_back(sigz); }
    void setStartingVertex(double z_start, double sigz_start);

    void addTrack(int index, SRecTrack& _track);
    void addTrack(int index, KalmanTrack& _track);
    void addTrack(int index, TrkPar& _trkpar);

    int processOnePair();

    int findVertex();
    double findDimuonVertexFast(SRecTrack& track1, SRecTrack& track2);
    double findSingleMuonVertex(SRecTrack& _track);
    double findSingleMuonVertex(Node& _node_start);
    double findSingleMuonVertex(TrkPar& _trkpar_start);

    double getVertexZ0() { return _vtxpar_curr._r[2][0]; }
    double getVXChisq() { return _chisq_vertex; }
    double getKFChisq() { return _chisq_kalman; }
    int getNTracks() { return _trkpar_curr.size(); }

    void updateVertex();

    void bookEvaluation(std::string evalFileName = "vtx_eval.root");
    void fillEvaluation();

    void print();

  const std::string& get_eval_file_name() const {
    return evalFileName;
  }

  void set_eval_file_name(const std::string& evalFileName) {
    this->evalFileName = evalFileName;
  }

private:

    int InitField(PHCompositeNode *topNode);

    int InitGeom(PHCompositeNode *topNode);

    int MakeNodes(PHCompositeNode *topNode);

    int GetNodes(PHCompositeNode *topNode);

    std::vector<TrkPar> _trkpar_curr;

    VtxPar _vtxpar_curr;

    Node _node_vertex;

    KalmanFilter* _kmfit;

    double _chisq_vertex;
    double _chisq_kalman;

    std::vector<double> z_start;
    std::vector<double> sig_z_start;

    std::vector<double> z_vertex;
    std::vector<double> r_vertex;
    std::vector<double> chisq_km;
    std::vector<double> chisq_vx;

    int _max_iteration;
    double _tolerance;

    GenFitExtrapolator _extrapolator;

    bool optimize;

    TGeoManager * _t_geo_manager;

    SRecEvent* _recEvent;

    std::string evalFileName;
    TFile* evalFile;
    TTree* evalTree;

    int runID;
    int eventID;
    int targetPos;
    int nPos;
    int nNeg;
    int p_idx_eval;
    int m_idx_eval;
    int choice_eval;
    int choice_by_kf_eval;
    int choice_by_vx_eval;

    int nStart;
    double z_start_eval[50];
    int nIter_eval[50];
    double chisq_kf_eval[50];
    double chisq_vx_eval[50];
    double z_vertex_eval[50];
    double r_vertex_eval[50];

    double m_chisq_kf_eval;
    double s_chisq_kf_eval;
    double m_z_vertex_eval;
    double s_z_vertex_eval;
};

#endif