HF-Jet/HFDiJetMomImbalance/BDiJetModule.C

0001 /*!
0002  *  \file   BDiJetModule.C
0003  *  \brief    Gether information from Truth Jet and Tracks, DST -> NTuple
0004  *  \details  Gether information from Truth Jet and Tracks, DST -> NTuple
0005  *  \author   Dennis V. Perepelitsa
0006  */
0007
0008 #include "BDiJetModule.h"
0009
0010 #include <phool/getClass.h>
0011 #include <phool/phool.h>
0012 #include <phool/PHCompositeNode.h>
0013 #include <phool/PHIODataNode.h>
0014 #include <phool/PHNodeIterator.h>
0015
0016 #include <phgeom/PHGeomUtility.h>
0017 #include <phfield/PHFieldUtility.h>
0018
0019 #include <fun4all/Fun4AllServer.h>
0020 #include <fun4all/Fun4AllReturnCodes.h>
0021 #include <fun4all/PHTFileServer.h>
0022
0023 #include <g4main/PHG4TruthInfoContainer.h>
0024 #include <g4main/PHG4Particle.h>
0025 #include <g4main/PHG4Hit.h>
0026 #include <g4main/PHG4HitContainer.h>
0027
0028 #include <g4hough/SvtxCluster.h>
0029 #include <g4hough/SvtxClusterMap.h>
0030 #include <g4hough/SvtxHitMap.h>
0031 #include <g4hough/SvtxTrack.h>
0032 #include <g4hough/SvtxVertex.h>
0033 #include <g4hough/SvtxTrackMap.h>
0034 #include <g4hough/SvtxVertexMap.h>
0035
0036 #include <g4jets/JetMap.h>
0037 #include <g4jets/Jet.h>
0038
0039 #include <g4detectors/PHG4CellContainer.h>
0040 #include <g4detectors/PHG4CylinderGeomContainer.h>
0041 #include <g4detectors/PHG4Cell.h>
0042 #include <g4detectors/PHG4CylinderGeom_MAPS.h>
0043 #include <g4detectors/PHG4CylinderGeom_Siladders.h>
0044
0045 #include <phgenfit/Fitter.h>
0046 #include <phgenfit/PlanarMeasurement.h>
0047 #include <phgenfit/Track.h>
0048 #include <phgenfit/SpacepointMeasurement.h>
0049
0050 //GenFit
0051 #include <GenFit/FieldManager.h>
0052 #include <GenFit/GFRaveConverters.h>
0053 #include <GenFit/GFRaveVertex.h>
0054 #include <GenFit/GFRaveVertexFactory.h>
0055 #include <GenFit/MeasuredStateOnPlane.h>
0056 #include <GenFit/RKTrackRep.h>
0057 #include <GenFit/StateOnPlane.h>
0058 #include <GenFit/Track.h>
0059
0060 //Rave
0061 #include <rave/Version.h>
0062 #include <rave/Track.h>
0063 #include <rave/VertexFactory.h>
0064 #include <rave/ConstantMagneticField.h>
0065
0066 #include <phhepmc/PHHepMCGenEvent.h>
0067 #include <HepMC/GenEvent.h>
0068 #include <HepMC/GenVertex.h>
0069
0070 #include <HFJetTruthGeneration/HFJetDefs.h>
0071
0072 #include "TTree.h"
0073 #include "TFile.h"
0074 #include "TLorentzVector.h"
0075
0076
0077 #include <iostream>
0078 #include <map>
0079 #include <utility>
0080 #include <vector>
0081 #include <memory>
0082
0083
0084
0085 #define LogDebug(exp)   std::cout<<"DEBUG: "  <<__FILE__<<": "<<__LINE__<<": "<< exp <<"\n"
0086 #define LogError(exp)   std::cout<<"ERROR: "  <<__FILE__<<": "<<__LINE__<<": "<< exp <<"\n"
0087 #define LogWarning(exp) std::cout<<"WARNING: "  <<__FILE__<<": "<<__LINE__<<": "<< exp <<"\n"
0088
0089
0090 BDiJetModule::BDiJetModule(const std::string &name,
0091                            const std::string &ofName) :
0092   SubsysReco(name),
0093   _jetmap_truth(NULL),
0094   _clustermap(NULL),
0095   _trackmap(NULL),
0096   _vertexmap(NULL),
0097   _verbose(false),
0098   _write_tree(true),
0099   _ana_truth(true),
0100   _ana_reco(true),
0101   _do_evt_display(false),
0102   _use_ladder_geom(false),
0103   _cut_jet(true),
0104   _cut_Ncluster(false),
0105   _foutname(ofName),
0106   _f(NULL),
0107   _tree(NULL),
0108   _primary_pid_guess(211),
0109   _cut_min_pT(0.1),
0110   _cut_chi2_ndf(5.0),
0111   _cut_jet_pT(10.0),
0112   _cut_jet_eta(0.7),
0113   _cut_jet_R(0.4),
0114   _track_fitting_alg_name("DafRef"),
0115   _fitter(NULL),
0116   _vertexing_method("avf-smoothing:1"),
0117   _vertex_finder(NULL)
0118 {
0119
0120
0121 }
0122
0123 int BDiJetModule::Init(PHCompositeNode *topNode)
0124 {
0125
0126   //-- set up the tree
0127   if ( _write_tree )
0128   {
0129     ResetVariables();
0130     InitTree();
0131   }
0132
0133   return 0;
0134
0135 }
0136
0137 int BDiJetModule::InitRun(PHCompositeNode *topNode)
0138 {
0139
0140   TGeoManager* tgeo_manager = PHGeomUtility::GetTGeoManager(topNode);
0141   PHField * field = PHFieldUtility::GetFieldMapNode(nullptr, topNode);
0142
0143   _fitter = PHGenFit::Fitter::getInstance(tgeo_manager, field,
0144                                           _track_fitting_alg_name, "RKTrackRep", _do_evt_display);
0145
0146   if (!_fitter) {
0147     std::cerr << PHWHERE << std::endl;
0148     return Fun4AllReturnCodes::ABORTRUN;
0149   }
0150
0151   _vertex_finder = new genfit::GFRaveVertexFactory(verbosity);
0152   _vertex_finder->setMethod(_vertexing_method.data());
0153
0154   if (!_vertex_finder) {
0155     std::cerr << PHWHERE << std::endl;
0156     return Fun4AllReturnCodes::ABORTRUN;
0157   }
0158
0159
0160   //-- set counters
0161   _ievent = -1; // since we count at the beginning of the event, start at -1
0162   _b_event = -1;
0163
0164
0165   return 0;
0166
0167 }
0168
0169 int BDiJetModule::process_event(PHCompositeNode *topNode)
0170 {
0171
0172   // count event
0173   _ievent++;
0174   _b_event = _ievent;
0175
0176   //-- reset tree variables
0177   if ( _write_tree )
0178     ResetVariables();
0179
0180   //------
0181   //-- Get the nodes
0182   //------
0183   int retval = GetNodes(topNode);
0184   if ( retval != Fun4AllReturnCodes::EVENT_OK )
0185     return retval;
0186
0187   //------
0188   //-- Analyze truth jets
0189   //------
0190   if ( _ana_truth )
0191   {
0192     _truthjet_n = 0;
0193
0194     int ijet_t = 0;
0195     for (JetMap::Iter iter = _jetmap_truth->begin(); iter != _jetmap_truth->end(); ++iter)
0196     {
0197       Jet* this_jet = iter->second;
0198
0199       float this_pt = this_jet->get_pt();
0200       float this_phi = this_jet->get_phi();
0201       float this_eta = this_jet->get_eta();
0202
0203       // this_jet->get_property(static_cast<Jet::PROPERTY>(prop_JetPartonFlavor)) != 5)
0204       if (this_jet->get_pt() < _cut_jet_pT || fabs(this_eta) > _cut_jet_eta)
0205         continue;
0206
0207       if ( _write_tree )
0208       {
0209         _truthjet_parton_flavor[_truthjet_n] = this_jet->get_property(static_cast<Jet::PROPERTY>(prop_JetPartonFlavor));
0210         _truthjet_hadron_flavor[_truthjet_n] = this_jet->get_property(static_cast<Jet::PROPERTY>(prop_JetHadronFlavor));
0211         _truthjet_pt[_truthjet_n] = this_pt;
0212         _truthjet_phi[_truthjet_n] = this_phi;
0213         _truthjet_eta[_truthjet_n] = this_eta;
0214       }
0215
0216       if (_verbose)
0217         std::cout << " truth jet #" << ijet_t << ", pt / eta / phi = "
0218                   << this_pt << " / " << this_eta << " / " << this_phi
0219                   << std::endl;
0220
0221       _truthjet_n++;
0222     }
0223
0224     // only care if the event has a jet of interest
0225     if ( _truthjet_n < 1 )
0226       return Fun4AllReturnCodes::ABORTEVENT;
0227
0228   } // _ana_truth
0229
0230   //------
0231   //-- Analyze reconstructed information
0232   //------
0233   if ( _ana_reco )
0234   {
0235
0236     //! stands for Refit_GenFit_Tracks
0237     std::vector<genfit::Track*> rf_gf_tracks;
0238     rf_gf_tracks.clear();
0239     std::vector<PHGenFit::Track*> rf_phgf_tracks;
0240     rf_phgf_tracks.clear();
0241
0242     std::map<unsigned int, unsigned int> svtxtrk_gftrk_map;
0243     std::map<unsigned int, unsigned int> svtxtrk_nmvtx_map;
0244     std::map<unsigned int, unsigned int> svtxtrk_nintt_map;
0245
0246     SvtxVertex *vertex = _vertexmap->get(0);
0247
0248     for (SvtxTrackMap::Iter iter = _trackmap->begin();
0249          iter != _trackmap->end();
0250          ++iter)
0251     {
0252       SvtxTrack* svtx_track = iter->second;
0253       if (!svtx_track)
0254         continue;
0255       if (!(svtx_track->get_pt() > _cut_min_pT))
0256         continue;
0257       if ((svtx_track->get_chisq() / svtx_track->get_ndf()) > _cut_chi2_ndf)
0258         continue;
0259
0260       int n_MVTX = 0, n_INTT = 0;
0261       for (SvtxTrack::ConstClusterIter iter2 = svtx_track->begin_clusters();
0262            iter2 != svtx_track->end_clusters();
0263            iter2++)
0264       {
0265
0266         SvtxCluster* cluster = _clustermap->get(*iter2);
0267         unsigned int layer = cluster->get_layer();
0268
0269         if (layer < 3) n_MVTX++;
0270         else if (layer < 7) n_INTT++;
0271       } // iter2 over SvtxTrack::ConstClusterIter
0272
0273       //cout << "n MVTX: " << n_MVTX << ", n INTT: " << n_INTT << std::endl;
0274
0275       /*
0276       //if (n_MVTX<2 || n_INTT<2)
0277       if (n_MVTX<2)
0278         continue;
0279       */
0280
0281       PHGenFit::Track* rf_phgf_track = MakeGenFitTrack(topNode, svtx_track, vertex);
0282
0283       if (rf_phgf_track)
0284       {
0285         svtxtrk_gftrk_map.insert(std::pair<int, int>(svtx_track->get_id(), rf_phgf_tracks.size()));
0286         svtxtrk_nmvtx_map.insert(std::pair<int, int>(svtx_track->get_id(), n_MVTX));
0287         svtxtrk_nintt_map.insert(std::pair<int, int>(svtx_track->get_id(), n_INTT));
0288         rf_phgf_tracks.push_back(rf_phgf_track);
0289
0290         if ( _cut_Ncluster )
0291         {
0292           if ( n_MVTX >= 1 && n_INTT >= 2 )
0293             rf_gf_tracks.push_back(rf_phgf_track->getGenFitTrack());
0294         }
0295         else
0296         {
0297           rf_gf_tracks.push_back(rf_phgf_track->getGenFitTrack());
0298         }
0299
0300       }
0301     } // iter over SvtxTrackMap
0302
0303     //! find vertex using tracks
0304     std::vector<genfit::GFRaveVertex*> rave_vertices;
0305     rave_vertices.clear();
0306     //!
0307
0308     if ( _vertexing_method.compare("avr-smoothing:1") != 0 )
0309       _vertex_finder->setMethod(_vertexing_method.data());
0310     if (rf_gf_tracks.size() >= 2)
0311     {
0312       try {
0313         _vertex_finder->findVertices(&rave_vertices, rf_gf_tracks);
0314       } catch (...) {
0315         std::cout << PHWHERE << "GFRaveVertexFactory::findVertices failed!";
0316       }
0317     }
0318
0319     FillVertexMap(rave_vertices, rf_gf_tracks);
0320
0321     //! change the vertex finding method
0322     if ( _vertexing_method.compare("avr-smoothing:1") != 0 )
0323       _vertex_finder->setMethod("avr-smoothing:1");
0324     std::vector<genfit::GFRaveVertex*> rave_vertices_jet;
0325     rave_vertices_jet.clear();
0326
0327     //! scan jetmap
0328     for (JetMap::ConstIter iter = _jetmap_truth->begin(); iter != _jetmap_truth->end(); iter++)
0329     {
0330       Jet *jet = iter->second;
0331
0332       float jet_pT = jet->get_pt();
0333       float jet_eta = jet->get_eta();
0334       float jet_phi = jet->get_phi();
0335
0336       if ( jet_pT < _cut_jet_pT ) continue;
0337       if ( fabs(jet_eta) > _cut_jet_eta ) continue;
0338
0339       int counter_r10 = 0, counter_miss = 0;
0340
0341       std::vector<genfit::Track*> rf_gf_tracks_jet;
0342       rf_gf_tracks_jet.clear();
0343       // vector<genfit::Track*> rf_gf_tracks_jet_pT05;
0344       // rf_gf_tracks_jet_pT05.clear();
0345       // vector<genfit::Track*> rf_gf_tracks_jet_pT10;
0346       // rf_gf_tracks_jet_pT10.clear();
0347       // vector<genfit::Track*> rf_gf_tracks_jet_pT15;
0348       // rf_gf_tracks_jet_pT15.clear();
0349       // vector<genfit::Track*> rf_gf_tracks_jet_pT20;
0350       // rf_gf_tracks_jet_pT20.clear();
0351
0352       //cout << "JET ETA: " << jet_eta << ", JET PHI: " << jet_phi << std::endl;
0353
0354       for (SvtxTrackMap::ConstIter iter2 = _trackmap->begin(); iter2 != _trackmap->end(); iter2++)
0355       {
0356         SvtxTrack* svtx_track = iter2->second;
0357
0358         float trk_phi = svtx_track->get_phi();
0359         float trk_eta = svtx_track->get_eta();
0360         // float trk_pT = svtx_track->get_pt();
0361
0362         float sin_phi = sin(jet_phi - trk_phi);
0363         float cos_phi = cos(jet_phi - trk_phi);
0364         float dphi = atan2(sin_phi, cos_phi);
0365
0366         if (sqrt((jet_eta - trk_eta) * (jet_eta - trk_eta) + dphi * dphi) < 1.0)
0367         {
0368           counter_r10++;
0369         }
0370
0371         if (sqrt((jet_eta - trk_eta) * (jet_eta - trk_eta) + dphi * dphi) > _cut_jet_R) continue;
0372
0373         if (svtxtrk_gftrk_map.find(svtx_track->get_id()) != svtxtrk_gftrk_map.end()) {
0374           unsigned int trk_index = svtxtrk_gftrk_map[svtx_track->get_id()];
0375           unsigned int nclus_mvtx = svtxtrk_nmvtx_map[svtx_track->get_id()];
0376           unsigned int nclus_intt = svtxtrk_nintt_map[svtx_track->get_id()];
0377
0378           //cout << "NCLUS MVTX: " << nclus_mvtx << std::endl;
0379
0380           PHGenFit::Track* rf_phgf_track = rf_phgf_tracks.at(trk_index);
0381
0382           if ( _cut_Ncluster ) {
0383             if ( nclus_mvtx >= 1 && nclus_intt >= 2 ) {
0384               rf_gf_tracks_jet.push_back(rf_phgf_track->getGenFitTrack());
0385               // if (trk_pT > 0.5) rf_gf_tracks_jet_pT05.push_back(rf_phgf_track->getGenFitTrack());
0386               // if (trk_pT > 1.0) rf_gf_tracks_jet_pT10.push_back(rf_phgf_track->getGenFitTrack());
0387               // if (trk_pT > 1.5) rf_gf_tracks_jet_pT15.push_back(rf_phgf_track->getGenFitTrack());
0388               // if (trk_pT > 2.0) rf_gf_tracks_jet_pT20.push_back(rf_phgf_track->getGenFitTrack());
0389             }
0390           } else {
0391             rf_gf_tracks_jet.push_back(rf_phgf_track->getGenFitTrack());
0392             // if (trk_pT > 0.5) rf_gf_tracks_jet_pT05.push_back(rf_phgf_track->getGenFitTrack());
0393             // if (trk_pT > 1.0) rf_gf_tracks_jet_pT10.push_back(rf_phgf_track->getGenFitTrack());
0394             // if (trk_pT > 1.5) rf_gf_tracks_jet_pT15.push_back(rf_phgf_track->getGenFitTrack());
0395             // if (trk_pT > 2.0) rf_gf_tracks_jet_pT20.push_back(rf_phgf_track->getGenFitTrack());
0396           }
0397         } else {
0398           counter_miss++;
0399         }
0400       }//trackmap
0401
0402       //! SV reco
0403       //!
0404       if (rf_gf_tracks_jet.size() > 1) {
0405         try {
0406           _vertex_finder->findVertices(&rave_vertices_jet, rf_gf_tracks_jet);
0407         } catch (...) {
0408           std::cout << PHWHERE << "GFRaveVertexFactory::findVertices failed!";
0409         }
0410       }
0411
0412       //cout << "N MISS: " << counter_miss << std::endl;
0413       if ( verbosity > 0 )
0414       {
0415         std::cout << "JET PT: " << jet_pT
0416                   << ", N TRK: " << int(jet->size_comp())
0417                   << ", CUT10: " << counter_r10
0418                   << ", CUT04: " << rf_gf_tracks_jet.size()
0419                   << ", N VTX: " << rave_vertices_jet.size()
0420                   << std::endl;
0421       }
0422
0423       if ( _write_tree )
0424       {
0425         rv_sv_pT00_nvtx[rv_sv_njets] = rave_vertices_jet.size();
0426         for (int ivtx = 0; ivtx < int(rave_vertices_jet.size()); ivtx++) {
0427           genfit::GFRaveVertex* rave_vtx = rave_vertices_jet[ivtx];
0428           rv_sv_pT00_vtx_x[rv_sv_njets][ivtx] = rave_vtx->getPos().X();
0429           rv_sv_pT00_vtx_y[rv_sv_njets][ivtx] = rave_vtx->getPos().Y();
0430           rv_sv_pT00_vtx_z[rv_sv_njets][ivtx] = rave_vtx->getPos().Z();
0431
0432           rv_sv_pT00_vtx_ex[rv_sv_njets][ivtx] = sqrt(rave_vtx->getCov()[0][0]);
0433           rv_sv_pT00_vtx_ey[rv_sv_njets][ivtx] = sqrt(rave_vtx->getCov()[1][1]);
0434           rv_sv_pT00_vtx_ez[rv_sv_njets][ivtx] = sqrt(rave_vtx->getCov()[2][2]);
0435
0436           rv_sv_pT00_vtx_ntrk[rv_sv_njets][ivtx] = (int)rave_vtx->getNTracks();
0437
0438           float vtx_mass, vtx_px, vtx_py, vtx_pz;
0439           rv_sv_pT00_vtx_ntrk_good[rv_sv_njets][ivtx] = GetSVMass_mom(rave_vtx, vtx_mass, vtx_px, vtx_py, vtx_pz);
0440
0441           //cout << "N TRK: " << rv_sv_pT00_vtx_ntrk[rv_sv_njets][ivtx] << ", GOOD: " << rv_sv_pT00_vtx_ntrk_good[rv_sv_njets][ivtx] << std::endl;
0442
0443           TVector3 vec1(vtx_px, vtx_py, vtx_pz);
0444           TVector3 vec2(rv_sv_pT00_vtx_x[rv_sv_njets][ivtx] - rv_prim_vtx[0], rv_sv_pT00_vtx_y[rv_sv_njets][ivtx] - rv_prim_vtx[1], rv_sv_pT00_vtx_z[rv_sv_njets][ivtx] - rv_prim_vtx[2]);
0445           float theta = vec1.Angle(vec2);
0446           float A = vec1.Mag() * sin(theta);
0447           float vtx_mass_corr = sqrt(vtx_mass * vtx_mass + A * A) + A;
0448
0449           rv_sv_pT00_vtx_mass[rv_sv_njets][ivtx] = vtx_mass;
0450           rv_sv_pT00_vtx_mass_corr[rv_sv_njets][ivtx] = vtx_mass_corr;
0451           rv_sv_pT00_vtx_pT[rv_sv_njets][ivtx] = sqrt(vtx_px * vtx_px + vtx_py * vtx_py);
0452         }
0453       }
0454
0455       for (genfit::GFRaveVertex* vertex : rave_vertices_jet) {
0456         delete vertex;
0457       }
0458       rave_vertices_jet.clear();
0459
0460
0461       //! jet information
0462       rv_sv_jet_id[rv_sv_njets] = jet->get_id();
0463       rv_sv_jet_pT[rv_sv_njets] = jet_pT;
0464       rv_sv_jet_phi[rv_sv_njets] = jet_phi;
0465       if (jet->has_property(static_cast<Jet::PROPERTY>(prop_JetPartonFlavor))) {
0466         rv_sv_jet_prop[rv_sv_njets][0] = int(jet->get_property(static_cast<Jet::PROPERTY>(prop_JetPartonFlavor)));
0467         rv_sv_jet_prop[rv_sv_njets][1] = int(jet->get_property(static_cast<Jet::PROPERTY>(prop_JetHadronFlavor)));
0468       }
0469
0470       rv_sv_njets++;
0471
0472     }
0473
0474
0475
0476
0477   } // _ana_reco
0478
0479
0480   //-- Cleanup
0481
0482   // fill tree
0483   if ( _write_tree )
0484     _tree->Fill();
0485
0486   //-- End
0487   return 0;
0488
0489 }
0490
0491 int BDiJetModule::End(PHCompositeNode * topNode)
0492 {
0493   if ( _write_tree )
0494   {
0495     _f->Write();
0496     _f->Close();
0497   }
0498
0499   return 0;
0500 }
0501
0502 /*
0503  * Initialize the TTree
0504  */
0505 void BDiJetModule::InitTree()
0506 {
0507   _f = new TFile(_foutname.c_str(), "RECREATE");
0508
0509   _tree = new TTree("ttree", "a withered deciduous tree");
0510
0511   _tree->Branch("event", &_b_event, "event/I");
0512
0513   if ( _ana_truth )
0514   {
0515     _tree->Branch("truthjet_n", &_truthjet_n, "truthjet_n/I");
0516     _tree->Branch("truthjet_parton_flavor", _truthjet_parton_flavor, "truthjet_parton_flavor[truthjet_n]/I");
0517     _tree->Branch("truthjet_hadron_flavor", _truthjet_hadron_flavor, "truthjet_hadron_flavor[truthjet_n]/I");
0518     _tree->Branch("truthjet_pt", _truthjet_pt, "truthjet_pt[truthjet_n]/F");
0519     _tree->Branch("truthjet_eta", _truthjet_eta,
0520                   "truthjet_eta[truthjet_n]/F");
0521     _tree->Branch("truthjet_phi", _truthjet_phi,
0522                   "truthjet_phi[truthjet_n]/F");
0523   }
0524
0525   if ( _ana_reco )
0526   {
0527     _tree->Branch("gf_prim_vtx", gf_prim_vtx, "gf_prim_vtx[3]/F");
0528     _tree->Branch("gf_prim_vtx_err", gf_prim_vtx_err, "gf_prim_vtx_err[3]/F");
0529     _tree->Branch("gf_prim_vtx_ntrk", &gf_prim_vtx_ntrk, "gf_prim_vtx_ntrk/I");
0530     _tree->Branch("rv_prim_vtx", rv_prim_vtx, "rv_prim_vtx[3]/F");
0531     _tree->Branch("rv_prim_vtx_err", rv_prim_vtx_err, "rv_prim_vtx_err[3]/F");
0532     _tree->Branch("rv_prim_vtx_ntrk", &rv_prim_vtx_ntrk, "rv_prim_vtx_ntrk/I");
0533
0534     _tree->Branch("rv_sv_njets", &rv_sv_njets, "rv_sv_njets/I");
0535     _tree->Branch("rv_sv_jet_id", rv_sv_jet_id, "rv_sv_jet_id[rv_sv_njets]/I");
0536     _tree->Branch("rv_sv_jet_prop", rv_sv_jet_prop, "rv_sv_jet_prop[rv_sv_njets][2]/I");
0537     _tree->Branch("rv_sv_jet_pT", rv_sv_jet_pT, "rv_sv_jet_pT[rv_sv_njets]/F");
0538     _tree->Branch("rv_sv_jet_phi", rv_sv_jet_phi, "rv_sv_jet_phi[rv_sv_njets]/F");
0539
0540     _tree->Branch("rv_sv_pT00_nvtx", rv_sv_pT00_nvtx, "rv_sv_pT00_nvtx[rv_sv_njets]/I");
0541     _tree->Branch("rv_sv_pT00_vtx_x", rv_sv_pT00_vtx_x, "rv_sv_pT00_vtx_x[rv_sv_njets][30]/F");
0542     _tree->Branch("rv_sv_pT00_vtx_y", rv_sv_pT00_vtx_y, "rv_sv_pT00_vtx_y[rv_sv_njets][30]/F");
0543     _tree->Branch("rv_sv_pT00_vtx_z", rv_sv_pT00_vtx_z, "rv_sv_pT00_vtx_z[rv_sv_njets][30]/F");
0544     _tree->Branch("rv_sv_pT00_vtx_ex", rv_sv_pT00_vtx_ex, "rv_sv_pT00_vtx_ex[rv_sv_njets][30]/F");
0545     _tree->Branch("rv_sv_pT00_vtx_ey", rv_sv_pT00_vtx_ey, "rv_sv_pT00_vtx_ey[rv_sv_njets][30]/F");
0546     _tree->Branch("rv_sv_pT00_vtx_ez", rv_sv_pT00_vtx_ez, "rv_sv_pT00_vtx_ez[rv_sv_njets][30]/F");
0547     _tree->Branch("rv_sv_pT00_vtx_ntrk", rv_sv_pT00_vtx_ntrk, "rv_sv_pT00_vtx_ntrk[rv_sv_njets][30]/I");
0548     _tree->Branch("rv_sv_pT00_vtx_ntrk_good", rv_sv_pT00_vtx_ntrk_good, "rv_sv_pT00_vtx_ntrk_good[rv_sv_njets][30]/I");
0549     _tree->Branch("rv_sv_pT00_vtx_mass", rv_sv_pT00_vtx_mass, "rv_sv_pT00_vtx_mass[rv_sv_njets][30]/F");
0550     _tree->Branch("rv_sv_pT00_vtx_mass_corr", rv_sv_pT00_vtx_mass_corr, "rv_sv_pT00_vtx_mass_corr[rv_sv_njets][30]/F");
0551     _tree->Branch("rv_sv_pT00_vtx_pT", rv_sv_pT00_vtx_pT, "rv_sv_pT00_vtx_pT[rv_sv_njets][30]/F");
0552   }
0553
0554 }
0555
0556
0557 /*
0558  * Reset all the TTree variables
0559  * Should be called for each event
0560  */
0561 void BDiJetModule::ResetVariables()
0562 {
0563
0564   //-- Truth jet info
0565   _truthjet_n = 0;
0566   for (int n = 0; n < 10; n++)
0567   {
0568     _truthjet_parton_flavor[n] = -9999;
0569     _truthjet_hadron_flavor[n] = -9999;
0570
0571     _truthjet_pt[n] = -1;
0572     _truthjet_eta[n] = -1;
0573     _truthjet_phi[n] = -1;
0574   }
0575
0576   //-- Reco info
0577   gf_prim_vtx_ntrk = rv_prim_vtx_ntrk = 0;
0578   for (int i = 0; i < 3; i++)
0579   {
0580     gf_prim_vtx[i] = gf_prim_vtx_err[i] = -999;
0581     rv_prim_vtx[i] = rv_prim_vtx_err[i] = -999;
0582   }//i
0583
0584   rv_sv_njets = 0;
0585
0586   for (int ijet = 0; ijet < 10; ijet++)
0587   {
0588     rv_sv_jet_id[ijet] = -999;
0589     rv_sv_jet_prop[ijet][0] = -999;
0590     rv_sv_jet_prop[ijet][1] = -999;
0591     rv_sv_jet_pT[ijet] = -999;
0592     rv_sv_jet_phi[ijet] = -999;
0593
0594     rv_sv_pT00_nvtx[ijet] = 0;
0595
0596     for (int ivtx = 0; ivtx < 30; ivtx++) {
0597       rv_sv_pT00_vtx_ntrk[ijet][ivtx] = 0;
0598       rv_sv_pT00_vtx_ntrk_good[ijet][ivtx] = 0;
0599
0600       rv_sv_pT00_vtx_mass[ijet][ivtx] = -999;
0601       rv_sv_pT00_vtx_mass_corr[ijet][ivtx] = -999;
0602       rv_sv_pT00_vtx_pT[ijet][ivtx] = -999;
0603
0604       rv_sv_pT00_vtx_x[ijet][ivtx] = -999;
0605       rv_sv_pT00_vtx_y[ijet][ivtx] = -999;
0606       rv_sv_pT00_vtx_z[ijet][ivtx] = -999;
0607       rv_sv_pT00_vtx_ex[ijet][ivtx] = -999;
0608       rv_sv_pT00_vtx_ey[ijet][ivtx] = -999;
0609       rv_sv_pT00_vtx_ez[ijet][ivtx] = -999;
0610
0611     }//ivtx
0612   }//ijet
0613
0614 }
0615
0616 /*
0617  * Get all the necessary nodes from the node tree
0618  */
0619 int BDiJetModule::GetNodes(PHCompositeNode * topNode)
0620 {
0621
0622   // Truth jet node
0623   _jetmap_truth = findNode::getClass<JetMap>(topNode,
0624                   "AntiKt_Truth_r04");
0625   if ( !_jetmap_truth && _ana_truth && _ievent < 2)
0626   {
0627     std::cout << PHWHERE << " AntiKt_Truth_r04 node not found ... aborting" << std::endl;
0628     return Fun4AllReturnCodes::ABORTEVENT;
0629   }
0630
0631   // Input Svtx Clusters
0632   _clustermap = findNode::getClass<SvtxClusterMap>(topNode, "SvtxClusterMap");
0633   if (!_clustermap && _ana_reco && _ievent < 2) {
0634     std::cout << PHWHERE << " SvtxClusterMap node not found on node tree"
0635               << std::endl;
0636     return Fun4AllReturnCodes::ABORTEVENT;
0637   }
0638
0639   // Input Svtx Tracks
0640   _trackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0641   if (!_trackmap && _ana_reco && _ievent < 2) {
0642     std::cout << PHWHERE << " SvtxClusterMap node not found on node tree"
0643               << std::endl;
0644     return Fun4AllReturnCodes::ABORTEVENT;
0645   }
0646
0647   // Input Svtx Vertices
0648   _vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
0649   if (!_vertexmap && _ana_reco && _ievent < 2) {
0650     std::cout << PHWHERE << " SvtxVertexrMap node not found on node tree"
0651               << std::endl;
0652     return Fun4AllReturnCodes::ABORTEVENT;
0653   }
0654
0655
0656
0657   return Fun4AllReturnCodes::EVENT_OK;
0658
0659 }
0660
0661
0662 /*
0663  * Load reconstructed track into GenFit track for input to Rave
0664  */
0665 PHGenFit::Track* BDiJetModule::MakeGenFitTrack(PHCompositeNode * topNode,
0666     const SvtxTrack * intrack,
0667     const SvtxVertex * vertex)
0668 {
0669
0670
0671   if (!intrack) {
0672     std::cerr << PHWHERE << " Input SvtxTrack is NULL!" << std::endl;
0673     return NULL;
0674   }
0675
0676   SvtxHitMap* hitsmap = NULL;
0677   hitsmap = findNode::getClass<SvtxHitMap>(topNode, "SvtxHitMap");
0678   if (!hitsmap) {
0679     std::cout << PHWHERE << "ERROR: Can't find node SvtxHitMap" << std::endl;
0680     return NULL;
0681   }
0682
0683   PHG4CellContainer* cells_svtx = findNode::getClass<PHG4CellContainer>(topNode, "G4CELL_SVTX");
0684   PHG4CellContainer* cells_intt = findNode::getClass<PHG4CellContainer>(topNode, "G4CELL_SILICON_TRACKER");
0685   PHG4CellContainer* cells_maps = findNode::getClass<PHG4CellContainer>(topNode, "G4CELL_MAPS");
0686
0687   if (_use_ladder_geom and !cells_svtx and !cells_intt and !cells_maps) {
0688     std::cout << PHWHERE << "No PHG4CellContainer found!" << std::endl;
0689     return NULL;
0690   }
0691
0692   PHG4CylinderGeomContainer* geom_container_intt = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_SILICON_TRACKER");
0693   PHG4CylinderGeomContainer* geom_container_maps = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MAPS");
0694
0695   if (_use_ladder_geom and !geom_container_intt and !geom_container_maps) {
0696     std::cout << PHWHERE << "No PHG4CylinderGeomContainer found!" << std::endl;
0697     return NULL;
0698   }
0699
0700   // prepare seed
0701   TVector3 seed_mom(100, 0, 0);
0702   TVector3 seed_pos(0, 0, 0);
0703   TMatrixDSym seed_cov(6);
0704   for (int i = 0; i < 6; i++) {
0705     for (int j = 0; j < 6; j++) {
0706       seed_cov[i][j] = 100.;
0707     }
0708   }
0709
0710   /*
0711   TVector3 seed_pos(intrack->get_x(), intrack->get_y(), intrack->get_z());
0712   TVector3 seed_mom(intrack->get_px(), intrack->get_py(), intrack->get_pz());
0713   TMatrixDSym seed_cov(6);
0714   for (int i=0; i<6; i++){
0715     for (int j=0; j<6; j++){
0716       seed_cov[i][j] = intrack->get_error(i,j);
0717     }
0718   }
0719   */
0720
0721   // Create measurements
0722   std::vector<PHGenFit::Measurement*> measurements;
0723
0724
0725   std::map<float, unsigned int> m_r_cluster_id;
0726
0727   for (auto iter = intrack->begin_clusters(); iter != intrack->end_clusters(); ++iter) {
0728     unsigned int cluster_id = *iter;
0729     SvtxCluster* cluster = _clustermap->get(cluster_id);
0730     float x = cluster->get_x();
0731     float y = cluster->get_y();
0732     float r = sqrt(x * x + y * y);
0733     m_r_cluster_id.insert(std::pair<float, unsigned int>(r, cluster_id));
0734   }
0735
0736   for (auto iter = m_r_cluster_id.begin(); iter != m_r_cluster_id.end(); ++iter) {
0737     //for (SvtxTrack::ConstClusterIter iter = intrack->begin_clusters(); iter != intrack->end_clusters(); ++iter){
0738     unsigned int cluster_id = iter->second;
0739     //unsigned int cluster_id = *iter;
0740     SvtxCluster* cluster = _clustermap->get(cluster_id);
0741     if (!cluster) {
0742       LogError("No cluster Found!");
0743       continue;
0744     }
0745
0746     TVector3 pos(cluster->get_x(), cluster->get_y(), cluster->get_z());
0747
0748     seed_mom.SetPhi(pos.Phi());
0749     seed_mom.SetTheta(pos.Theta());
0750
0751     //TODO use u, v explicitly?
0752     TVector3 n(cluster->get_x(), cluster->get_y(), 0);
0753
0754     unsigned int begin_hit_id = *(cluster->begin_hits());
0755     SvtxHit* svtxhit = hitsmap->find(begin_hit_id)->second;
0756
0757     PHG4Cell* cell_svtx = nullptr;
0758     PHG4Cell* cell_intt = nullptr;
0759     PHG4Cell* cell_maps = nullptr;
0760
0761     if (_use_ladder_geom and cells_svtx) cell_svtx = cells_svtx->findCell(svtxhit->get_cellid());
0762     if (_use_ladder_geom and cells_intt) cell_intt = cells_intt->findCell(svtxhit->get_cellid());
0763     if (_use_ladder_geom and cells_maps) cell_maps = cells_maps->findCell(svtxhit->get_cellid());
0764     if (_use_ladder_geom and !(cell_svtx or cell_intt or cell_maps)) {
0765       if (verbosity >= 0)
0766         LogError("!(cell_svtx or cell_intt or cell_maps)");
0767       continue;
0768     }
0769
0770     //float phi_tilt[7] = {0.304, 0.304, 0.304, 0.244, 0.244, 0.209, 0.201};
0771     unsigned int layer = cluster->get_layer();
0772
0773     //NEW
0774     if (_use_ladder_geom) {
0775       if (cell_maps) {
0776         PHG4Cell* cell = cell_maps;
0777
0778         int stave_index = cell->get_stave_index();
0779         int half_stave_index = cell->get_half_stave_index();
0780         int module_index = cell->get_module_index();
0781         int chip_index = cell->get_chip_index();
0782
0783         double ladder_location[3] = { 0.0, 0.0, 0.0 };
0784         PHG4CylinderGeom_MAPS *geom = (PHG4CylinderGeom_MAPS*) geom_container_maps->GetLayerGeom(layer);
0785         // returns the center of the sensor in world coordinates - used to get the ladder phi location
0786         geom->find_sensor_center(stave_index, half_stave_index, module_index, chip_index, ladder_location);
0787         n.SetXYZ(ladder_location[0], ladder_location[1], 0);
0788         n.RotateZ(geom->get_stave_phi_tilt());
0789       } else if (cell_intt) {
0790         PHG4Cell* cell = cell_intt;
0791         PHG4CylinderGeom_Siladders* geom = (PHG4CylinderGeom_Siladders*) geom_container_intt->GetLayerGeom(layer);
0792         double hit_location[3] = { 0.0, 0.0, 0.0 };
0793         geom->find_segment_center(cell->get_ladder_z_index(), cell->get_ladder_phi_index(), hit_location);
0794
0795         n.SetXYZ(hit_location[0], hit_location[1], 0);
0796         n.RotateZ(geom->get_strip_phi_tilt());
0797       }
0798     }
0799
0800     /*
0801     //OLD
0802     if ((cells_maps and geom_container_maps) and layer < 3) {
0803       unsigned int begin_hit_id = *(cluster->begin_hits());
0804       SvtxHit* hit = hitsmap->find(begin_hit_id)->second;
0805       PHG4Cell* cell = cells_maps->findCell(hit->get_cellid());
0806       int stave_index = cell->get_stave_index();
0807       int half_stave_index = cell->get_half_stave_index();
0808       int module_index = cell->get_module_index();
0809       int chip_index = cell->get_chip_index();
0810
0811       double ladder_location[3] = {0.0, 0.0, 0.0};
0812       PHG4CylinderGeom_MAPS *geom = (PHG4CylinderGeom_MAPS*) geom_container_maps->GetLayerGeom(layer);
0813       geom->find_sensor_center(stave_index, half_stave_index, module_index, chip_index, ladder_location);
0814       n.SetXYZ(ladder_location[0], ladder_location[1], 0);
0815       n.RotateZ(phi_tilt[layer]);
0816     } else if ((cells_intt and geom_container_intt) and pos.Perp() < 30.) {
0817       unsigned int begin_hit_id = *(cluster->begin_hits());
0818       SvtxHit* hit = hitsmap->find(begin_hit_id)->second;
0819       PHG4Cell* cell = cells_intt->findCell(hit->get_cellid());
0820       PHG4CylinderGeom_Siladders* geom = (PHG4CylinderGeom_Siladders*) geom_container_intt->GetLayerGeom(layer);
0821       double hit_location[3] = { 0.0, 0.0, 0.0 };
0822       geom->find_segment_center(cell->get_ladder_z_index(),cell->get_ladder_phi_index(), hit_location);
0823
0824       n.SetXYZ(hit_location[0], hit_location[1], 0);
0825       n.RotateZ(phi_tilt[layer]);
0826     }
0827     */
0828
0829     PHGenFit::Measurement* meas = new PHGenFit::PlanarMeasurement(pos, n,
0830         cluster->get_rphi_error(), cluster->get_z_error());
0831
0832     measurements.push_back(meas);
0833   }
0834
0835
0836   //TODO Add multiple TrackRep choices.
0837   genfit::AbsTrackRep* rep = new genfit::RKTrackRep(_primary_pid_guess);
0838   PHGenFit::Track* track(new PHGenFit::Track(rep, seed_pos, seed_mom, seed_cov));
0839   track->addMeasurements(measurements);
0840
0841   if (_fitter->processTrack(track, false) != 0) {
0842     if (verbosity >= 1)
0843       LogWarning("Track fitting failed");
0844     return NULL;
0845   }
0846
0847   TVector3 vertex_position(0, 0, 0);
0848   if (vertex) {
0849     vertex_position.SetXYZ(vertex->get_x(), vertex->get_y(), vertex->get_z());
0850   }
0851
0852   std::shared_ptr<genfit::MeasuredStateOnPlane> gf_state_vertex_ca = NULL;
0853   try {
0854     gf_state_vertex_ca = std::shared_ptr < genfit::MeasuredStateOnPlane> (track->extrapolateToPoint(vertex_position));
0855   } catch (...) {
0856     if (verbosity >= 2)
0857       LogWarning("extrapolateToPoint failed!");
0858     return NULL;
0859   }
0860
0861   TVector3 mom = gf_state_vertex_ca->getMom();
0862   TMatrixDSym cov = gf_state_vertex_ca->get6DCov();
0863
0864   //cout << "OUT Ex: " << sqrt(cov[0][0]) << ", Ey: " << sqrt(cov[1][1]) << ", Ez: " << sqrt(cov[2][2]) << std::endl;
0865
0866   //cout << "IN Px: " << intrack->get_px() << ", Py: " << intrack->get_py() << ", Pz: " << intrack->get_pz() << std::endl;
0867   //cout << "OUT Px: " << mom.X() << ", Py: " << mom.Y() << ", Pz: " << mom.Z() << std::endl;
0868
0869   return track;
0870 }
0871
0872 /*
0873  * Fill SvtxVertexMap from GFRaveVertexes and Tracks
0874  */
0875 void BDiJetModule::FillVertexMap(
0876   const std::vector<genfit::GFRaveVertex*>& rave_vertices,
0877   const std::vector<genfit::Track*>& gf_tracks) {
0878
0879   for (unsigned int ivtx = 0; ivtx < rave_vertices.size(); ++ivtx)
0880   {
0881     genfit::GFRaveVertex* rave_vtx = rave_vertices[ivtx];
0882
0883     //cout << "V0 x: " << rave_vtx->getPos().X() << ", y: " << rave_vtx->getPos().Y() << ", z: " << rave_vtx->getPos().Z() << std::endl;
0884     rv_prim_vtx[0] = rave_vtx->getPos().X();
0885     rv_prim_vtx[1] = rave_vtx->getPos().Y();
0886     rv_prim_vtx[2] = rave_vtx->getPos().Z();
0887
0888     rv_prim_vtx_err[0] = sqrt(rave_vtx->getCov()[0][0]);
0889     rv_prim_vtx_err[1] = sqrt(rave_vtx->getCov()[1][1]);
0890     rv_prim_vtx_err[2] = sqrt(rave_vtx->getCov()[2][2]);
0891
0892     rv_prim_vtx_ntrk = rave_vtx->getNTracks();
0893   } // ivtx
0894
0895   for (SvtxVertexMap::Iter iter = _vertexmap->begin();
0896        iter != _vertexmap->end();
0897        ++iter)
0898   {
0899     SvtxVertex *svtx_vertex = iter->second;
0900
0901     //cout << "V1 x: " << svtx_vertex->get_x() << ", y: " << svtx_vertex->get_y() << ", z: " << svtx_vertex->get_z() << std::endl;
0902     gf_prim_vtx[0] = svtx_vertex->get_x();
0903     gf_prim_vtx[1] = svtx_vertex->get_y();
0904     gf_prim_vtx[2] = svtx_vertex->get_z();
0905
0906     gf_prim_vtx_err[0] = sqrt(svtx_vertex->get_error(0, 0));
0907     gf_prim_vtx_err[1] = sqrt(svtx_vertex->get_error(1, 1));
0908     gf_prim_vtx_err[2] = sqrt(svtx_vertex->get_error(2, 2));
0909
0910     gf_prim_vtx_ntrk = int(svtx_vertex->size_tracks());
0911   } // iter on SvtxVertexMap
0912
0913   return;
0914 }
0915
0916 int BDiJetModule::GetSVMass_mom(
0917   const genfit::GFRaveVertex* rave_vtx,
0918   float & vtx_mass,
0919   float & vtx_px,
0920   float & vtx_py,
0921   float & vtx_pz
0922 ) {
0923
0924   float sum_E = 0, sum_px = 0, sum_py = 0, sum_pz = 0;
0925
0926   int N_good = 0;
0927
0928   for (unsigned int itrk = 0; itrk < rave_vtx->getNTracks(); itrk++) {
0929     TVector3 mom_trk = rave_vtx->getParameters(itrk)->getMom();
0930
0931     double w_trk = rave_vtx->getParameters(itrk)->getWeight();
0932
0933     sum_px += mom_trk.X();
0934     sum_py += mom_trk.Y();
0935     sum_pz += mom_trk.Z();
0936     sum_E += sqrt(mom_trk.Mag2() + 0.140 * 0.140);
0937
0938     //cout << "W: " << w_trk << std::endl;
0939     if ( w_trk > 0.6 ) {
0940       N_good++;
0941     }
0942
0943   }//itrk
0944
0945   vtx_mass =  sqrt(sum_E * sum_E - sum_px * sum_px - sum_py * sum_py - sum_pz * sum_pz);
0946   vtx_px = sum_px;
0947   vtx_py = sum_py;
0948   vtx_pz = sum_pz;
0949
0950   //cout << "Mass: " << vtx_mass << ", pT: " << vtx_pT << std::endl;
0951   return N_good;
0952 }
0953
0954
0955