Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-03 08:15:04 

0001 #include "SVReco.h"
0002 /*#include <trackbase_historic/SvtxCluster.h>
0003 #include <trackbase_historic/SvtxClusterMap.h>
0004 #include <trackbase_historic/SvtxHitMap.h>*/
0005 #include <trackbase_historic/SvtxTrack.h>
0006 #include <trackbase_historic/SvtxTrack_v1.h>
0007 #include <trackbase_historic/SvtxTrackState_v1.h>
0008 #include <trackbase_historic/SvtxVertex.h>
0009 #include <trackbase_historic/SvtxVertex_v1.h>
0010 #include <trackbase_historic/SvtxTrackMap.h>
0011 #include <trackbase_historic/SvtxVertexMap.h>
0012 
0013 #include <trackbase/TrkrClusterContainer.h>
0014 #include <trackbase/TrkrClusterv1.h>
0015 #include <trackbase/TrkrHitSet.h>
0016 #include <trackbase/TrkrHitSetContainer.h>
0017 #include <trackbase/TrkrClusterHitAssoc.h>
0018 #include <mvtx/MvtxDefs.h>
0019 #include <intt/InttDefs.h>
0020 
0021 #include <g4jets/JetMap.h>
0022 
0023 #include <fun4all/Fun4AllReturnCodes.h>
0024 #include <fun4all/PHTFileServer.h>
0025 
0026 //#include <g4detectors/PHG4CellContainer.h>
0027 #include <g4detectors/PHG4Cell.h>
0028 #include <g4detectors/PHG4CylinderGeomContainer.h>
0029 #include <mvtx/CylinderGeom_Mvtx.h>
0030 #include <intt/CylinderGeomIntt.h>
0031 
0032 #include <g4main/PHG4Hit.h>
0033 #include <g4main/PHG4HitContainer.h>
0034 #include <g4main/PHG4TruthInfoContainer.h>
0035 #include <g4main/PHG4Particle.h>
0036 
0037 #include <phgenfit/Fitter.h>
0038 #include <phgenfit/PlanarMeasurement.h>
0039 #include <phgenfit/Track.h>
0040 #include <phgenfit/SpacepointMeasurement.h>
0041 
0042 #include <phool/getClass.h>
0043 #include <phool/phool.h>
0044 #include <phool/PHCompositeNode.h>
0045 #include <phool/PHIODataNode.h>
0046 #include <phool/PHNodeIterator.h>
0047 
0048 #include <phgeom/PHGeomUtility.h>
0049 #include <phfield/PHFieldUtility.h>
0050 
0051 #include <GenFit/FieldManager.h>
0052 #include <GenFit/GFRaveVertex.h>
0053 #include <GenFit/GFRaveVertexFactory.h>
0054 #include <GenFit/MeasuredStateOnPlane.h>
0055 #include <GenFit/RKTrackRep.h>
0056 #include <GenFit/StateOnPlane.h>
0057 #include <GenFit/Track.h>
0058 #include <GenFit/KalmanFitterInfo.h>
0059 
0060 //#include <HFJetTruthGeneration/HFJetDefs.h>
0061 
0062 #include <TClonesArray.h>
0063 #include <TMatrixDSym.h>
0064 #include <TTree.h>
0065 #include <TVector3.h>
0066 #include <TRandom3.h>
0067 #include <TRotation.h>
0068 
0069 #include <iostream>
0070 #include <utility>
0071 #include <memory>
0072 
0073 #define LogDebug(exp)       std::cout<<"DEBUG: "  <<__FILE__<<": "<<__LINE__<<": "<< exp <<std::endl
0074 #define LogError(exp)       std::cout<<"ERROR: "  <<__FILE__<<": "<<__LINE__<<": "<< exp <<std::endl
0075 #define LogWarning(exp) std::cout<<"WARNING: "<<__FILE__<<": "<<__LINE__<<": "<< exp <<std::endl
0076 
0077 using namespace std;
0078 
0079 /*Rave
0080 #include <rave/Version.h>
0081 //#include <rave/Track.h>
0082 #include <rave/VertexFactory.h>
0083 #include <rave/ConstantMagneticField.h>
0084 */
0085 //GenFit
0086 //#include <GenFit/GFRaveConverters.h>
0087 
0088 
0089 /*
0090  * Constructor
0091  */
0092 SVReco::SVReco(const string &name) :
0093     _mag_field_file_name("/phenix/upgrades/decadal/fieldmaps/sPHENIX.2d.root"),
0094     _mag_field_re_scaling_factor(1.4 / 1.5), //what is this and why?
0095     _reverse_mag_field(false),
0096     _fitter(NULL),
0097     _track_fitting_alg_name("DafRef"),
0098     _n_maps_layer(3),
0099     _n_intt_layer(4),
0100     _primary_pid_guess(11), //for the tracks
0101     _cut_Ncluster(false),
0102     _cut_min_pT(0.1),
0103     _cut_dca(5.0), //probably in mm
0104     _cut_chi2_ndf(5),
0105     _use_ladder_geom(false),
0106     _vertex_finder(NULL),
0107     _vertexing_method("avf-smoothing:1"), /*need a list of these and their proper domains*/
0108     _clustermap(NULL),
0109     _trackmap(NULL),
0110     _vertexmap(NULL),
0111     _do_eval(false),
0112     _verbosity(10),
0113     _do_evt_display(false)
0114 {
0115 
0116 }
0117 
0118 
0119 int SVReco::InitEvent(PHCompositeNode *topNode) {
0120   cout<<"getting SVR nodes"<<endl;
0121     GetNodes(topNode);
0122     //cout<<"got vertexing nodes"<<endl;
0123     //! stands for Refit_GenFit_Tracks
0124     for(auto p : _main_rf_phgf_tracks) delete p;
0125     _main_rf_phgf_tracks.clear();
0126 
0127     _svtxtrk_gftrk_map.clear();
0128     _svtxtrk_wt_map.clear();
0129     _svtxtrk_id.clear();
0130 
0131     //iterate over all tracks to find priary vertex and make rave/genfit objects
0132   cout<<"start SVR track loop"<<endl;
0133     for (SvtxTrackMap::Iter iter = _trackmap->begin(); iter != _trackmap->end();
0134             ++iter) {
0135         SvtxTrack* svtx_track = iter->second;
0136         // do track cuts
0137         if (!svtx_track || svtx_track->get_ndf()<40 || svtx_track->get_pt()<_cut_min_pT ||
0138                 (svtx_track->get_chisq()/svtx_track->get_ndf())>_cut_chi2_ndf ||
0139                 fabs(svtx_track->get_dca3d_xy())>_cut_dca || fabs(svtx_track->get_dca3d_z())>_cut_dca)
0140             continue;
0141 
0142         int n_MVTX = 0, n_INTT = 0, n_TPC = 0;
0143         //cout<<"Keys:";
0144         for (SvtxTrack::ConstClusterKeyIter iter2 = svtx_track->begin_cluster_keys(); iter2!=svtx_track->end_cluster_keys(); iter2++) {
0145             TrkrDefs::cluskey cluster_key = *iter2;
0146             //cout<<cluster_key<<',';
0147             //count where the hits are
0148             float layer = (float) TrkrDefs::getLayer(cluster_key);
0149             if (layer<_n_maps_layer) n_MVTX++;
0150             else if (layer<_n_maps_layer+_n_intt_layer) n_INTT++;
0151             else n_TPC++;
0152         }//cluster loop
0153         //cluster cuts
0154         //cout<<"\n cluster loop with n_MVTX="<<n_MVTX<<" n_INTT="<<n_INTT<<" and nTPC="<<n_TPC<<endl;
0155         if ( _cut_Ncluster && (n_MVTX<2 || n_INTT<2 || n_TPC<25) ){
0156             continue;
0157         }
0158         //cout << (svtx_track->get_chisq()/svtx_track->get_ndf()) << ", " << n_TPC << ", " << svtx_track->get_pt() << endl;
0159         //cout << svtx_track->get_ndf() << ", " << svtx_track->size_clusters() << endl;
0160         //cout<<"making genfit"<<endl;
0161         PHGenFit::Track* rf_phgf_track = MakeGenFitTrack(svtx_track); //convert SvtxTrack to GenFit Track
0162         //cout<<"made genfit"<<endl;
0163 
0164         if (rf_phgf_track) {
0165             //svtx_track->identify();
0166             //make a map to connect SvtxTracks to their respective GenFit Tracks
0167             _svtxtrk_id.push_back(svtx_track->get_id());
0168             _svtxtrk_gftrk_map[svtx_track->get_id()] = _main_rf_phgf_tracks.size();
0169             _main_rf_phgf_tracks.push_back(rf_phgf_track); //to be used by findSecondaryVerticies
0170         }
0171     }
0172     cout<<"exit track loop ntracks="<<_main_rf_phgf_tracks.size()<<'\n';
0173 
0174     return Fun4AllReturnCodes::EVENT_OK;
0175 }
0176 
0177 PHGenFit::Track* SVReco::getPHGFTrack(SvtxTrack* svtxtrk){
0178  if(svtxtrk)return _main_rf_phgf_tracks[_svtxtrk_gftrk_map[svtxtrk->get_id()]];
0179  else return NULL;
0180 }
0181 
0182 
0183 int SVReco::InitRun(PHCompositeNode *topNode) {
0184 
0185     CreateNodes(topNode);
0186 
0187     TGeoManager* tgeo_manager = PHGeomUtility::GetTGeoManager(topNode);
0188     PHField * field = PHFieldUtility::GetFieldMapNode(nullptr, topNode);
0189 
0190     _fitter = PHGenFit::Fitter::getInstance(tgeo_manager, field,
0191             _track_fitting_alg_name, "RKTrackRep", _do_evt_display);
0192 
0193     if (!_fitter) {
0194         cerr << PHWHERE << endl;
0195         return Fun4AllReturnCodes::ABORTRUN;
0196     }
0197 
0198     _vertex_finder = new genfit::GFRaveVertexFactory(_verbosity);
0199     _vertex_finder->setMethod(_vertexing_method.data());
0200 
0201     if (!_vertex_finder) {
0202         std::cerr<< PHWHERE<<" bad run init no SVR"<<endl;
0203         return Fun4AllReturnCodes::ABORTRUN;
0204     }
0205 
0206     return Fun4AllReturnCodes::EVENT_OK;
0207 }
0208 
0209 //note that this might only work for conversion like events
0210 genfit::GFRaveVertex* SVReco::findSecondaryVertex(SvtxTrack* track1, SvtxTrack* track2) {
0211     //_vertex_finder->setMethod("avr-smoothing:1");
0212     _vertex_finder->setMethod("avr");
0213     //_vertex_finder->setMethod("avf-smoothing:1");
0214     //_vertex_finder->setMethod("kalman");
0215     vector<genfit::GFRaveVertex*> rave_vertices_conversion;
0216     vector<genfit::Track*> rf_gf_tracks_conversion;
0217 
0218     if (_svtxtrk_gftrk_map.find(track1->get_id())!=_svtxtrk_gftrk_map.end()&&
0219             _svtxtrk_gftrk_map.find(track2->get_id())!=_svtxtrk_gftrk_map.end())
0220     {
0221         unsigned int trk_index = _svtxtrk_gftrk_map[track1->get_id()];
0222         PHGenFit::Track* rf_phgf_track = _main_rf_phgf_tracks[trk_index];
0223         rf_gf_tracks_conversion.push_back(rf_phgf_track->getGenFitTrack());
0224         //check the genfit track is working
0225         /*cout<<"Track Comparison Original:\n";
0226             track1->identify();*/
0227         //printGenFitTrackKinematics(rf_phgf_track);
0228 
0229         trk_index = _svtxtrk_gftrk_map[track2->get_id()];
0230         rf_phgf_track = _main_rf_phgf_tracks[trk_index];
0231         //    track2->identify();
0232         // printGenFitTrackKinematics(rf_phgf_track);
0233         rf_gf_tracks_conversion.push_back(rf_phgf_track->getGenFitTrack());
0234     }
0235     if (rf_gf_tracks_conversion.size()>1){
0236         try{
0237             _vertex_finder->findVertices(&rave_vertices_conversion, rf_gf_tracks_conversion);
0238         }catch (...){
0239             std::cout << PHWHERE << "GFRaveVertexFactory::findVertices failed!";
0240         }
0241         /*      if(TMath::Abs(rave_vertices_conversion[0]->getChi2()-1.)>.3){
0242                         cout<<"PRINTING:\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\n";
0243                         rave_vertices_conversion[0]->Print();
0244                         track1->identify();
0245                         track2->identify();
0246                         cout<<"\n\n\n\n\n";
0247                         }*/
0248         //if a vertex is found return ownership 
0249     //TODO check rave_verticies_conversion for mem leak
0250         if(rave_vertices_conversion.size()>0) return rave_vertices_conversion[0];
0251         else return NULL;
0252     }//more than 1 track 
0253     else return NULL;
0254 }
0255 
0256 SVReco::~SVReco(){
0257     cout<<PHWHERE<<"delete"<<endl;
0258     if(_fitter){
0259         delete _fitter;
0260         _fitter=NULL;
0261     }
0262     cout<<PHWHERE<<"delete"<<endl;
0263     if(_vertex_finder){
0264         delete _vertex_finder;
0265         _vertex_finder=NULL;
0266     }
0267     cout<<PHWHERE<<"delete"<<endl;
0268     for (std::vector<PHGenFit::Track*>::iterator i = _main_rf_phgf_tracks.begin(); i != _main_rf_phgf_tracks.end(); ++i)
0269     {
0270         if(*i)
0271             delete *i;
0272         *i=NULL;
0273     }
0274 }
0275 
0276 //prints the track using the trival verex for point extrapolation
0277 void SVReco::printGenFitTrackKinematics(PHGenFit::Track* track){
0278 
0279     TVector3 vertex_position(0, 0, 0);
0280 
0281     std::shared_ptr<genfit::MeasuredStateOnPlane> gf_state_vertex_ca = NULL;
0282     try {
0283         gf_state_vertex_ca = std::shared_ptr < genfit::MeasuredStateOnPlane> (track->extrapolateToPoint(vertex_position));
0284     } catch (...) {
0285         if (_verbosity >= 2)
0286             LogWarning("extrapolateToPoint failed!");
0287     }
0288     TVector3 mom = gf_state_vertex_ca->getMom();
0289     TMatrixDSym cov = gf_state_vertex_ca->get6DCov();
0290 
0291     cout << "OUT Ex: " << sqrt(cov[0][0]) << ", Ey: " << sqrt(cov[1][1]) << ", Ez: " << sqrt(cov[2][2]) << endl;
0292     cout << "OUT Px: " << mom.X() << ", Py: " << mom.Y() << ", Pz: " << mom.Z() << endl; 
0293 }
0294 
0295 
0296 
0297 //should be deprecated
0298 void SVReco::reset_eval_variables(){
0299     gf_prim_vtx_ntrk = rv_prim_vtx_ntrk = 0;
0300     for (int i=0; i<3; i++){
0301         gf_prim_vtx[i] = gf_prim_vtx_err[i] = -999;
0302         rv_prim_vtx[i] = rv_prim_vtx_err[i] = -999;
0303     }//i
0304 
0305     rv_sv_njets = 0;
0306 
0307     for (int ijet=0; ijet<10; ijet++){
0308         rv_sv_jet_id[ijet] = -999;
0309         rv_sv_jet_prop[ijet][0] = rv_sv_jet_prop[ijet][1] = -999;
0310         rv_sv_jet_pT[ijet] = -999;
0311         rv_sv_jet_px[ijet] = rv_sv_jet_py[ijet] = rv_sv_jet_pz[ijet] = -999;
0312 
0313         rv_sv_pT00_nvtx[ijet] = rv_sv_pT05_nvtx[ijet] = rv_sv_pT10_nvtx[ijet] = rv_sv_pT15_nvtx[ijet] = 0;
0314 
0315         for (int ivtx=0; ivtx<30; ivtx++){
0316             rv_sv_pT00_vtx_ntrk[ijet][ivtx] = rv_sv_pT05_vtx_ntrk[ijet][ivtx] = rv_sv_pT10_vtx_ntrk[ijet][ivtx] = rv_sv_pT15_vtx_ntrk[ijet][ivtx] = 0;
0317             rv_sv_pT00_vtx_ntrk_good[ijet][ivtx] = rv_sv_pT05_vtx_ntrk_good[ijet][ivtx] = rv_sv_pT10_vtx_ntrk_good[ijet][ivtx] = rv_sv_pT15_vtx_ntrk_good[ijet][ivtx] = 0;
0318             rv_sv_pT00_vtx_ntrk_good_pv[ijet][ivtx] = rv_sv_pT05_vtx_ntrk_good_pv[ijet][ivtx] = rv_sv_pT10_vtx_ntrk_good_pv[ijet][ivtx] = rv_sv_pT15_vtx_ntrk_good_pv[ijet][ivtx] = 0;
0319 
0320             rv_sv_pT00_vtx_mass[ijet][ivtx] = rv_sv_pT00_vtx_mass_corr[ijet][ivtx] = rv_sv_pT00_vtx_pT[ijet][ivtx] = -999;
0321             rv_sv_pT05_vtx_mass[ijet][ivtx] = rv_sv_pT05_vtx_mass_corr[ijet][ivtx] = rv_sv_pT05_vtx_pT[ijet][ivtx] = -999;
0322             rv_sv_pT10_vtx_mass[ijet][ivtx] = rv_sv_pT10_vtx_mass_corr[ijet][ivtx] = rv_sv_pT10_vtx_pT[ijet][ivtx] = -999;
0323             rv_sv_pT15_vtx_mass[ijet][ivtx] = rv_sv_pT15_vtx_mass_corr[ijet][ivtx] = rv_sv_pT15_vtx_pT[ijet][ivtx] = -999;
0324 
0325             rv_sv_pT00_vtx_x[ijet][ivtx] = rv_sv_pT00_vtx_y[ijet][ivtx] = rv_sv_pT00_vtx_z[ijet][ivtx] = -999;
0326             rv_sv_pT00_vtx_ex[ijet][ivtx] = rv_sv_pT00_vtx_ey[ijet][ivtx] = rv_sv_pT00_vtx_ez[ijet][ivtx] = -999;
0327             rv_sv_pT05_vtx_x[ijet][ivtx] = rv_sv_pT05_vtx_y[ijet][ivtx] = rv_sv_pT05_vtx_z[ijet][ivtx] = -999;
0328             rv_sv_pT05_vtx_ex[ijet][ivtx] = rv_sv_pT05_vtx_ey[ijet][ivtx] = rv_sv_pT05_vtx_ez[ijet][ivtx] = -999;
0329             rv_sv_pT10_vtx_x[ijet][ivtx] = rv_sv_pT10_vtx_y[ijet][ivtx] = rv_sv_pT10_vtx_z[ijet][ivtx] = -999;
0330             rv_sv_pT10_vtx_ex[ijet][ivtx] = rv_sv_pT10_vtx_ey[ijet][ivtx] = rv_sv_pT10_vtx_ez[ijet][ivtx] = -999;
0331             rv_sv_pT15_vtx_x[ijet][ivtx] = rv_sv_pT15_vtx_y[ijet][ivtx] = rv_sv_pT15_vtx_z[ijet][ivtx] = -999;
0332             rv_sv_pT15_vtx_ex[ijet][ivtx] = rv_sv_pT15_vtx_ey[ijet][ivtx] = rv_sv_pT15_vtx_ez[ijet][ivtx] = -999;
0333 
0334             rv_sv_pT00_vtx_jet_theta[ijet][ivtx] = rv_sv_pT00_vtx_pT[ijet][ivtx] = -999;
0335             rv_sv_pT05_vtx_jet_theta[ijet][ivtx] = rv_sv_pT05_vtx_pT[ijet][ivtx] = -999;
0336             rv_sv_pT10_vtx_jet_theta[ijet][ivtx] = rv_sv_pT10_vtx_pT[ijet][ivtx] = -999;
0337             rv_sv_pT15_vtx_jet_theta[ijet][ivtx] = rv_sv_pT15_vtx_pT[ijet][ivtx] = -999;
0338 
0339             rv_sv_pT00_vtx_chi2[ijet][ivtx] = rv_sv_pT00_vtx_ndf[ijet][ivtx] = -999;
0340             rv_sv_pT05_vtx_chi2[ijet][ivtx] = rv_sv_pT05_vtx_ndf[ijet][ivtx] = -999;
0341             rv_sv_pT10_vtx_chi2[ijet][ivtx] = rv_sv_pT10_vtx_ndf[ijet][ivtx] = -999;
0342             rv_sv_pT15_vtx_chi2[ijet][ivtx] = rv_sv_pT15_vtx_ndf[ijet][ivtx] = -999;
0343 
0344         }//ivtx
0345     }//ijet
0346 
0347     return;
0348 }
0349 
0350 int SVReco::CreateNodes(PHCompositeNode *topNode){
0351 
0352     return Fun4AllReturnCodes::EVENT_OK;
0353 }
0354 
0355 /*
0356  * GetNodes():
0357  *  Get all the all the required nodes off the node tree
0358  */
0359 int SVReco::GetNodes(PHCompositeNode * topNode){
0360     _clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0361     if (!_clustermap){
0362         cout << PHWHERE << " TRKR_CLUSTERS node not found on node tree"
0363             << endl;
0364         return Fun4AllReturnCodes::ABORTEVENT;
0365     }
0366     // Input Svtx Tracks
0367     _trackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0368     if (!_trackmap){
0369         cout << PHWHERE << " SvtxTrackMap node not found on node tree"
0370             << endl;
0371         return Fun4AllReturnCodes::ABORTEVENT;
0372     }
0373     // Input Svtx Vertices
0374     _vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
0375     if (!_vertexmap){
0376         cout << PHWHERE << " SvtxVertexrMap node not found on node tree"
0377             << endl;
0378         return Fun4AllReturnCodes::ABORTEVENT;
0379     }
0380 
0381     _geom_container_intt = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_INTT");
0382     if (!_geom_container_intt)
0383     {
0384         cout << PHWHERE << "CYLINDERGEOM_INTT node not found on node tree"
0385             << endl;
0386         return Fun4AllReturnCodes::ABORTEVENT;
0387     }
0388     _geom_container_maps = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MVTX");
0389         if (!_geom_container_maps)
0390         {
0391             cout << PHWHERE << "CYLINDERGEOM_MVTX node not found on node tree"
0392                 << endl;
0393             return Fun4AllReturnCodes::ABORTEVENT;
0394         }
0395     if(!_vertex_finder){
0396         std::cerr<< PHWHERE<<" bad run init no SVR"<<endl;
0397         return Fun4AllReturnCodes::ABORTRUN;
0398     } 
0399 
0400 
0401     return Fun4AllReturnCodes::EVENT_OK;
0402 }
0403 
0404 //Edited version original from @sh-lim
0405 PHGenFit::Track* SVReco::MakeGenFitTrack(const SvtxTrack* intrack){
0406     if (!intrack){
0407         cerr << PHWHERE << " Input SvtxTrack is NULL!" << endl;
0408         return NULL;
0409     }
0410 
0411     if (_use_ladder_geom and !_geom_container_intt and !_geom_container_maps) {
0412         cout << PHWHERE << "No PHG4CylinderGeomContainer found!" << endl;
0413         return NULL;
0414     }
0415 
0416     TVector3 seed_pos(intrack->get_x(), intrack->get_y(), intrack->get_z());
0417     TVector3 seed_mom(intrack->get_px(), intrack->get_py(), intrack->get_pz()); //mom stands for momentum
0418     TMatrixDSym seed_cov(6);
0419     for (int i=0; i<6; i++){
0420         for (int j=0; j<6; j++){
0421             seed_cov[i][j] = intrack->get_error(i,j);
0422         }
0423     }
0424 
0425     // Create measurements
0426     std::vector<PHGenFit::Measurement*> measurements;
0427 
0428     for (auto iter = intrack->begin_cluster_keys(); iter != intrack->end_cluster_keys(); ++iter){
0429         //    unsigned int cluster_id = *iter;
0430         TrkrCluster* cluster = _clustermap->findCluster(*iter);
0431         if (!cluster) {
0432             LogError("No cluster Found!");
0433             continue;
0434         }
0435         float x = cluster->getPosition(0);
0436         float y = cluster->getPosition(1);
0437         float radius = sqrt(x*x+y*y);
0438         TVector3 pos(cluster->getPosition(0), cluster->getPosition(1), cluster->getPosition(2));
0439         seed_mom.SetPhi(pos.Phi());
0440         seed_mom.SetTheta(pos.Theta());
0441 
0442         TVector3 n(cluster->getPosition(0), cluster->getPosition(1), 0);
0443         //cout<<"Cluster with {"<<cluster->getPosition(0)<<','<<cluster->getPosition(0)<<"}\n";
0444 
0445         if (_use_ladder_geom){ //I don't understand this bool
0446             unsigned int trkrid = TrkrDefs::getTrkrId(*iter);
0447             unsigned int layer = TrkrDefs::getLayer(*iter);
0448             if (trkrid == TrkrDefs::mvtxId) {
0449                 int stave_index = MvtxDefs::getStaveId(*iter);
0450                 int chip_index = MvtxDefs::getChipId(*iter);
0451 
0452                 double ladder_location[3] = { 0.0, 0.0, 0.0 };
0453                 //not exactly sure where the cylinder geoms are currently objectified. check this 
0454                 CylinderGeom_Mvtx *geom = (CylinderGeom_Mvtx*) _geom_container_maps->GetLayerGeom(layer);
0455                 // returns the center of the sensor in world coordinates - used to get the ladder phi location
0456                 geom->find_sensor_center(stave_index, 0, 0, chip_index, ladder_location);//the mvtx module and half stave are 0
0457                 n.SetXYZ(ladder_location[0], ladder_location[1], 0);
0458                 n.RotateZ(geom->get_stave_phi_tilt());
0459             } 
0460             else if (trkrid == TrkrDefs::inttId) {
0461                 //this may bug but it looks ok for now
0462                 CylinderGeomIntt* geom = (CylinderGeomIntt*) _geom_container_intt->GetLayerGeom(layer);
0463                 double hit_location[3] = { 0.0, 0.0, 0.0 };
0464                 geom->find_segment_center(InttDefs::getLadderZId(*iter),InttDefs::getLadderPhiId(*iter), hit_location);
0465 
0466                 n.SetXYZ(hit_location[0], hit_location[1], 0);
0467                 n.RotateZ(geom->get_strip_phi_tilt());
0468             }
0469         }//if use_ladder_geom
0470         PHGenFit::Measurement* meas = new PHGenFit::PlanarMeasurement(pos, n,radius*cluster->getPhiError(), cluster->getZError());
0471         measurements.push_back(meas);
0472     }//cluster loop
0473     genfit::AbsTrackRep* rep = new genfit::RKTrackRep(_primary_pid_guess);
0474     PHGenFit::Track* track(new PHGenFit::Track(rep, seed_pos, seed_mom, seed_cov));
0475     track->addMeasurements(measurements);
0476 
0477     if (_fitter->processTrack(track, false) != 0) {
0478         if (_verbosity >= 1)
0479             LogWarning("Track fitting failed");
0480         return NULL;
0481     }
0482 
0483     track->getGenFitTrack()->setMcTrackId(intrack->get_id());
0484 
0485     return track;
0486 }
0487 
0488 //inspired by PHG4TrackKalmanFitter
0489 PHGenFit::Track* SVReco::MakeGenFitTrack(const SvtxTrack* intrack, const SvtxVertex* invertex){
0490   if (!intrack){
0491     cerr << PHWHERE << " Input SvtxTrack is NULL!" << endl;
0492     return NULL;
0493   }
0494 
0495   if (_use_ladder_geom and !_geom_container_intt and !_geom_container_maps) {
0496     cout << PHWHERE << "No PHG4CylinderGeomContainer found!" << endl;
0497     return NULL;
0498   }
0499 
0500   // Create measurements
0501   std::vector<PHGenFit::Measurement*> measurements;
0502 
0503   //create space point measurement from vtx
0504   if (invertex and invertex->size_tracks() > 1) {
0505     TVector3 pos(invertex->get_x(), invertex->get_y(), invertex->get_z());
0506     TMatrixDSym cov(3);
0507     cov.Zero();
0508     bool is_vertex_cov_sane = true;
0509     cout<<"Making covarience for vtx measurement"<<endl;
0510     for (unsigned int i = 0; i < 3; i++)
0511       for (unsigned int j = 0; j < 3; j++) {
0512         cov(i, j) = invertex->get_error(i, j);
0513       }
0514 
0515     cout<<"Made covarience"<<endl;
0516     if (is_vertex_cov_sane) {
0517       PHGenFit::Measurement* meas = new PHGenFit::SpacepointMeasurement(
0518           pos, cov);
0519       measurements.push_back(meas);
0520     }
0521 
0522     //convert SvtxTrack to matricies
0523     TVector3 seed_pos(intrack->get_x(), intrack->get_y(), intrack->get_z());
0524     TVector3 seed_mom(intrack->get_px(), intrack->get_py(), intrack->get_pz()); //mom stands for momentum
0525     TMatrixDSym seed_cov(6);
0526     for (int i=0; i<6; i++){
0527       for (int j=0; j<6; j++){
0528         seed_cov[i][j] = intrack->get_error(i,j);
0529       }
0530     } 
0531     cout<<"Making track cluster measurments"<<endl;
0532     //make measurements from the track clusters
0533     for (auto iter = intrack->begin_cluster_keys(); iter != intrack->end_cluster_keys(); ++iter){
0534       //    unsigned int cluster_id = *iter;
0535       TrkrCluster* cluster = _clustermap->findCluster(*iter);
0536       if (!cluster) {
0537         LogError("No cluster Found!");
0538         continue;
0539       }
0540       float x = cluster->getPosition(0);
0541       float y = cluster->getPosition(1);
0542       float radius = sqrt(x*x+y*y);
0543       TVector3 pos(cluster->getPosition(0), cluster->getPosition(1), cluster->getPosition(2));
0544       seed_mom.SetPhi(pos.Phi());
0545       seed_mom.SetTheta(pos.Theta());
0546 
0547       TVector3 n(cluster->getPosition(0), cluster->getPosition(1), 0);
0548       //cout<<"Cluster with {"<<cluster->getPosition(0)<<','<<cluster->getPosition(0)<<"}\n";
0549 
0550       if (_use_ladder_geom){ //I don't understand this bool
0551         unsigned int trkrid = TrkrDefs::getTrkrId(*iter);
0552         unsigned int layer = TrkrDefs::getLayer(*iter);
0553         if (trkrid == TrkrDefs::mvtxId) {
0554           int stave_index = MvtxDefs::getStaveId(*iter);
0555           int chip_index = MvtxDefs::getChipId(*iter);
0556 
0557           double ladder_location[3] = { 0.0, 0.0, 0.0 };
0558           //not exactly sure where the cylinder geoms are currently objectified. check this 
0559           CylinderGeom_Mvtx *geom = (CylinderGeom_Mvtx*) _geom_container_maps->GetLayerGeom(layer);
0560           // returns the center of the sensor in world coordinates - used to get the ladder phi location
0561           geom->find_sensor_center(stave_index, 0, 0, chip_index, ladder_location);//the mvtx module and half stave are 0
0562           n.SetXYZ(ladder_location[0], ladder_location[1], 0);
0563           n.RotateZ(geom->get_stave_phi_tilt());
0564         } 
0565         else if (trkrid == TrkrDefs::inttId) {
0566           //this may bug but it looks ok for now
0567           CylinderGeomIntt* geom = (CylinderGeomIntt*) _geom_container_intt->GetLayerGeom(layer);
0568           double hit_location[3] = { 0.0, 0.0, 0.0 };
0569           geom->find_segment_center(InttDefs::getLadderZId(*iter),InttDefs::getLadderPhiId(*iter), hit_location);
0570 
0571           n.SetXYZ(hit_location[0], hit_location[1], 0);
0572           n.RotateZ(geom->get_strip_phi_tilt());
0573         }
0574       }//if use_ladder_geom
0575       PHGenFit::Measurement* meas = new PHGenFit::PlanarMeasurement(pos, n,radius*cluster->getPhiError(), cluster->getZError());
0576       measurements.push_back(meas);
0577     }//cluster loop
0578     genfit::AbsTrackRep* rep = new genfit::RKTrackRep(_primary_pid_guess);
0579     PHGenFit::Track* track(new PHGenFit::Track(rep, seed_pos, seed_mom, seed_cov));
0580     track->addMeasurements(measurements);
0581 
0582     if (_fitter->processTrack(track, false) != 0) {
0583       if (_verbosity >= 1)
0584         LogWarning("Track fitting failed");
0585       return NULL;
0586     }
0587     track->getGenFitTrack()->setMcTrackId(intrack->get_id());
0588     return track;
0589   }//valid vtx 
0590   else{
0591     cerr<<PHWHERE<<": invalid vertex"<<endl;
0592     return NULL;
0593   }
0594 
0595 }
0596 
0597 //inspired by PHG4TrackKalmanFitter
0598 PHGenFit::Track* SVReco::MakeGenFitTrack(const SvtxTrack* intrack, const genfit::GFRaveVertex* invertex){
0599   if (!intrack){
0600     cerr << PHWHERE << " Input SvtxTrack is NULL!" << endl;
0601     return NULL;
0602   }
0603 
0604   if (_use_ladder_geom and !_geom_container_intt and !_geom_container_maps) {
0605     cout << PHWHERE << "No PHG4CylinderGeomContainer found!" << endl;
0606     return NULL;
0607   }
0608 
0609   // Create measurements
0610   std::vector<PHGenFit::Measurement*> measurements;
0611 
0612   //create space point measurement from vtx
0613   //TODO check that getNTracks is properly initialized 
0614   if (invertex and invertex->getNTracks() > 1) {
0615     TVector3 pos=invertex->getPos();
0616     TMatrixDSym cov = invertex->getCov();
0617     PHGenFit::Measurement* meas = new PHGenFit::SpacepointMeasurement(
0618           pos, cov);
0619     measurements.push_back(meas);
0620 
0621     //convert SvtxTrack to matricies
0622     TVector3 seed_pos(intrack->get_x(), intrack->get_y(), intrack->get_z());
0623     TVector3 seed_mom(intrack->get_px(), intrack->get_py(), intrack->get_pz()); //mom stands for momentum
0624     TMatrixDSym seed_cov(6);
0625     for (int i=0; i<6; i++){
0626       for (int j=0; j<6; j++){
0627         seed_cov[i][j] = intrack->get_error(i,j);
0628       }
0629     } 
0630     cout<<"Making track cluster measurments"<<endl;
0631     //make measurements from the track clusters
0632     for (auto iter = intrack->begin_cluster_keys(); iter != intrack->end_cluster_keys(); ++iter){
0633       //    unsigned int cluster_id = *iter;
0634       TrkrCluster* cluster = _clustermap->findCluster(*iter);
0635       if (!cluster) {
0636         LogError("No cluster Found!");
0637         continue;
0638       }
0639       float x = cluster->getPosition(0);
0640       float y = cluster->getPosition(1);
0641       float radius = sqrt(x*x+y*y);
0642       TVector3 pos(cluster->getPosition(0), cluster->getPosition(1), cluster->getPosition(2));
0643       seed_mom.SetPhi(pos.Phi());
0644       seed_mom.SetTheta(pos.Theta());
0645 
0646       TVector3 n(cluster->getPosition(0), cluster->getPosition(1), 0);
0647       //cout<<"Cluster with {"<<cluster->getPosition(0)<<','<<cluster->getPosition(0)<<"}\n";
0648 
0649       if (_use_ladder_geom){ //I don't understand this bool
0650         unsigned int trkrid = TrkrDefs::getTrkrId(*iter);
0651         unsigned int layer = TrkrDefs::getLayer(*iter);
0652         if (trkrid == TrkrDefs::mvtxId) {
0653           int stave_index = MvtxDefs::getStaveId(*iter);
0654           int chip_index = MvtxDefs::getChipId(*iter);
0655 
0656           double ladder_location[3] = { 0.0, 0.0, 0.0 };
0657           //not exactly sure where the cylinder geoms are currently objectified. check this 
0658           CylinderGeom_Mvtx *geom = (CylinderGeom_Mvtx*) _geom_container_maps->GetLayerGeom(layer);
0659           // returns the center of the sensor in world coordinates - used to get the ladder phi location
0660           geom->find_sensor_center(stave_index, 0, 0, chip_index, ladder_location);//the mvtx module and half stave are 0
0661           n.SetXYZ(ladder_location[0], ladder_location[1], 0);
0662           n.RotateZ(geom->get_stave_phi_tilt());
0663         } 
0664         else if (trkrid == TrkrDefs::inttId) {
0665           //this may bug but it looks ok for now
0666           CylinderGeomIntt* geom = (CylinderGeomIntt*) _geom_container_intt->GetLayerGeom(layer);
0667           double hit_location[3] = { 0.0, 0.0, 0.0 };
0668           geom->find_segment_center(InttDefs::getLadderZId(*iter),InttDefs::getLadderPhiId(*iter), hit_location);
0669 
0670           n.SetXYZ(hit_location[0], hit_location[1], 0);
0671           n.RotateZ(geom->get_strip_phi_tilt());
0672         }
0673       }//if use_ladder_geom
0674       PHGenFit::Measurement* meas = new PHGenFit::PlanarMeasurement(pos, n,radius*cluster->getPhiError(), cluster->getZError());
0675       measurements.push_back(meas);
0676     }//cluster loop
0677     genfit::AbsTrackRep* rep = new genfit::RKTrackRep(_primary_pid_guess);
0678     PHGenFit::Track* track(new PHGenFit::Track(rep, seed_pos, seed_mom, seed_cov));
0679     track->addMeasurements(measurements);
0680 
0681     if (_fitter->processTrack(track, false) != 0) {
0682       if (_verbosity >= 1)
0683         LogWarning("Track fitting failed");
0684       return NULL;
0685     }
0686     track->getGenFitTrack()->setMcTrackId(intrack->get_id());
0687     return track;
0688   }//valid vtx 
0689   else{
0690     cerr<<PHWHERE<<": invalid vertex"<<endl;
0691     return NULL;
0692   }
0693 
0694 }
0695 
0696 /*inspired by PHG4TrackKalmanFitter
0697 PHGenFit::Track* SVReco::MakeGenFitTrack(const PHGenFit::Track* intrack, const genfit::GFRaveVertex* invertex){
0698   if (!intrack){
0699     cerr << PHWHERE << " Input PHGenFit::Track is NULL!" << endl;
0700     return NULL;
0701   }
0702 
0703   if (_use_ladder_geom and !_geom_container_intt and !_geom_container_maps) {
0704     cout << PHWHERE << "No PHG4CylinderGeomContainer found!" << endl;
0705     return NULL;
0706   }
0707 
0708   // Create measurements
0709   std::vector<PHGenFit::Measurement*> measurements;
0710 
0711   //create space point measurement from vtx
0712   //TODO check that getNTracks is properly initialized 
0713   if (invertex and invertex->getNTracks() > 1) {
0714     TVector3 pos=invertex->getPos();
0715     TMatrixDSym cov = invertex->getCov();
0716     PHGenFit::Measurement* meas = new PHGenFit::SpacepointMeasurement(
0717           pos, cov);
0718     measurements.push_back(meas);
0719 
0720     //convert SvtxTrack to matricies
0721     TVector3 seed_pos = intrack->getGenFitTrack()->getPos();
0722     TVector3 seed_mom = intrack->get_mom(); //mom stands for momentum
0723     TMatrixDSym seed_cov = intrack->getGenFitTrack()->getCov();
0724     
0725     cout<<"Making track cluster measurments"<<endl;
0726     //make measurements from the track clusters
0727     for (auto iter = intrack->get_cluster_keys().begin(); iter != intrack->get_cluster_keys().end(); ++iter){
0728       //    unsigned int cluster_id = *iter;
0729       TrkrCluster* cluster = _clustermap->findCluster(*iter);
0730       if (!cluster) {
0731         LogError("No cluster Found!");
0732         continue;
0733       }
0734       float x = cluster->getPosition(0);
0735       float y = cluster->getPosition(1);
0736       float radius = sqrt(x*x+y*y);
0737       TVector3 pos(cluster->getPosition(0), cluster->getPosition(1), cluster->getPosition(2));
0738       seed_mom.SetPhi(pos.Phi());
0739       seed_mom.SetTheta(pos.Theta());
0740 
0741       TVector3 n(cluster->getPosition(0), cluster->getPosition(1), 0);
0742       //cout<<"Cluster with {"<<cluster->getPosition(0)<<','<<cluster->getPosition(0)<<"}\n";
0743 
0744       if (_use_ladder_geom){ //I don't understand this bool
0745         unsigned int trkrid = TrkrDefs::getTrkrId(*iter);
0746         unsigned int layer = TrkrDefs::getLayer(*iter);
0747         if (trkrid == TrkrDefs::mvtxId) {
0748           int stave_index = MvtxDefs::getStaveId(*iter);
0749           int chip_index = MvtxDefs::getChipId(*iter);
0750 
0751           double ladder_location[3] = { 0.0, 0.0, 0.0 };
0752           //not exactly sure where the cylinder geoms are currently objectified. check this 
0753           CylinderGeom_Mvtx *geom = (CylinderGeom_Mvtx*) _geom_container_maps->GetLayerGeom(layer);
0754           // returns the center of the sensor in world coordinates - used to get the ladder phi location
0755           geom->find_sensor_center(stave_index, 0, 0, chip_index, ladder_location);//the mvtx module and half stave are 0
0756           n.SetXYZ(ladder_location[0], ladder_location[1], 0);
0757           n.RotateZ(geom->get_stave_phi_tilt());
0758         } 
0759         else if (trkrid == TrkrDefs::inttId) {
0760           //this may bug but it looks ok for now
0761           CylinderGeomIntt* geom = (CylinderGeomIntt*) _geom_container_intt->GetLayerGeom(layer);
0762           double hit_location[3] = { 0.0, 0.0, 0.0 };
0763           geom->find_segment_center(InttDefs::getLadderZId(*iter),InttDefs::getLadderPhiId(*iter), hit_location);
0764 
0765           n.SetXYZ(hit_location[0], hit_location[1], 0);
0766           n.RotateZ(geom->get_strip_phi_tilt());
0767         }
0768       }//if use_ladder_geom
0769       PHGenFit::Measurement* meas = new PHGenFit::PlanarMeasurement(pos, n,radius*cluster->getPhiError(), cluster->getZError());
0770       measurements.push_back(meas);
0771     }//cluster loop
0772     genfit::AbsTrackRep* rep = new genfit::RKTrackRep(_primary_pid_guess);
0773     PHGenFit::Track* track(new PHGenFit::Track(rep, seed_pos, seed_mom, seed_cov));
0774     track->addMeasurements(measurements);
0775 
0776     if (_fitter->processTrack(track, false) != 0) {
0777       if (_verbosity >= 1)
0778         LogWarning("Track fitting failed");
0779       return NULL;
0780     }
0781     track->getGenFitTrack()->setMcTrackId(intrack->get_id());
0782     return track;
0783   }//valid vtx 
0784   else{
0785     cerr<<PHWHERE<<": invalid vertex"<<endl;
0786     return NULL;
0787   }
0788 
0789 }*/
0790 
0791 //From PHG4TrackKalmanFitter
0792 SvtxVertex* SVReco::GFRVVtxToSvtxVertex(genfit::GFRaveVertex* rave_vtx)const{
0793   SvtxVertex* svtx_vtx= new SvtxVertex_v1();
0794 
0795     svtx_vtx->set_chisq(rave_vtx->getChi2());
0796     svtx_vtx->set_ndof(rave_vtx->getNdf());
0797     svtx_vtx->set_position(0, rave_vtx->getPos().X());
0798     svtx_vtx->set_position(1, rave_vtx->getPos().Y());
0799     svtx_vtx->set_position(2, rave_vtx->getPos().Z());
0800 
0801     for (int i = 0; i < 3; i++)
0802       for (int j = 0; j < 3; j++)
0803         svtx_vtx->set_error(i, j, rave_vtx->getCov()[i][j]);
0804 
0805     for (unsigned int i = 0; i < rave_vtx->getNTracks(); i++) {
0806       //TODO Assume id's are sync'ed between _trackmap_refit and gf_tracks, need to change?
0807       const genfit::Track* rave_track =
0808           rave_vtx->getParameters(i)->getTrack();
0809       for (unsigned int j = 0; j < _main_rf_phgf_tracks.size(); j++) {
0810         if (rave_track == _main_rf_phgf_tracks[j]->getGenFitTrack())
0811           svtx_vtx->insert_track(j);
0812       }
0813     }
0814     return svtx_vtx;
0815 }
0816 
0817 /*
0818  * Fill SvtxVertexMap from GFRaveVertexes and Tracks
0819  */
0820 void SVReco::FillVertexMap(
0821         const std::vector<genfit::GFRaveVertex*>& rave_vertices,
0822         const std::vector<genfit::Track*>& gf_tracks){
0823 
0824     for (unsigned int ivtx=0; ivtx<rave_vertices.size(); ++ivtx){
0825         genfit::GFRaveVertex* rave_vtx = rave_vertices[ivtx];
0826         rv_prim_vtx_ntrk = rave_vtx->getNTracks();
0827 
0828         //TVector3 vertex_position(rv_prim_vtx[0], rv_prim_vtx[1], rv_prim_vtx[2]);
0829 
0830         //cout << "N TRK gf: " << gf_tracks.size() << ", rv: " << rv_prim_vtx_ntrk << endl;
0831 
0832         for (int itrk=0; itrk< rv_prim_vtx_ntrk; itrk++){
0833 
0834             TVector3 rvtrk_mom = rave_vtx->getParameters(itrk)->getMom();
0835             float rvtrk_w = rave_vtx->getParameters(itrk)->getWeight();
0836 
0837             unsigned int rvtrk_mc_id = rave_vtx->getParameters(itrk)->getTrack()->getMcTrackId();
0838             _svtxtrk_wt_map[rvtrk_mc_id] = rvtrk_w;
0839 
0840             //cout << "w: " << rvtrk_w << ", mc id: " << rvtrk_mc_id << endl;
0841             /*
0842                  SvtxTrack* svtx_track = _trackmap->get(rvtrk_mc_id);
0843 
0844                  cout << "rave trk, px: " << rvtrk_mom.Px() << ", py: " << rvtrk_mom.Py() << ", pz: " << rvtrk_mom.Pz() << endl;
0845                  cout << "svtx trk, px: " << svtx_track->get_px() << ", py: " << svtx_track->get_py() << ", pz: " << svtx_track->get_pz() << endl;
0846                  */
0847 
0848             /*
0849                  for (SvtxTrackMap::ConstIter iter3=_trackmap->begin(); iter3!=_trackmap->end(); iter3++)
0850                  {
0851                  SvtxTrack* svtx_track = iter3->second;
0852 
0853                  if ( fabs((svtx_track->get_px()-rvtrk_mom.Px())/svtx_track->get_px())<0.10
0854                  && fabs((svtx_track->get_py()-rvtrk_mom.Py())/svtx_track->get_py())<0.10
0855                  && fabs((svtx_track->get_pz()-rvtrk_mom.Pz())/svtx_track->get_pz())<0.10 ){
0856                  cout << "rave trk, px: " << rvtrk_mom.Px() << ", py: " << rvtrk_mom.Py() << ", pz: " << rvtrk_mom.Pz() << endl;
0857             //cout << "ggff trk, px: " << gftrk_mom.Px() << ", py: " << gftrk_mom.Py() << ", pz: " << gftrk_mom.Pz() << endl;
0858             cout << "svtx trk, px: " << svtx_track->get_px() << ", py: " << svtx_track->get_py() << ", pz: " << svtx_track->get_pz() << endl;
0859             }
0860             }//iter3
0861             */
0862 
0863             //unsigned int trk_id = svtxtrk_id[itrk];
0864 
0865             /*
0866                  cout << "W: " << w_trk 
0867                  << ", id: " << rave_vtx->getParameters(itrk)->GetUniqueID()
0868                  << ", px: " << rave_vtx->getParameters(itrk)->getMom().Px() 
0869                  << ", py: " << rave_vtx->getParameters(itrk)->getMom().Py() 
0870                  << ", pz: " << rave_vtx->getParameters(itrk)->getMom().Pz() 
0871                  << endl;
0872                  */
0873 
0874             //if (svtxtrk_gftrk_map.find(svtx_track->get_id())!=svtxtrk_gftrk_map.end()){
0875             //}
0876 
0877             //TVector3 mom_trk = rave_vtx->getParameters(itrk)->getMom();
0878         }//primvtx_tracks loop
0879     }//rave verticies loop
0880     return;
0881 }
0882 
0883 int SVReco::GetSVMass_mom(
0884         const genfit::GFRaveVertex* rave_vtx,
0885         float & vtx_mass,
0886         float & vtx_px,
0887         float & vtx_py,
0888         float & vtx_pz,
0889         int & ntrk_good_pv
0890         ){
0891 
0892     float sum_E = 0, sum_px = 0, sum_py = 0, sum_pz = 0;
0893 
0894     int N_good = 0, N_good_pv = 0;
0895 
0896     for (unsigned int itrk=0; itrk<rave_vtx->getNTracks(); itrk++){
0897         TVector3 mom_trk = rave_vtx->getParameters(itrk)->getMom(); 
0898 
0899         double w_trk = rave_vtx->getParameters(itrk)->getWeight();
0900 
0901         sum_px += mom_trk.X();
0902         sum_py += mom_trk.Y();
0903         sum_pz += mom_trk.Z();
0904         sum_E += sqrt(mom_trk.Mag2() + 0.140*0.140);
0905 
0906         //cout << "W: " << w_trk << endl;
0907         if ( w_trk>0.7 ){
0908             N_good++;
0909 
0910             unsigned int rvtrk_mc_id = rave_vtx->getParameters(itrk)->getTrack()->getMcTrackId();
0911             //cout << "mc_id: " << rvtrk_mc_id << ", wt: " << svtxtrk_wt_map[rvtrk_mc_id] << endl;
0912             if ( _svtxtrk_wt_map[rvtrk_mc_id]>0.7 ){
0913                 N_good_pv++;
0914             }
0915         }//
0916 
0917     }//itrk
0918 
0919     vtx_mass =  sqrt(sum_E*sum_E - sum_px*sum_px - sum_py*sum_py - sum_pz*sum_pz);
0920     vtx_px = sum_px;
0921     vtx_py = sum_py;
0922     vtx_pz = sum_pz;
0923 
0924     ntrk_good_pv = N_good_pv;
0925 
0926     //cout << "Mass: " << vtx_mass << ", pT: " << vtx_pT << endl;
0927     return N_good;
0928 }
0929 
0930 //Seems deprecated
0931 void SVReco::FillSVMap(
0932         const std::vector<genfit::GFRaveVertex*>& rave_vertices,
0933         const std::vector<genfit::Track*>& gf_tracks){
0934 
0935     return;
0936 }
0937 
0938 PHGenFit::Track*  SVReco::refitTrack(SvtxVertex* vtx, SvtxTrack* svtxtrk){
0939   PHGenFit::Track* gftrk = MakeGenFitTrack(svtxtrk,vtx);
0940   if(gftrk) {
0941     cout<<"good genfit refit"<<endl;
0942     //MakeSvtxTrack(svtxtrk,gftrk,vtx);
0943     return gftrk;
0944   }
0945   else {
0946     cout<<"No refit possible"<<endl;
0947     return NULL;
0948   }
0949 }
0950 
0951 PHGenFit::Track*  SVReco::refitTrack(genfit::GFRaveVertex* vtx, SvtxTrack* svtxtrk){
0952   PHGenFit::Track* gftrk = MakeGenFitTrack(svtxtrk,vtx);
0953   if(gftrk) {
0954     cout<<"good genfit refit"<<endl;
0955     //MakeSvtxTrack(svtxtrk,gftrk,vtx);
0956     return gftrk;
0957   }
0958   else {
0959     cout<<"No refit possible"<<endl;
0960     return NULL;
0961   }
0962 }
0963 
0964 /*PHGenFit::Track*  SVReco::refitTrack(genfit::GFRaveVertex* vtx, PHGenFit::Track* phgf_track){
0965   PHGenFit::Track* gftrk = MakeGenFitTrack(phgf_track,vtx);
0966   if(gftrk) {
0967     cout<<"good genfit refit"<<endl;
0968     //MakeSvtxTrack(svtxtrk,gftrk,vtx);
0969     return gftrk;
0970   }
0971   else {
0972     cout<<"No refit possible"<<endl;
0973     return NULL;
0974   }
0975 }*/
0976 
0977 
0978 /*FIXME this code is broken I have made zero attempt to find out why
0979 * I decided not to use the SvtxTrack after they are refit
0980 * may need to make the phgf_track pointer shared 
0981 */
0982 std::shared_ptr<SvtxTrack> SVReco::MakeSvtxTrack(const SvtxTrack* svtx_track,
0983         const PHGenFit::Track* phgf_track, const SvtxVertex* vertex) {
0984 
0985     double chi2 = phgf_track->get_chi2();
0986     double ndf = phgf_track->get_ndf();
0987 
0988     TVector3 vertex_position(0, 0, 0);
0989     TMatrixF vertex_cov(3,3);
0990     double dvr2 = 0;
0991     double dvz2 = 0;
0992 
0993     if (vertex) {
0994         vertex_position.SetXYZ(vertex->get_x(), vertex->get_y(),
0995                 vertex->get_z());
0996         dvr2 = vertex->get_error(0, 0) + vertex->get_error(1, 1);
0997         dvz2 = vertex->get_error(2, 2);
0998 
0999         for(int i=0;i<3;i++)
1000             for(int j=0;j<3;j++)
1001                 vertex_cov[i][j] = vertex->get_error(i,j);
1002     }
1003     else{
1004         cerr<<PHWHERE<<": No vertex to make SvtxTrack"<<endl;
1005     }
1006 
1007     //genfit::MeasuredStateOnPlane* gf_state_beam_line_ca = nullptr;
1008     std::shared_ptr<genfit::MeasuredStateOnPlane> gf_state_beam_line_ca = nullptr;
1009     try {
1010         gf_state_beam_line_ca = std::shared_ptr<genfit::MeasuredStateOnPlane>(phgf_track->extrapolateToLine(vertex_position,
1011                     TVector3(0., 0., 1.)));
1012     } catch (...) {
1013         if (_verbosity >= 2)
1014             LogWarning("extrapolateToLine failed!");
1015     }
1016     if(!gf_state_beam_line_ca) return nullptr;
1017 
1018     /*!
1019      *  1/p, u'/z', v'/z', u, v
1020      *  u is defined as momentum X beam line at POCA of the beam line
1021      *  v is alone the beam line
1022      *  so u is the dca2d direction
1023      */
1024 
1025     double u = gf_state_beam_line_ca->getState()[3];
1026     double v = gf_state_beam_line_ca->getState()[4];
1027 
1028     double du2 = gf_state_beam_line_ca->getCov()[3][3];
1029     double dv2 = gf_state_beam_line_ca->getCov()[4][4];
1030 
1031     //delete gf_state_beam_line_ca;
1032 
1033     //const SvtxTrack_v1* temp_track = static_cast<const SvtxTrack_v1*> (svtx_track);
1034     //  SvtxTrack_v1* out_track = new SvtxTrack_v1(
1035     //      *static_cast<const SvtxTrack_v1*>(svtx_track));
1036     std::shared_ptr < SvtxTrack_v1 > out_track = std::shared_ptr < SvtxTrack_v1
1037         > (new SvtxTrack_v1(*static_cast<const SvtxTrack_v1*>(svtx_track)));
1038 
1039     out_track->set_dca2d(u);
1040     out_track->set_dca2d_error(sqrt(du2 + dvr2));
1041 
1042     std::shared_ptr<genfit::MeasuredStateOnPlane> gf_state_vertex_ca = nullptr;
1043     try {
1044         gf_state_vertex_ca = std::shared_ptr < genfit::MeasuredStateOnPlane
1045             > (phgf_track->extrapolateToPoint(vertex_position));
1046     } catch (...) {
1047         if (_verbosity >= 2)
1048             LogWarning("extrapolateToPoint failed!");
1049     }
1050     if (!gf_state_vertex_ca) {
1051         //delete out_track;
1052         return nullptr;
1053     }
1054 
1055     TVector3 mom = gf_state_vertex_ca->getMom();
1056     TVector3 pos = gf_state_vertex_ca->getPos();
1057     TMatrixDSym cov = gf_state_vertex_ca->get6DCov();
1058 
1059     //  genfit::MeasuredStateOnPlane* gf_state_vertex_ca =
1060     //      phgf_track->extrapolateToLine(vertex_position,
1061     //          TVector3(0., 0., 1.));
1062 
1063     u = gf_state_vertex_ca->getState()[3];
1064     v = gf_state_vertex_ca->getState()[4];
1065 
1066     du2 = gf_state_vertex_ca->getCov()[3][3];
1067     dv2 = gf_state_vertex_ca->getCov()[4][4];
1068 
1069 
1070     double dca3d = sqrt(u * u + v * v);
1071     double dca3d_error = sqrt(du2 + dv2 + dvr2 + dvz2);
1072 
1073     out_track->set_dca(dca3d);
1074     out_track->set_dca_error(dca3d_error);
1075 
1076     /*!
1077      * dca3d_xy, dca3d_z
1078      */
1079 
1080 
1081     /*
1082     // Rotate from u,v,n to r: n X Z, Z': n X r, n using 5D state/cov
1083     // commented on 2017-10-09
1084     TMatrixF pos_in(3,1);
1085     TMatrixF cov_in(3,3);
1086     pos_in[0][0] = gf_state_vertex_ca->getState()[3];
1087     pos_in[1][0] = gf_state_vertex_ca->getState()[4];
1088     pos_in[2][0] = 0.;
1089     cov_in[0][0] = gf_state_vertex_ca->getCov()[3][3];
1090     cov_in[0][1] = gf_state_vertex_ca->getCov()[3][4];
1091     cov_in[0][2] = 0.;
1092     cov_in[1][0] = gf_state_vertex_ca->getCov()[4][3];
1093     cov_in[1][1] = gf_state_vertex_ca->getCov()[4][4];
1094     cov_in[1][2] = 0.;
1095     cov_in[2][0] = 0.;
1096     cov_in[2][1] = 0.;
1097     cov_in[2][2] = 0.;
1098     TMatrixF pos_out(3,1);
1099     TMatrixF cov_out(3,3);
1100     TVector3 vu = gf_state_vertex_ca->getPlane().get()->getU();
1101     TVector3 vv = gf_state_vertex_ca->getPlane().get()->getV();
1102     TVector3 vn = vu.Cross(vv);
1103     pos_cov_uvn_to_rz(vu, vv, vn, pos_in, cov_in, pos_out, cov_out);
1104     //! vertex cov in (u',v',n')
1105     TMatrixF vertex_cov_out(3,3);
1106     get_vertex_error_uvn(vu,vv,vn, vertex_cov, vertex_cov_out);
1107     float dca3d_xy = pos_out[0][0];
1108     float dca3d_z  = pos_out[1][0];
1109     float dca3d_xy_error = sqrt(cov_out[0][0] + vertex_cov_out[0][0]);
1110     float dca3d_z_error  = sqrt(cov_out[1][1] + vertex_cov_out[1][1]);
1111     //Begin DEBUG
1112     //  LogDebug("rotation debug---------- ");
1113     //  gf_state_vertex_ca->Print();
1114     //  LogDebug("dca rotation---------- ");
1115     //  pos_out = pos_in;
1116     //  cov_out = cov_in;
1117     //  pos_in.Print();
1118     //  cov_in.Print();
1119     //  pos_out.Print();
1120     //  cov_out.Print();
1121     //  cout
1122     //    <<"dca3d_xy: "<<dca3d_xy <<" +- "<<dca3d_xy_error*dca3d_xy_error
1123     //    <<"; dca3d_z: "<<dca3d_z<<" +- "<< dca3d_z_error*dca3d_z_error
1124     //    <<"\n";
1125     //  gf_state_vertex_ca->get6DCov().Print();
1126     //  LogDebug("vertex rotation---------- ");
1127     //  vertex_position.Print();
1128     //  vertex_cov.Print();
1129     //  vertex_cov_out.Print();
1130     //End DEBUG
1131     */
1132 
1133     //
1134     // in: X, Y, Z; out; r: n X Z, Z X r, Z
1135 
1136     float dca3d_xy = NAN;
1137     float dca3d_z  = NAN;
1138     float dca3d_xy_error = NAN;
1139     float dca3d_z_error  = NAN;
1140 
1141     try{
1142         TMatrixF pos_in(3,1);
1143         TMatrixF cov_in(3,3);
1144         TMatrixF pos_out(3,1);
1145         TMatrixF cov_out(3,3);
1146 
1147         TVectorD state6(6); // pos(3), mom(3)
1148         TMatrixDSym cov6(6,6); //
1149 
1150         gf_state_vertex_ca->get6DStateCov(state6, cov6);
1151 
1152         TVector3 vn(state6[3], state6[4], state6[5]);
1153 
1154         // mean of two multivariate gaussians Pos - Vertex
1155         pos_in[0][0] = state6[0] - vertex_position.X();
1156         pos_in[1][0] = state6[1] - vertex_position.Y();
1157         pos_in[2][0] = state6[2] - vertex_position.Z();
1158 
1159 
1160         for(int i=0;i<3;++i){
1161             for(int j=0;j<3;++j){
1162                 cov_in[i][j] = cov6[i][j] + vertex_cov[i][j];
1163             }
1164         }
1165 
1166         dca3d_xy = pos_out[0][0];
1167         dca3d_z  = pos_out[2][0];
1168         dca3d_xy_error = sqrt(cov_out[0][0]);
1169         dca3d_z_error  = sqrt(cov_out[2][2]);
1170 
1171 #ifdef _DEBUG_
1172         cout<<__LINE__<<": Vertex: ----------------"<<endl;
1173         vertex_position.Print();
1174         vertex_cov.Print();
1175 
1176         cout<<__LINE__<<": State: ----------------"<<endl;
1177         state6.Print();
1178         cov6.Print();
1179 
1180         cout<<__LINE__<<": Mean: ----------------"<<endl;
1181         pos_in.Print();
1182         cout<<"===>"<<endl;
1183         pos_out.Print();
1184 
1185         cout<<__LINE__<<": Cov: ----------------"<<endl;
1186         cov_in.Print();
1187         cout<<"===>"<<endl;
1188         cov_out.Print();
1189 
1190         cout<<endl;
1191 #endif
1192 
1193     } catch (...) {
1194         if (_verbosity > 0)
1195             LogWarning("DCA calculationfailed!");
1196     }
1197 
1198     out_track->set_dca3d_xy(dca3d_xy);
1199     out_track->set_dca3d_z(dca3d_z);
1200     out_track->set_dca3d_xy_error(dca3d_xy_error);
1201     out_track->set_dca3d_z_error(dca3d_z_error);
1202 
1203     //if(gf_state_vertex_ca) delete gf_state_vertex_ca;
1204 
1205     out_track->set_chisq(chi2);
1206     out_track->set_ndf(ndf);
1207     out_track->set_charge(phgf_track->get_charge());
1208 
1209     out_track->set_px(mom.Px());
1210     out_track->set_py(mom.Py());
1211     out_track->set_pz(mom.Pz());
1212 
1213     out_track->set_x(pos.X());
1214     out_track->set_y(pos.Y());
1215     out_track->set_z(pos.Z());
1216 
1217     for (int i = 0; i < 6; i++) {
1218         for (int j = i; j < 6; j++) {
1219             out_track->set_error(i, j, cov[i][j]);
1220         }
1221     }
1222 
1223     //  for (SvtxTrack::ConstClusterIter iter = svtx_track->begin_clusters();
1224     //      iter != svtx_track->end_clusters(); ++iter) {
1225     //    unsigned int cluster_id = *iter;
1226     //    SvtxCluster* cluster = _clustermap->get(cluster_id);
1227     //    if (!cluster) {
1228     //      LogError("No cluster Found!");
1229     //      continue;
1230     //    }
1231     //    //cluster->identify(); //DEBUG
1232     //
1233     //    //unsigned int l = cluster->get_layer();
1234     //
1235     //    TVector3 pos(cluster->get_x(), cluster->get_y(), cluster->get_z());
1236     //
1237     //    double radius = pos.Pt();
1238     //
1239     //    std::shared_ptr<genfit::MeasuredStateOnPlane> gf_state = nullptr;
1240     //    try {
1241     //      gf_state = std::shared_ptr < genfit::MeasuredStateOnPlane
1242     //          > (phgf_track->extrapolateToCylinder(radius,
1243     //              TVector3(0, 0, 0), TVector3(0, 0, 1), 0));
1244     //    } catch (...) {
1245     //      if (Verbosity() >= 2)
1246     //        LogWarning("Exrapolation failed!");
1247     //    }
1248     //    if (!gf_state) {
1249     //      if (Verbosity() > 1)
1250     //        LogWarning("Exrapolation failed!");
1251     //      continue;
1252     //    }
1253     //
1254     //    //SvtxTrackState* state = new SvtxTrackState_v1(radius);
1255     //    std::shared_ptr<SvtxTrackState> state = std::shared_ptr<SvtxTrackState> (new SvtxTrackState_v1(radius));
1256     //    state->set_x(gf_state->getPos().x());
1257     //    state->set_y(gf_state->getPos().y());
1258     //    state->set_z(gf_state->getPos().z());
1259     //
1260     //    state->set_px(gf_state->getMom().x());
1261     //    state->set_py(gf_state->getMom().y());
1262     //    state->set_pz(gf_state->getMom().z());
1263     //
1264     //    //gf_state->getCov().Print();
1265     //
1266     //    for (int i = 0; i < 6; i++) {
1267     //      for (int j = i; j < 6; j++) {
1268     //        state->set_error(i, j, gf_state->get6DCov()[i][j]);
1269     //      }
1270     //    }
1271     //
1272     //    out_track->insert_state(state.get());
1273     //
1274     //#ifdef _DEBUG_
1275     //    cout
1276     //    <<__LINE__
1277     //    <<": " << radius <<" => "
1278     //    <<sqrt(state->get_x()*state->get_x() + state->get_y()*state->get_y())
1279     //    <<endl;
1280     //#endif
1281     //  }
1282 
1283 #ifdef _DEBUG_
1284     cout << __LINE__ << endl;
1285 #endif
1286 
1287     const genfit::Track *gftrack = phgf_track->getGenFitTrack();
1288     const genfit::AbsTrackRep *rep = gftrack->getCardinalRep();
1289     for(unsigned int id = 0; id< gftrack->getNumPointsWithMeasurement();++id) {
1290         genfit::TrackPoint *trpoint = gftrack->getPointWithMeasurementAndFitterInfo(id, gftrack->getCardinalRep());
1291 
1292         if(!trpoint) {
1293             if (_verbosity > 1)
1294                 LogWarning("!trpoint");
1295             continue;
1296         }
1297 
1298         genfit::KalmanFitterInfo* kfi = static_cast<genfit::KalmanFitterInfo*>( trpoint->getFitterInfo(rep) );
1299         if(!kfi) {
1300             if (_verbosity > 1)
1301                 LogWarning("!kfi");
1302             continue;
1303         }
1304 
1305         std::shared_ptr<const genfit::MeasuredStateOnPlane> gf_state = nullptr;
1306         try {
1307             //gf_state = std::shared_ptr <genfit::MeasuredStateOnPlane> (const_cast<genfit::MeasuredStateOnPlane*> (&(kfi->getFittedState(true))));
1308             const genfit::MeasuredStateOnPlane* temp_state = &(kfi->getFittedState(true));
1309             gf_state = std::shared_ptr <genfit::MeasuredStateOnPlane> (new genfit::MeasuredStateOnPlane(*temp_state));
1310         } catch (...) {
1311             if (_verbosity > 1)
1312                 LogWarning("Exrapolation failed!");
1313         }
1314         if (!gf_state) {
1315             if (_verbosity > 1)
1316                 LogWarning("Exrapolation failed!");
1317             continue;
1318         }
1319         genfit::MeasuredStateOnPlane temp;
1320         float pathlength = -phgf_track->extrapolateToPoint(temp,vertex_position,id);
1321 
1322         std::shared_ptr<SvtxTrackState> state = std::shared_ptr<SvtxTrackState> (new SvtxTrackState_v1(pathlength));
1323         state->set_x(gf_state->getPos().x());
1324         state->set_y(gf_state->getPos().y());
1325         state->set_z(gf_state->getPos().z());
1326 
1327         state->set_px(gf_state->getMom().x());
1328         state->set_py(gf_state->getMom().y());
1329         state->set_pz(gf_state->getMom().z());
1330 
1331         //gf_state->getCov().Print();
1332 
1333         for (int i = 0; i < 6; i++) {
1334             for (int j = i; j < 6; j++) {
1335                 state->set_error(i, j, gf_state->get6DCov()[i][j]);
1336             }
1337         }
1338 
1339         out_track->insert_state(state.get());
1340 
1341 #ifdef _DEBUG_
1342         cout
1343             <<__LINE__
1344             <<": " << id
1345             <<": " << pathlength <<" => "
1346             <<sqrt(state->get_x()*state->get_x() + state->get_y()*state->get_y())
1347             <<endl;
1348 #endif
1349     }
1350 
1351     return out_track;
1352 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!