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; }
 //Fitting in the target center (Abi)
    void enable_fit_target_center() {fit_target_center = 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;
    bool fit_target_center;
 
    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