- // ROOT
- #include "TMath.h"
-
- // Gaudi
- #include "GaudiKernel/PhysicalConstants.h"
- #include "GaudiUtils/Aida2ROOT.h"
-
- // Tb
- #include "TbKernel/ITbGeometrySvc.h"
- #include "Event/TbCluster.h"
- #include "Event/TbTrack.h"
- #include "Event/TbState.h"
-
- #include "Event/TbKalmanTrack.h"
- #include "Event/TbKalmanNode.h"
- #include "Event/TbKalmanPixelMeasurement.h"
-
- #include <math.h>
-
- // Local
- #include "TbTrackFitter.h"
-
- /** @file TbTrackFitter.cpp
- *
- * Implementation of class : TbTrackFitter
- * Author: Panagiotis Tsopelas
- */
-
- DECLARE_ALGORITHM_FACTORY(TbTrackFitter)
-
- //=============================================================================
- /// Standard constructor
- //=============================================================================
- TbTrackFitter::TbTrackFitter(const std::string& name, ISvcLocator* pSvcLocator)
- : GaudiTupleAlg(name, pSvcLocator), m_trackFit(NULL) {
-
- // Global default values.
- declareProperty("onOff_hists", m_onOff_hists = true);
- declareProperty("scat2", m_scat2 = 1.35e-8);
- declareProperty("hiterror2", m_hiterror2 = 1.6e-5);
- declareProperty("theta0", m_theta0 = 1.e-4);
- declareProperty("direction", m_direction = -1);
- // Local default values.
- m_nPlanes = 8;
- sigmax = 0.004 * Gaudi::Units::mm;
- sigmay = 0.004 * Gaudi::Units::mm;
- }
-
- //=============================================================================
- /// Initialization
- //=============================================================================
- StatusCode TbTrackFitter::initialize() {
-
- StatusCode sc = GaudiAlgorithm::initialize();
- if (sc.isFailure()) return sc;
- setHistoTopDir("Tb/");
- if (m_onOff_hists) setup_hists();
-
- // Setup the track fit tool.
- m_trackFit = tool<ITbTrackFit>("TbTrackFit", "Fitter", this);
-
- // Initialize the random number generator.
- sc = m_gauss.initialize(randSvc(), Rndm::Gauss(0.0, 1.0));
- if (!sc) {
- error() << "Cannot initialize Gaussian random number generator." << endmsg;
- return sc;
- }
-
- return StatusCode::SUCCESS;
- }
-
- //=============================================================================
- /// Main execution
- //=============================================================================
- StatusCode TbTrackFitter::execute() {
- if (msgLevel(MSG::DEBUG)) debug() << " ==> execute()" << endmsg;
-
- // TbTrack container
- std::vector<LHCb::TbTrack*> tracks_vec;
- // TbKalmanTrack container
- std::vector<LHCb::TbKalmanTrack*> ktracks_vec;
-
- const unsigned int nPlanes = m_nPlanes;
- // const double sigmax = 0.004 * Gaudi::Units::mm;
- // const double sigmay = 0.004 * Gaudi::Units::mm;
-
- // Gap from next/previous event
- // std::cout << std::endl << std::endl << std::endl ;
-
- // make random track parameters -> True straight track - NO scattering
- double Rx0 = m_gauss() * 15. * Gaudi::Units::mm ;
- double Rtx = m_gauss() * 1.e-4 ;
- double Ry0 = m_gauss() * 15. * Gaudi::Units::mm ;
- double Rty = m_gauss() * 1.e-4 ;
- // print random initial track parameters
- // std::cout << "--> Random track with x0 = " << Rx0 << " and tx = " << Rtx << std::endl;
-
- // theta0 for MCS through small angles
- double theta0 = TMath::ATan(m_theta0) ;
- // store true slope and intercept from scattering
- std::vector<double> scatRtx(nPlanes);
- std::vector<double> scatRty(nPlanes);
-
- // Declare TbTrack
- LHCb::TbTrack* track = new LHCb::TbTrack();
-
- // make random clusters ..
- std::vector<LHCb::TbCluster> clusters(nPlanes);
-
- // .. but fix their Z positions according to the Alignment File
- clusters[0].setZ( 0. ); clusters[4].setZ( 306. );
- clusters[1].setZ( 31. ); clusters[5].setZ( 336. );
- clusters[2].setZ( 60. ); clusters[6].setZ( 366. );
- clusters[3].setZ( 89. ); clusters[7].setZ( 397. );
-
- const double weight = 1. / (0.004 * 0.004);
- // initialize 1st cluster
- clusters[0].setX( ( Rtx + Rx0 ) + m_gauss()*sigmax ); // random straight line
- clusters[0].setY( ( Rty + Ry0 ) + m_gauss()*sigmay ); // random straight line
- clusters[0].setXErr( 0.004 );
- clusters[0].setYErr( 0.004 );
- clusters[0].setWx(weight);
- clusters[0].setWy(weight);
- clusters[0].setPlane(0);
- track->addToClusters(&clusters[0]);
-
- // std::cout << " Clusters in X: " ;
- // std::cout << clusters[0].x() << ", " ;
- // update x0 and tx due to scattering
- scatRx0[0] = Rx0 ;
- scatRtx[0] = Rtx + m_gauss()*theta0;
- scatRy0[0] = Ry0;
- scatRty[0] = Rty + m_gauss()*theta0;
-
- // fill the rest of the clusters
- for (unsigned int i = 1; i < nPlanes; ++i) {
-
- // update x0 and tx due to scattering
- scatRx0[i] = ( clusters[i].z() - clusters[i-1].z() ) * scatRtx[i-1] + scatRx0[i-1] ;
- scatRtx[i] = scatRtx[i-1] + m_gauss() * theta0;
- scatRy0[i] = ( clusters[i].z() - clusters[i-1].z() ) * scatRty[i-1] + scatRy0[i-1];
- scatRty[i] = scatRty[i-1] + m_gauss() * theta0;
-
- clusters[i].setX( scatRx0[i] + m_gauss()*sigmax );
- clusters[i].setY( scatRy0[i] + m_gauss()*sigmay );
- clusters[i].setXErr( 0.004 );
- clusters[i].setYErr( 0.004 );
- clusters[i].setWx(weight);
- clusters[i].setWy(weight);
- clusters[i].setPlane(i);
- track->addToClusters(&clusters[i]);
- if (i==nPlanes) continue;
-
- }
-
- // const SmartRefVector<LHCb::TbCluster>& clustersL = track->clusters();
- // std::cout << std::endl << "Track with " << track->size() << " clusters. " << std::endl;
- // // print true positions of track
- // std::cout << " True X positions of (scattered) track: ";
- // std::cout << Rx0 << ", " ;
- // for (int i=0; i<nPlanes; i++){
- // std::cout << scatRx0[i] << ", " ;
- // }
- // std::cout << std::endl;
- // std::cout << " Cluster X positions : ";
- // for (auto it = clustersL.cbegin(), end = clustersL.cend(); it != end; ++it) {
- // std::cout << (*it)->x() << ", ";
- // }
- // std::cout << " Cluster X errors : ";
- // for (auto it = clustersL.cbegin(), end = clustersL.cend(); it != end; ++it) {
- // std::cout << (*it)->xErr() << ", ";
- // }
- //// for (int i=0; i<nPlanes; i++){
- //// std::cout << clusters[i].x() << ", " ;
- //// }
- //
- // std::cout << std::endl;
- // std::cout << " Cluster Z positions : ";
- // for (auto it = clustersL.cbegin(), end = clustersL.cend(); it != end; ++it) {
- // std::cout << (*it)->z() << ", ";
- // }
- //// for (int i=0; i<nPlanes; i++){
- //// std::cout << clusters[i].z() << ", " ;
- //// }
- // std::cout << std::endl;
-
-
- // fit the track
- m_trackFit->fit(track);
-
-
-
- // print results of simple fit
- // std::cout << std::endl << " * chi2 on x-z plane: " << track->chi2() << " with ndof = " << track->ndof() << std::endl
- // << " x0: " << track->firstState().x() << " +/- " << TMath::Sqrt(track->firstState().errX2() )
- // << ", y0: " << track->firstState().y() << " +/- " << TMath::Sqrt(track->firstState().errY2() )
- // << std::endl
- // << " tx: " << track->firstState().tx() << " +/- " << TMath::Sqrt( track->firstState().errTx2() )
- // << ", ty: " << track->firstState().ty() << " +/- " << TMath::Sqrt( track->firstState().errTy2() )
- // << std::endl ;
-
-
-
-
-
- //---------------------------------------------------------------------------
- // Wouter's new code
- //---------------------------------------------------------------------------
-
- // node : combination of a measurement and a track state is referred to as a node
- // std::cout << std::endl << " ================================================== " << std::endl;
-
- // create a fit track object (which is actually also a TbTrack)
- LHCb::TbKalmanTrack* ktrack = new LHCb::TbKalmanTrack(*track, m_hiterror2, m_scat2) ;
- // fit it
- ktrack->fit() ;
-
- // dump some info about the fit
- // ktrack->print() ;
- // std::cout << " ================================================== " << std::endl;
- //---------------------------------------------------------------------------
-
-
- // store the track in the track vector
- tracks_vec.push_back(track);
- // store the ktrack in the track vector
- ktracks_vec.push_back(ktrack);
-
- // Check whether to fill histograms
- if (m_onOff_hists) {
- fill_hists(tracks_vec);
- fill_khists(ktracks_vec);
- }
-
- // delete the pointer to the track
- delete track;
-
- return StatusCode::SUCCESS;
- }
-
- //=============================================================================
- /// F U N C T I O N S
- //=============================================================================
-
- //=============================================================================
- /// Fill Histograms for TbTracks
- //=============================================================================
- void TbTrackFitter::fill_hists(std::vector<LHCb::TbTrack*>& tracks) {
-
- std::vector<LHCb::TbTrack*>::iterator ictra;
- for (ictra = tracks.begin(); ictra != tracks.end(); ictra++) {
-
- // if ( TMath::Prob((*ictra)->chi2() , (*ictra)->ndof()) < 0.05) continue;
-
- // Fill the track histos
- m_slopesX->Fill((*ictra)->firstState().tx());
- m_slopesY->Fill((*ictra)->firstState().ty());
- m_Sfit_chi2->Fill((*ictra)->chi2PerNdof());
- m_Sfit_prob->Fill( TMath::Prob((*ictra)->chi2() , (*ictra)->ndof()) );
-
- // Get the clusters of this TbTrack
- const SmartRefVector<LHCb::TbCluster> clusters = (*ictra)->clusters();
-
- // Loop through the clusters of this TbTrack
- SmartRefVector<LHCb::TbCluster>::const_iterator icclu;
- for (icclu = clusters.begin(); icclu != clusters.end(); ++icclu) {
- int ichip = (*icclu)->plane();
-
- m_clustersX->Fill((*icclu)->x());
- m_clustersY->Fill((*icclu)->y());
-
- double xTraAtZclu = (*ictra)->firstState().x() +
- (*ictra)->firstState().tx() * (double)(*icclu)->z();
- double yTraAtZclu = (*ictra)->firstState().y() +
- (*ictra)->firstState().ty() * (double)(*icclu)->z();
-
- double resx = (*icclu)->x() - xTraAtZclu ;
- double resy = (*icclu)->y() - yTraAtZclu ;
-
- // Fill residuals
- m_resSfit_X[ichip]->Fill( resx );
- m_resSfit_Y[ichip]->Fill( resy );
-
- // Fill residual pulls
- m_respullSfit_X[ichip]->Fill((double) resx / sigmax);
- m_respullSfit_Y[ichip]->Fill((double) resy / sigmay);
-
- // Fill differences from true track
- m_trueSfit_X[ichip]->Fill( scatRx0[ichip] - xTraAtZclu );
- m_trueSfit_Y[ichip]->Fill( scatRy0[ichip] - yTraAtZclu );
-
- // Fill true pulls
- m_trackpullSfit_X[ichip]->Fill( (scatRx0[ichip] - xTraAtZclu) / std::sqrt((*ictra)->firstState().covariance()(0,0)) );
- m_trackpullSfit_Y[ichip]->Fill( (scatRy0[ichip] - yTraAtZclu) / std::sqrt((*ictra)->firstState().covariance()(1,1)) );
- // m_trackpullSfit_X[ichip]->Fill( (scatRx0[ichip] - xTraAtZclu) / sigmax );
- // m_trackpullSfit_Y[ichip]->Fill( (scatRy0[ichip] - yTraAtZclu) / sigmay );
-
-
-
- } // end of cluster loop
-
- } // end of track loop
- }
-
- //=============================================================================
- /// Fill Histograms for TbKalmanTracks
- //=============================================================================
- void TbTrackFitter::fill_khists(std::vector<LHCb::TbKalmanTrack*>& ktracks) {
-
- std::vector<LHCb::TbKalmanTrack*>::iterator icktra;
- for (icktra = ktracks.begin(); icktra != ktracks.end(); icktra++) {
-
- // Fill the track histos
- m_Kfit_chi2->Fill((*icktra)->chi2());
- m_Kfit_prob->Fill( TMath::Prob((*icktra)->chi2(), (*icktra)->ndof()) );
-
- // node->residualX() / std::sqrt(pixelhit->residualCovX()
-
- // Get the nodes of this TbKalmanTrack
- //const std::vector<LHCb::TbKalmanNode*>& knodes = (*icktra)->nodes();
-
- // Loop through the nodes of this TbKalmanTrack
- for( auto node : (*icktra)->nodes() ) {
- auto pixnode = dynamic_cast< LHCb::TbKalmanPixelMeasurement*>( node) ;
- if( pixnode ) {
- int ichip = pixnode->cluster().plane();
-
- // Fill unbiased residuals
- m_XunresKfit[ichip]->Fill( pixnode->residualX() * pixnode->covX() / pixnode->residualCovX() );
- m_YunresKfit[ichip]->Fill( pixnode->residualY() * pixnode->covY() / pixnode->residualCovY() );
-
- // Fill residuals
- m_resKfit_X[ichip]->Fill( pixnode->residualX() );
- m_resKfit_Y[ichip]->Fill( pixnode->residualY() );
-
- // Fill residual errors
- m_reserrKfit_X[ichip]->Fill( std::sqrt(pixnode->residualCovX()) );
- m_reserrKfit_Y[ichip]->Fill( std::sqrt(pixnode->residualCovY()) );
-
- // Fill residual pulls
- m_respullKfit_X[ichip]->Fill( pixnode->residualX() / std::sqrt(pixnode->residualCovX() ) );
- m_respullKfit_Y[ichip]->Fill( pixnode->residualY() / std::sqrt(pixnode->residualCovY() ) );
-
- // Fill differences from true tracks
- m_trueKfit_X[ichip]->Fill( scatRx0[ichip] - pixnode->state().x() );
- m_trueKfit_Y[ichip]->Fill( scatRy0[ichip] - pixnode->state().y() );
-
- // Fill differences errors from true tracks
- m_trueerrKfit_X[ichip]->Fill( TMath::Sqrt(pixnode->state().covariance()(0,0)) );
- m_trueerrKfit_Y[ichip]->Fill( TMath::Sqrt(pixnode->state().covariance()(1,1)) );
-
- // Fill true pulls
- m_trackpullKfit_X[ichip]->Fill( (scatRx0[ichip] - pixnode->state().x() ) / TMath::Sqrt(pixnode->state().covariance()(0,0)) );
- m_trackpullKfit_Y[ichip]->Fill( (scatRy0[ichip] - pixnode->state().y() ) / TMath::Sqrt(pixnode->state().covariance()(1,1)) );
-
- }
- }
- } // end of Ktrack loop
- }
-
-
- //=============================================================================
- /// Setup Histograms
- //=============================================================================
- void TbTrackFitter::setup_hists() {
-
- // TbTrack parameters plots
- m_Sfit_chi2 = Gaudi::Utils::Aida2ROOT::aida2root(
- book1D("StraightLineFit/chi2perndof", "Chi2/ndof", -0.5, 49.5, 100));
-
- m_Sfit_prob = Gaudi::Utils::Aida2ROOT::aida2root(
- book1D("StraightLineFit/probability", "Chi2 prob of S fit", 0.0, 1.0, 100));
-
- m_slopesX = Gaudi::Utils::Aida2ROOT::aida2root(
- book1D("StraightLineFit/slopesX", "slopes in X of TbTracks", -1.e-3, 1.e-3, 100));
-
- m_slopesY = Gaudi::Utils::Aida2ROOT::aida2root(
- book1D("StraightLineFit/slopesY", "slopes in Y of TbTracks", -1.e-3, 1.e-3, 100));
-
- m_Kfit_chi2 = Gaudi::Utils::Aida2ROOT::aida2root(
- book1D("KalmanFit/chi2", "Chi2", -0.5, 49.5, 100));
-
- m_Kfit_prob = Gaudi::Utils::Aida2ROOT::aida2root(
- book1D("KalmanFit/probability", "Chi2 prob of K fit", 0.0, 1.0, 100));
-
- m_clustersX = Gaudi::Utils::Aida2ROOT::aida2root(
- book1D("clustersX", "cluster X positions ", -1.e-2, 1.e-2, 100));
-
- m_clustersY = Gaudi::Utils::Aida2ROOT::aida2root(
- book1D("clustersY", "cluster Y positions ", -1.e-2, 1.e-2, 100));
-
-
-
- // Residual plots
- std::string hist_name;
- for (unsigned int i = 0; i < m_nPlanes; i++) {
- std::stringstream ss_chip;
- ss_chip << i;
-
- // straight line fit
- hist_name = "StraightLineFit/Residuals_on_X/plane_" + ss_chip.str();
- m_resSfit_X.push_back(Gaudi::Utils::Aida2ROOT::aida2root(
- book1D(hist_name.c_str(), hist_name.c_str(), -1.e-2, 1.e-2, 100)));
-
- hist_name = "StraightLineFit/Residuals_on_Y/plane_" + ss_chip.str();
- m_resSfit_Y.push_back(Gaudi::Utils::Aida2ROOT::aida2root(
- book1D(hist_name.c_str(), hist_name.c_str(), -1.e-2, 1.e-2, 100)));
-
-
- hist_name = "StraightLineFit/ResidualPull_on_X/plane_" + ss_chip.str();
- m_respullSfit_X.push_back(Gaudi::Utils::Aida2ROOT::aida2root(
- book1D(hist_name.c_str(), hist_name.c_str(), -10., 10., 100)));
-
- hist_name = "StraightLineFit/ResidualPull_on_Y/plane_" + ss_chip.str();
- m_respullSfit_Y.push_back(Gaudi::Utils::Aida2ROOT::aida2root(
- book1D(hist_name.c_str(), hist_name.c_str(), -10., 10., 100)));
-
-
- hist_name = "StraightLineFit/Differences_from_true_on_X/plane_" + ss_chip.str();
- m_trueSfit_X.push_back(Gaudi::Utils::Aida2ROOT::aida2root(
- book1D(hist_name.c_str(), hist_name.c_str(), -1.e-2, 1.e-2, 100)));
-
- hist_name = "StraightLineFit/Differences_from_true_on_Y/plane_" + ss_chip.str();
- m_trueSfit_Y.push_back(Gaudi::Utils::Aida2ROOT::aida2root(
- book1D(hist_name.c_str(), hist_name.c_str(), -1.e-2, 1.e-2, 100)));
-
-
- hist_name = "StraightLineFit/TrackPull_on_X/plane_" + ss_chip.str();
- m_trackpullSfit_X.push_back(Gaudi::Utils::Aida2ROOT::aida2root(
- book1D(hist_name.c_str(), hist_name.c_str(), -10., 10., 100)));
-
- hist_name = "StraightLineFit/TrackPull_on_Y/plane_" + ss_chip.str();
- m_trackpullSfit_Y.push_back(Gaudi::Utils::Aida2ROOT::aida2root(
- book1D(hist_name.c_str(), hist_name.c_str(), -10., 10., 100)));
-
-
-
- // Kalman fit
- hist_name = "KalmanFit/UnbiasedResidualsX/plane_" + ss_chip.str();
- m_XunresKfit.push_back(Gaudi::Utils::Aida2ROOT::aida2root(
- book1D(hist_name.c_str(), hist_name.c_str(), -1.e-2, 1.e-2, 100)));
-
-
- hist_name = "KalmanFit/UnbiasedResidualsY/plane_" + ss_chip.str();
- m_YunresKfit.push_back(Gaudi::Utils::Aida2ROOT::aida2root(
- book1D(hist_name.c_str(), hist_name.c_str(), -1.e-2, 1.e-2, 100)));
-
- //
-
- hist_name = "KalmanFit/Residuals_on_X/plane_" + ss_chip.str();
- m_resKfit_X.push_back(Gaudi::Utils::Aida2ROOT::aida2root(
- book1D(hist_name.c_str(), hist_name.c_str(), -1.e-2, 1.e-2, 100)));
-
-
- hist_name = "KalmanFit/Residuals_on_Y/plane_" + ss_chip.str();
- m_resKfit_Y.push_back(Gaudi::Utils::Aida2ROOT::aida2root(
- book1D(hist_name.c_str(), hist_name.c_str(), -1.e-2, 1.e-2, 100)));
-
- hist_name = "KalmanFit/Residual_errors_on_X/plane_" + ss_chip.str();
- m_reserrKfit_X.push_back(Gaudi::Utils::Aida2ROOT::aida2root(
- book1D(hist_name.c_str(), hist_name.c_str(), 0., 5.e-3, 1000)));
-
-
- hist_name = "KalmanFit/Residual_errors_on_Y/plane_" + ss_chip.str();
- m_reserrKfit_Y.push_back(Gaudi::Utils::Aida2ROOT::aida2root(
- book1D(hist_name.c_str(), hist_name.c_str(), 0., 5.e-3, 1000)));
-
-
- hist_name = "KalmanFit/ResidualPull_on_X/plane_" + ss_chip.str();
- m_respullKfit_X.push_back(Gaudi::Utils::Aida2ROOT::aida2root(
- book1D(hist_name.c_str(), hist_name.c_str(), -10., 10., 100)));
-
-
- hist_name = "KalmanFit/ResidualPull_on_Y/plane_" + ss_chip.str();
- m_respullKfit_Y.push_back(Gaudi::Utils::Aida2ROOT::aida2root(
- book1D(hist_name.c_str(), hist_name.c_str(), -10., 10., 100)));
-
-
- hist_name = "KalmanFit/Differences_from_true_on_X/plane_" + ss_chip.str();
- m_trueKfit_X.push_back(Gaudi::Utils::Aida2ROOT::aida2root(
- book1D(hist_name.c_str(), hist_name.c_str(), -1.e-2, 1.e-2, 100)));
-
-
- hist_name = "KalmanFit/Differences_from_true_on_Y/plane_" + ss_chip.str();
- m_trueKfit_Y.push_back(Gaudi::Utils::Aida2ROOT::aida2root(
- book1D(hist_name.c_str(), hist_name.c_str(), -1.e-2, 1.e-2, 100)));
-
- hist_name = "KalmanFit/Differences_errors_from_true_on_X/plane_" + ss_chip.str();
- m_trueerrKfit_X.push_back(Gaudi::Utils::Aida2ROOT::aida2root(
- book1D(hist_name.c_str(), hist_name.c_str(), 0., 5.e-3, 1000)));
-
-
- hist_name = "KalmanFit/Differences_errors_from_true_on_Y/plane_" + ss_chip.str();
- m_trueerrKfit_Y.push_back(Gaudi::Utils::Aida2ROOT::aida2root(
- book1D(hist_name.c_str(), hist_name.c_str(), 0., 5.e-3, 1000)));
-
-
- hist_name = "KalmanFit/TrackPull_on_X/plane_" + ss_chip.str();
- m_trackpullKfit_X.push_back(Gaudi::Utils::Aida2ROOT::aida2root(
- book1D(hist_name.c_str(), hist_name.c_str(), -10., 10., 100)));
-
-
- hist_name = "KalmanFit/TrackPull_on_Y/plane_" + ss_chip.str();
- m_trackpullKfit_Y.push_back(Gaudi::Utils::Aida2ROOT::aida2root(
- book1D(hist_name.c_str(), hist_name.c_str(), -10., 10., 100)));
- }
- }
-
- //=============================================================================
- /// End
- //=============================================================================