Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:12:46 

0001 /*!
0002  *  \file       BJetModule.cc
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 "BJetModule.h"
0009 
0010 #include <hfjettruthgeneration/HFJetDefs.h>
0011 
0012 #include <Acts/Definitions/Algebra.hpp>
0013 
0014 #include <fun4all/Fun4AllReturnCodes.h>
0015 #include <fun4all/Fun4AllServer.h>
0016 
0017 #include <phool/PHCompositeNode.h>
0018 #include <phool/getClass.h>
0019 
0020 #include <g4main/PHG4Hit.h>
0021 #include <g4main/PHG4Particle.h>
0022 #include <g4main/PHG4TruthInfoContainer.h>
0023 #include <g4main/PHG4VtxPoint.h>
0024 
0025 #include <g4jets/Jet.h>
0026 #include <g4jets/JetMap.h>
0027 
0028 #include <trackbase_historic/SvtxTrack.h>
0029 #include <trackbase_historic/SvtxTrackMap.h>
0030 #include <trackbase_historic/SvtxVertex.h>
0031 #include <trackbase_historic/SvtxVertexMap.h>
0032 
0033 #include <g4eval/JetEvalStack.h>
0034 #include <g4eval/SvtxEvalStack.h>
0035 
0036 #include <HepMC/GenEvent.h>
0037 #include <HepMC/GenVertex.h>
0038 #include <phhepmc/PHHepMCGenEvent.h>
0039 #include <phhepmc/PHHepMCGenEventMap.h>
0040 
0041 #include "TDatabasePDG.h"
0042 #include "TLorentzVector.h"
0043 
0044 #include <iostream>
0045 #include <memory>
0046 
0047 //#define _DEBUG_
0048 
0049 #define LogDebug(exp) std::cout << "DEBUG: " << __FILE__ << ": " << __LINE__ << ": " << exp << "\n"
0050 #define LogError(exp) std::cout << "ERROR: " << __FILE__ << ": " << __LINE__ << ": " << exp << "\n"
0051 #define LogWarning(exp) std::cout << "WARNING: " << __FILE__ << ": " << __LINE__ << ": " << exp << "\n"
0052 
0053 using namespace std;
0054 
0055 BJetModule::BJetModule(const string &name, const string &out)
0056   : SubsysReco(name)
0057   , _foutname(out)
0058   , _truthjetmap_name("AntiKt_Truth_r04")
0059   , _recojetmap_name("AntiKt_Tower_r04")
0060   , _trackmap_name("SvtxTrackMap")
0061   , _vertexmap_name("SvtxVertexMap")
0062   , _embedding_id(1)
0063 {
0064 }
0065 
0066 int BJetModule::Init(PHCompositeNode *topNode)
0067 {
0068   _ievent = 0;
0069 
0070   _b_event = -1;
0071 
0072   reset_tree_vars();
0073 
0074   _f = new TFile(_foutname.c_str(), "RECREATE");
0075 
0076   _tree = new TTree("T", "a withered deciduous tree");
0077 
0078   _tree->Branch("event", &_b_event, "event/I");
0079   setBranches();
0080 
0081   return 0;
0082 }
0083 
0084 double calc_dca(const TVector3 &track_point, const TVector3 &track_direction, const TVector3 &vertex)
0085 {
0086   TVector3 VP = vertex - track_point;
0087   double d = -99;
0088 
0089   if (track_direction.Mag() > 0)
0090   {
0091     d = (track_direction.Cross(VP)).Mag() / track_direction.Mag();
0092   }
0093 
0094   return d;
0095 }
0096 
0097 bool calc_dca3d_line(
0098     double &dca_xy, double &dca_z, /*cm*/
0099     const TVector3 &track_point, const TVector3 &track_direction, const TVector3 &vertex)
0100 {
0101   dca_xy = NAN;
0102   dca_z = NAN;
0103 
0104   if (track_direction.Mag() == 0) return false;
0105 
0106   TVector3 PV = track_point - vertex;
0107   if (PV.Mag() < 0.000001)
0108   {
0109     dca_xy = 0;
0110     dca_z = 0;
0111     return true;
0112   }
0113 
0114   TVector3 PVxMom = track_direction.Cross(PV);
0115 
0116   TVector3 PCA = track_direction;
0117 
0118   if (PVxMom.Mag() < 0.000001)
0119   {
0120     std::cout << __FILE__ << ": " << __LINE__ << ": PVxMom.Mag2() < 0.000001" << std::endl;
0121     return false;
0122   }
0123 
0124   PCA.Rotate(TMath::PiOver2(), PVxMom);  // direction
0125 
0126   PCA.SetMag(1.);
0127 
0128   PCA.SetMag(PV.Dot(PCA));
0129 
0130   TVector3 r = track_direction.Cross(TVector3(0, 0, 1));
0131   r.SetMag(1.);
0132 
0133   dca_xy = PCA.Dot(r);
0134 
0135   dca_z = PCA.Z();
0136 
0137   return true;
0138 }
0139 
0140 bool calc_dca_circle_line(
0141     double &dca_xy, double &dca_z,                                                                   /*cm*/
0142     const TVector3 &track_point /*cm*/, const TVector3 &track_mom /*GeV/c*/, const TVector3 &vertex, /*cm*/
0143     const double &q = 1. /*e*/, const double &B = 1.4 /*Tesla*/)
0144 {
0145   double z0 = track_point.Z() - vertex.Z();
0146   double r0 = (track_point - vertex).Perp();
0147   double pz = track_mom.Z();
0148   double pt = track_mom.Perp();
0149 
0150   if (pt == 0 || B == 0 || q == 0)
0151   {
0152     LogWarning("pt == 0 || B ==0 || q==0");
0153     dca_z = NAN;
0154     dca_xy = NAN;
0155     return false;
0156   }
0157 
0158   dca_z = z0 - pz / pt * r0;
0159 
0160   const double ten_over_c = 333.564095198;  // GeV/c in cm*e*T
0161   double R = pt / B / q * ten_over_c;       // radius in cm
0162 
0163   // make 2d versions
0164   TVector3 track_point_2d(track_point.X(), track_point.Y(), 0);
0165   TVector3 track_mom_2d(track_mom.X(), track_mom.Y(), 0);
0166   TVector3 vertex_2d(vertex.X(), vertex.Y(), 0);
0167 
0168   // calc track circle center
0169   TVector3 track_2d_circle_center = track_mom_2d;    //copy mom_2d
0170   track_2d_circle_center.RotateZ(TMath::PiOver2());  //rotate Pi/2
0171   track_2d_circle_center.SetMag(R);                  // scale to R, R could be negtive
0172   track_2d_circle_center += track_point_2d;          //shift base to point_2d
0173 
0174   dca_xy = TMath::Abs(R) - (track_2d_circle_center - vertex_2d).Mag();
0175 
0176   return true;
0177 }
0178 
0179 int BJetModule::process_event(PHCompositeNode *topNode)
0180 {
0181   _b_event = _ievent;
0182 
0183   reset_tree_vars();
0184 
0185   PHHepMCGenEventMap *geneventmap = findNode::getClass<PHHepMCGenEventMap>(
0186       topNode, "PHHepMCGenEventMap");
0187   if (!geneventmap)
0188   {
0189     std::cout << PHWHERE
0190               << " - Fatal error - missing node PHHepMCGenEventMap"
0191               << std::endl;
0192     return Fun4AllReturnCodes::ABORTRUN;
0193   }
0194 
0195   PHHepMCGenEvent *genevt = geneventmap->get(_embedding_id);
0196   if (!genevt)
0197   {
0198     std::cout << PHWHERE
0199               << " - Fatal error - node PHHepMCGenEventMap missing subevent with embedding ID "
0200               << _embedding_id;
0201     std::cout << ". Print PHHepMCGenEventMap:";
0202     geneventmap->identify();
0203     return Fun4AllReturnCodes::ABORTRUN;
0204   }
0205 
0206   //PHHepMCGenEvent *genevt = findNode::getClass<PHHepMCGenEvent>(topNode,"PHHepMCGenEvent");
0207   HepMC::GenEvent *theEvent = genevt->getEvent();
0208   //theEvent->print();
0209 
0210   PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0211 
0212   PHG4VtxPoint *first_point = truthinfo->GetPrimaryVtx(truthinfo->GetPrimaryVertexIndex());
0213   _b_truth_vertex_n = 0;
0214   _b_truth_vertex_x[_b_truth_vertex_n] = first_point->get_x();
0215   _b_truth_vertex_y[_b_truth_vertex_n] = first_point->get_y();
0216   _b_truth_vertex_z[_b_truth_vertex_n] = first_point->get_z();
0217   TVector3 truth_primary_vertex(first_point->get_x(), first_point->get_y(), first_point->get_z());
0218   ++_b_truth_vertex_n;
0219 
0220   auto jet_eval_stack = unique_ptr<JetEvalStack>(new JetEvalStack(topNode, _recojetmap_name, _truthjetmap_name));
0221   if (!jet_eval_stack)
0222   {
0223     LogError("!jet_eval_stack");
0224     return Fun4AllReturnCodes::ABORTRUN;
0225   }
0226 
0227   JetRecoEval *jet_reco_eval = jet_eval_stack->get_reco_eval();
0228   if (!jet_reco_eval)
0229   {
0230     LogError("!jet_reco_eval");
0231     return Fun4AllReturnCodes::ABORTRUN;
0232   }
0233 
0234   JetMap *truth_jets = findNode::getClass<JetMap>(topNode, _truthjetmap_name);
0235   JetMap *reco_jets = findNode::getClass<JetMap>(topNode, _recojetmap_name);
0236  
0237   _b_truthjet_n = 0;
0238   for (JetMap::Iter iter = truth_jets->begin(); iter != truth_jets->end();
0239        ++iter)
0240   {
0241     Jet *truth_jet = iter->second;
0242     if (truth_jet->get_pt() < 10 || fabs(truth_jet->get_eta()) > 2)
0243       continue;
0244 
0245     _b_truthjet_pt[_b_truthjet_n] = truth_jet->get_pt();
0246     _b_truthjet_phi[_b_truthjet_n] = truth_jet->get_phi();
0247     _b_truthjet_eta[_b_truthjet_n] = truth_jet->get_eta();
0248 
0249     int jet_flavor = -999;
0250 
0251     jet_flavor = truth_jet->get_property(static_cast<Jet::PROPERTY>(prop_JetPartonFlavor));
0252     if (abs(jet_flavor) < 100)
0253       _b_truthjet_parton_flavor[_b_truthjet_n] = jet_flavor;
0254 
0255     jet_flavor = truth_jet->get_property(static_cast<Jet::PROPERTY>(prop_JetHadronFlavor));
0256     if (abs(jet_flavor) < 100)
0257       _b_truthjet_hadron_flavor[_b_truthjet_n] = jet_flavor;
0258     //cout << "DEBUG: " << __LINE__ << endl;
0259     //auto reco_jet = unique_ptr<Jet>(jet_reco_eval->best_jet_from(truth_jet));
0260     Jet *reco_jet = jet_reco_eval->best_jet_from(truth_jet);
0261     
0262     //cout << "DEBUG: " << __LINE__ << endl;
0263     if (!reco_jet)
0264     {
0265       //// match by radius
0266       Jet* matchedjet = nullptr;
0267       float mindr = std::numeric_limits<float>::max();
0268       for (JetMap::Iter riter = reco_jets->begin(); riter != reco_jets->end();
0269        ++riter)
0270     {
0271       Jet* mjet = riter->second;
0272       if(mjet->get_pt() < 5)
0273         { continue; }
0274       float dr = dR(mjet->get_eta(), truth_jet->get_eta(),
0275             mjet->get_phi(), truth_jet->get_phi());
0276       if(dr < mindr)
0277         {
0278           mindr = dr;
0279           matchedjet = mjet;
0280         }
0281     }
0282        
0283       if(matchedjet)
0284     {
0285       _b_recojet_valid[_b_truthjet_n] = 1;
0286       _b_recojet_pt[_b_truthjet_n] = matchedjet->get_pt();
0287       _b_recojet_phi[_b_truthjet_n] = matchedjet->get_phi();
0288       _b_recojet_eta[_b_truthjet_n] = matchedjet->get_eta();
0289     }
0290       else
0291     {
0292       _b_recojet_valid[_b_truthjet_n] = 0;
0293     }
0294     }
0295     else
0296     {
0297       _b_recojet_valid[_b_truthjet_n] = 1;
0298       _b_recojet_pt[_b_truthjet_n] = reco_jet->get_pt();
0299       _b_recojet_phi[_b_truthjet_n] = reco_jet->get_phi();
0300       _b_recojet_eta[_b_truthjet_n] = reco_jet->get_eta();
0301     }
0302 
0303     _b_truthjet_n++;
0304     //cout << "DEBUG: " << __LINE__ << endl;
0305   }
0306 
0307   PHG4TruthInfoContainer::Range range = truthinfo->GetPrimaryParticleRange();
0308   //PHG4TruthInfoContainer::Range range = truthinfo->GetParticleRange();
0309 
0310   _b_particle_n = 0;
0311   for (auto iter = range.first; iter != range.second; ++iter)
0312   {
0313     PHG4Particle *g4particle = iter->second;  // You may ask yourself, why second?
0314 
0315     unsigned int embed_id = truthinfo->isEmbeded(g4particle->get_track_id());
0316 
0317     TLorentzVector t;
0318     t.SetPxPyPzE(g4particle->get_px(), g4particle->get_py(), g4particle->get_pz(), g4particle->get_e());
0319 
0320     float truth_pt = t.Pt();
0321     if (truth_pt < 0.5) continue;
0322     float truth_eta = t.Eta();
0323     if (fabs(truth_eta) > 1.1) continue;
0324     float truth_phi = t.Phi();
0325     int truth_pid = g4particle->get_pid();
0326 
0327     //original cuts
0328     /*
0329         if (truth_pid == 22 || truth_pid == 2112 || truth_pid == -2112
0330                 || truth_pid == 130)
0331             continue;
0332         if (truth_pid == 2112 || truth_pid == -2112 || truth_pid == 130)
0333             continue;
0334         // save high-pT photons
0335         //if (truth_pid == 22 && truth_pt < 20) continue;
0336         if (truth_pid == 12 || truth_pid == -12 || truth_pid == 16
0337                 || truth_pid == -16 || truth_pid == 14 || truth_pid == -14)
0338             continue;
0339             */
0340 
0341     //only keeps stable charged particles
0342     if (!(abs(truth_pid) == 211 || abs(truth_pid) == 321 || abs(truth_pid) == 2212 || abs(truth_pid) == 11 || abs(truth_pid) == 13))
0343       continue;
0344 
0345     _b_particle_pid[_b_particle_n] = truth_pid;
0346     _b_particle_pt[_b_particle_n] = truth_pt;
0347     _b_particle_eta[_b_particle_n] = truth_eta;
0348     _b_particle_phi[_b_particle_n] = truth_phi;
0349     _b_particle_embed[_b_particle_n] = embed_id;
0350 
0351     PHG4VtxPoint *point = truthinfo->GetVtx(g4particle->get_vtx_id());
0352     _b_particle_vertex_x[_b_particle_n] = point->get_x();
0353     _b_particle_vertex_y[_b_particle_n] = point->get_y();
0354     _b_particle_vertex_z[_b_particle_n] = point->get_z();
0355 
0356     TVector3 track_point(point->get_x(), point->get_y(), point->get_z());
0357     TVector3 track_mom(g4particle->get_px(), g4particle->get_py(), g4particle->get_pz());
0358 
0359     double truth_dca_xy = NAN;
0360     double truth_dca_z = NAN;
0361 
0362     //double charge = TDatabasePDG::Instance()->GetParticle(truth_pid)->Charge()/3.;
0363     //calc_dca_circle_line(truth_dca_xy, truth_dca_z, track_point, track_mom, truth_primary_vertex, charge, 1.4);
0364 
0365     calc_dca3d_line(truth_dca_xy, truth_dca_z, track_point, track_mom, truth_primary_vertex);
0366 
0367 #ifdef _DEBUG_
0368     {
0369       cout << "track_point: --------------" << endl;
0370       track_point.Print();
0371 
0372       cout << "track_mom: --------------" << endl;
0373       track_mom.Print();
0374 
0375       cout << "truth_primary_vertex: --------------" << endl;
0376       truth_primary_vertex.Print();
0377 
0378       cout
0379           << __LINE__
0380           << ": " << TDatabasePDG::Instance()->GetParticle(truth_pid)->GetName()
0381           << ": truth_dca_xy: " << truth_dca_xy
0382           << ": truth_dca_z: " << truth_dca_z
0383           << endl
0384           << endl;
0385     }
0386 #endif
0387 
0388     _b_particle_dca_xy[_b_particle_n] = truth_dca_xy;
0389     _b_particle_dca_z[_b_particle_n] = truth_dca_z;
0390 
0391     //      for (HepMC::GenEvent::particle_const_iterator p =
0392     //              theEvent->particles_begin(); p != theEvent->particles_end();
0393     //              ++p) {
0394     //
0395     //          if ((*p)->pdg_id() != truth_pid)
0396     //              continue;
0397     //
0398     //          if (fabs(truth_pt - (*p)->momentum().perp()) > 0.001)
0399     //              continue;
0400     //          if (fabs(truth_eta - (*p)->momentum().pseudoRapidity()) > 0.001)
0401     //              continue;
0402     //          if (fabs(truth_phi - (*p)->momentum().phi()) > 0.001)
0403     //              continue;
0404     //
0405     //          HepMC::GenVertex *production_vertex = (*p)->production_vertex();
0406     //          _b_particle_dca_xy[_b_particle_n] = production_vertex->point3d().perp();
0407     //
0408     //          {
0409     //              cout
0410     //              << __LINE__
0411     //              <<": From HepMC: {"
0412     //              << production_vertex->point3d().x() <<", "
0413     //              << production_vertex->point3d().y() <<", "
0414     //              << production_vertex->point3d().z() <<"}"
0415     //              <<endl<<endl;
0416     //          }
0417     //      }
0418 
0419     _b_particle_n++;
0420   }
0421 
0422   SvtxTrackMap *trackmap = findNode::getClass<SvtxTrackMap>(topNode, _trackmap_name.c_str());
0423 
0424   SvtxVertexMap *vertexmap = findNode::getClass<SvtxVertexMap>(topNode, _vertexmap_name.c_str());
0425 
0426   auto svtxevalstack = unique_ptr<SvtxEvalStack>(new SvtxEvalStack(topNode));
0427   if (!svtxevalstack)
0428   {
0429     LogError("!svtxevalstack");
0430     return Fun4AllReturnCodes::ABORTRUN;
0431   }
0432 
0433   svtxevalstack->next_event(topNode);
0434   SvtxTrackEval *trackeval = svtxevalstack->get_track_eval();
0435 
0436   float vertex_x = -99;
0437   float vertex_y = -99;
0438   float vertex_z = -99;
0439   //    float vertex_err_x = -99;
0440   //    float vertex_err_y = -99;
0441   //float vertex_err_z = -99;
0442 
0443   _b_track_n = 0;
0444 
0445   for (SvtxTrackMap::Iter iter = trackmap->begin(); iter != trackmap->end();
0446        ++iter)
0447   {
0448     SvtxTrack *track = iter->second;
0449 
0450     float track_pt = track->get_pt();
0451     float track_eta = track->get_eta();
0452     float track_phi = track->get_phi();
0453 
0454     if (track_pt < 0.5)
0455       continue;
0456     if (fabs(track_eta) > 1.1)
0457       continue;
0458 
0459     //std::set<PHG4Hit*> assoc_hits = trackeval->all_truth_hits(track);//TODO
0460     int nmaps = 0;
0461 
0462     unsigned int nclusters = track->size_cluster_keys();
0463     unsigned int nclusters_by_layer = 0;
0464 
0465     TrackSeed *siliconSeed = track->get_silicon_seed();
0466     for (auto iter = siliconSeed->begin_cluster_keys();
0467          iter != siliconSeed->end_cluster_keys(); ++iter)
0468     {
0469       TrkrDefs::cluskey key = *iter;
0470       unsigned int cluster_layer = TrkrDefs::getLayer(key);
0471 
0472       if (TrkrDefs::getTrkrId(key) == TrkrDefs::TrkrId::mvtxId)
0473       {
0474         ++nmaps;
0475       }
0476       nclusters_by_layer |= (0x3FFFFFFF & (0x1 << cluster_layer));
0477     }
0478 
0479     TrackSeed *tpcSeed = track->get_silicon_seed();
0480     for (auto iter = tpcSeed->begin_cluster_keys();
0481          iter != tpcSeed->end_cluster_keys(); ++iter)
0482     {
0483       unsigned int cluster_layer = TrkrDefs::getLayer(*iter);
0484       nclusters_by_layer |= (0x3FFFFFFF & (0x1 << cluster_layer));
0485     }
0486 
0487     PHG4Particle *g4particle = trackeval->max_truth_particle_by_nclusters(
0488         track);
0489 
0490     unsigned int truth_nclusters = trackeval->get_nclusters_contribution(
0491         track, g4particle);
0492     unsigned int truth_nclusters_by_layer = trackeval->get_nclusters_contribution_by_layer(
0493         track, g4particle);
0494 
0495     unsigned int truth_embed_id = truthinfo->isEmbeded(
0496         g4particle->get_track_id());
0497     bool truth_is_primary = truthinfo->is_primary(g4particle);
0498     int truth_pid = g4particle->get_pid();
0499 
0500     //TVector3 g4particle_mom(g4particle->get_px(),g4particle->get_py(),g4particle->get_pz());
0501     //_b_track_best_pt[_b_track_n] = g4particle_mom.Pt();
0502 
0503     //int truth_parent_id = g4particle->get_parent_id();
0504     //int truth_primary_id = g4particle->get_primary_id();
0505 
0506     TVector3 track_point(track->get_x(), track->get_y(), track->get_z());
0507     TVector3 track_direction(track->get_px(), track->get_py(), track->get_pz());
0508     SvtxVertex* svertex = vertexmap->get(track->get_vertex_id());
0509     vertex_x = svertex->get_x();
0510     vertex_y = svertex->get_y();
0511     vertex_z = svertex->get_z();
0512     TVector3 vertex(vertex_x, vertex_y, vertex_z);
0513 
0514     TVector3 track_point_2d(track->get_x(), track->get_y(), 0);
0515     TVector3 track_direction_2d(track->get_px(), track->get_py(), 0);
0516     TVector3 vertex_2d(vertex_x, vertex_y, 0);
0517     float dca3d_xy = NAN;
0518     float dca3d_xy_error = NAN;
0519     float dca3d_z = NAN;
0520     float dca3d_z_error = NAN;
0521     calc3DDCA(track, vertexmap, dca3d_xy, dca3d_xy_error, dca3d_z, dca3d_z_error);
0522 
0523     float dca2d_calc = calc_dca(track_point_2d, track_direction_2d, vertex_2d);
0524     float dca2d_calc_truth = calc_dca(track_point_2d, track_direction_2d, TVector3(0, 0, 0));
0525 
0526     float dca3d_calc = calc_dca(track_point, track_direction, vertex);
0527     float dca3d_calc_truth = calc_dca(track_point, track_direction, TVector3(0, 0, 0));
0528 
0529     //      float dca2d_calc_error = sqrt(
0530     //              dca2d_error*dca2d_error +
0531     //              vertex_err_x*vertex_err_x +
0532     //              vertex_err_y*vertex_err_y);
0533 
0534     _b_track_pca_x[_b_track_n] = track->get_x() - vertex_x;
0535     _b_track_pca_y[_b_track_n] = track->get_y() - vertex_y;
0536     _b_track_pca_z[_b_track_n] = track->get_z() - vertex_z;
0537 
0538     TVector3 dca_vector(
0539         _b_track_pca_x[_b_track_n],
0540         _b_track_pca_y[_b_track_n],
0541         _b_track_pca_z[_b_track_n]);
0542 
0543     _b_track_pca_phi[_b_track_n] = dca_vector.Phi();
0544 
0545     TLorentzVector t;
0546     t.SetPxPyPzE(g4particle->get_px(), g4particle->get_py(),
0547                  g4particle->get_pz(), g4particle->get_e());
0548 
0549     float truth_pt = t.Pt();
0550     float truth_eta = t.Eta();
0551     float truth_phi = t.Phi();
0552 
0553     PHG4VtxPoint *point = truthinfo->GetVtx(g4particle->get_vtx_id());
0554     _b_track_best_vertex_x[_b_track_n] = point->get_x();
0555     _b_track_best_vertex_y[_b_track_n] = point->get_y();
0556     _b_track_best_vertex_z[_b_track_n] = point->get_z();
0557 
0558     TVector3 track_best_point(point->get_x(), point->get_y(), point->get_z());
0559     TVector3 track_best_mom(g4particle->get_px(), g4particle->get_py(), g4particle->get_pz());
0560 
0561     double truth_dca_xy = NAN;
0562     double truth_dca_z = NAN;
0563 
0564     //double charge = TDatabasePDG::Instance()->GetParticle(truth_pid)->Charge()/3.;
0565     //calc_dca_circle_line(truth_dca_xy, truth_dca_z, track_best_point, track_best_mom, truth_primary_vertex, charge, 1.4);
0566 
0567     calc_dca3d_line(truth_dca_xy, truth_dca_z, track_best_point, track_best_mom, truth_primary_vertex);
0568 
0569     _b_track_best_dca_xy[_b_track_n] = truth_dca_xy;
0570     _b_track_best_dca_z[_b_track_n] = truth_dca_z;
0571 
0572     //_b_track_best_parent_pid[_b_track_n] = truthinfo->GetParticle(g4particle->get_parent_id())->get_pid();
0573     int truth_in = -1;
0574     int truth_out = -1;
0575     int nhepmc = 0;
0576     for (HepMC::GenEvent::particle_const_iterator p =
0577              theEvent->particles_begin();
0578          p != theEvent->particles_end();
0579          ++p)
0580     {
0581       if ((*p)->pdg_id() != truth_pid)
0582         continue;
0583 
0584       if (fabs(truth_pt - (*p)->momentum().perp()) > 0.001)
0585         continue;
0586       if (fabs(truth_eta - (*p)->momentum().pseudoRapidity()) > 0.001)
0587         continue;
0588       if (fabs(truth_phi - (*p)->momentum().phi()) > 0.001)
0589         continue;
0590 
0591       HepMC::GenVertex *production_vertex = (*p)->production_vertex();
0592 
0593       truth_in = production_vertex->particles_in_size();
0594       truth_out = production_vertex->particles_out_size();
0595 
0596       HepMC::GenVertex::particles_in_const_iterator first_parent =
0597           production_vertex->particles_in_const_begin();
0598       _b_track_best_parent_pid[_b_track_n] = (*first_parent)->pdg_id();
0599 
0600       nhepmc++;
0601     }
0602 
0603     _b_track_pt[_b_track_n] = track_pt;
0604     _b_track_eta[_b_track_n] = track_eta;
0605     _b_track_phi[_b_track_n] = track_phi;
0606 
0607     _b_track_dca3d_xy[_b_track_n] = dca3d_xy;
0608     _b_track_dca3d_xy_error[_b_track_n] = dca3d_xy_error;
0609 
0610     _b_track_dca3d_z[_b_track_n] = dca3d_z;
0611     _b_track_dca3d_z_error[_b_track_n] = dca3d_z_error;
0612 
0613     _b_track_dca2d_calc[_b_track_n] = dca2d_calc;
0614     _b_track_dca2d_calc_truth[_b_track_n] = dca2d_calc_truth;
0615 
0616     _b_track_dca3d_calc[_b_track_n] = dca3d_calc;
0617     _b_track_dca3d_calc_truth[_b_track_n] = dca3d_calc_truth;
0618 
0619     //      _b_track_pca_phi[_b_track_n] = dca_phi;
0620     //      _b_track_pca_x[_b_track_n] = dca_x;
0621     //      _b_track_pca_y[_b_track_n] = dca_y;
0622     //      _b_track_pca_z[_b_track_n] = dca_z;
0623 
0624     _b_track_quality[_b_track_n] = track->get_quality();
0625     _b_track_chisq[_b_track_n] = track->get_chisq();
0626     _b_track_ndf[_b_track_n] = track->get_ndf();
0627 
0628     _b_track_nmaps[_b_track_n] = nmaps;
0629 
0630     _b_track_nclusters[_b_track_n] = nclusters;
0631     _b_track_nclusters_by_layer[_b_track_n] = nclusters_by_layer;
0632 
0633     _b_track_best_nclusters[_b_track_n] = truth_nclusters;
0634     _b_track_best_nclusters_by_layer[_b_track_n] = truth_nclusters_by_layer;
0635     _b_track_best_embed[_b_track_n] = truth_embed_id;
0636     _b_track_best_primary[_b_track_n] = truth_is_primary;
0637     _b_track_best_pid[_b_track_n] = truth_pid;
0638     _b_track_best_pt[_b_track_n] = truth_pt;
0639 
0640     _b_track_best_in[_b_track_n] = truth_in;
0641     _b_track_best_out[_b_track_n] = truth_out;
0642 
0643     _b_track_n++;
0644   }
0645 
0646   _tree->Fill();
0647 
0648   _ievent++;
0649 
0650   //delete jet_eval_stack;
0651   //delete svtxevalstack;
0652 
0653   return 0;
0654 }
0655 
0656 void BJetModule::calc3DDCA(SvtxTrack* track, SvtxVertexMap* vertexmap, 
0657                float& dca3d_xy, float& dca3d_xy_error,
0658                float& dca3d_z, float& dca3d_z_error)
0659 {
0660   Acts::Vector3 pos(track->get_x(),
0661             track->get_y(),
0662             track->get_z());
0663   Acts::Vector3 mom(track->get_px(),
0664             track->get_py(),
0665             track->get_pz());
0666 
0667   auto vtxid = track->get_vertex_id();
0668   auto svtxVertex = vertexmap->get(vtxid);
0669   Acts::Vector3 vertex(svtxVertex->get_x(),
0670                svtxVertex->get_y(),
0671                svtxVertex->get_z());
0672 
0673   pos -= vertex;
0674 
0675   Acts::ActsSymMatrix<3> posCov;
0676   for(int i = 0; i < 3; ++i)
0677     {
0678       for(int j = 0; j < 3; ++j)
0679     {
0680       posCov(i, j) = track->get_error(i, j);
0681     } 
0682     }
0683   
0684   Acts::Vector3 r = mom.cross(Acts::Vector3(0.,0.,1.));
0685   float phi = atan2(r(1), r(0));
0686   
0687   Acts::RotationMatrix3 rot;
0688   Acts::RotationMatrix3 rot_T;
0689   rot(0,0) = cos(phi);
0690   rot(0,1) = -sin(phi);
0691   rot(0,2) = 0;
0692   rot(1,0) = sin(phi);
0693   rot(1,1) = cos(phi);
0694   rot(1,2) = 0;
0695   rot(2,0) = 0;
0696   rot(2,1) = 0;
0697   rot(2,2) = 1;
0698   
0699   rot_T = rot.transpose();
0700 
0701   Acts::Vector3 pos_R = rot * pos;
0702   Acts::ActsSymMatrix<3> rotCov = rot * posCov * rot_T;
0703 
0704   dca3d_xy= pos_R(0);
0705   dca3d_z = pos_R(2);
0706   dca3d_xy_error = sqrt(rotCov(0,0));
0707   dca3d_z_error = sqrt(rotCov(2,2));
0708 }
0709 
0710 int BJetModule::reset_tree_vars()
0711 {
0712   _b_truthjet_n = 0;
0713 
0714   for (int n = 0; n < 10; n++)
0715   {
0716     _b_truthjet_parton_flavor[n] = -9999;
0717     _b_truthjet_hadron_flavor[n] = -9999;
0718 
0719     _b_truthjet_pt[n] = -99;
0720     _b_truthjet_eta[n] = -99;
0721     _b_truthjet_phi[n] = -99;
0722 
0723     _b_recojet_valid[n] = 0;
0724     _b_truthjet_pt[n] = -99;
0725     _b_truthjet_eta[n] = -99;
0726     _b_truthjet_phi[n] = -99;
0727   }
0728 
0729   _b_particle_n = 0;
0730   _b_track_n = 0;
0731 
0732   for (int n = 0; n < 1000; n++)
0733   {
0734     _b_particle_pt[n] = -1;
0735     _b_particle_eta[n] = -1;
0736     _b_particle_phi[n] = -1;
0737     _b_particle_pid[n] = 0;
0738     _b_particle_embed[n] = 0;
0739 
0740     _b_particle_vertex_x[n] = -1;
0741     _b_particle_vertex_y[n] = -1;
0742     _b_particle_vertex_z[n] = -1;
0743 
0744     _b_particle_dca_xy[n] = -1;
0745     _b_particle_dca_z[n] = -1;
0746 
0747     _b_track_pt[n] = -1;
0748     _b_track_eta[n] = -1;
0749     _b_track_phi[n] = -1;
0750 
0751     _b_track_nmaps[n] = 0;
0752 
0753     _b_track_nclusters[n] = 0;
0754     _b_track_nclusters_by_layer[n] = 0;
0755 
0756     _b_track_dca2d[n] = -1;
0757     _b_track_dca2d_error[n] = -1;
0758     _b_track_dca3d_xy[n] = -1;
0759     _b_track_dca3d_xy_error[n] = -1;
0760     _b_track_dca3d_z[n] = -1;
0761     _b_track_dca3d_z_error[n] = -1;
0762     _b_track_dca2d_calc[n] = -1;
0763     _b_track_dca2d_calc_truth[n] = -1;
0764     _b_track_dca3d_calc[n] = -1;
0765     _b_track_dca3d_calc_truth[n] = -1;
0766 
0767     _b_track_pca_phi[n] = -99;
0768     _b_track_pca_x[n] = -99;
0769     _b_track_pca_y[n] = -99;
0770     _b_track_pca_z[n] = -99;
0771 
0772     _b_track_quality[n] = -99;
0773     _b_track_chisq[n] = -99;
0774     _b_track_ndf[n] = -99;
0775 
0776     _b_track_best_nclusters[n] = 0;
0777     _b_track_best_nclusters_by_layer[n] = 0;
0778     _b_track_best_embed[n] = 0;
0779     _b_track_best_primary[n] = false;
0780     _b_track_best_pid[n] = 0;
0781 
0782     _b_track_best_in[n] = 0;
0783     _b_track_best_out[n] = 0;
0784     _b_track_best_parent_pid[n] = 0;
0785 
0786     _b_track_best_vertex_x[n] = -99;
0787     _b_track_best_vertex_y[n] = -99;
0788     _b_track_best_vertex_z[n] = -99;
0789     _b_track_best_dca_xy[n] = -99;
0790     _b_track_best_dca_z[n] = -99;
0791   }
0792 
0793   return 0;
0794 }
0795 
0796 
0797 void BJetModule::setBranches()
0798 {
0799   
0800   _tree->Branch("truth_vertex_n", &_b_truth_vertex_n, "truth_vertex_n/I");
0801   _tree->Branch("truth_vertex_x", _b_truth_vertex_x, "truth_vertex_x[truth_vertex_n]/F");
0802   _tree->Branch("truth_vertex_y", _b_truth_vertex_y, "truth_vertex_y[truth_vertex_n]/F");
0803   _tree->Branch("truth_vertex_z", _b_truth_vertex_z, "truth_vertex_z[truth_vertex_n]/F");
0804 
0805   _tree->Branch("truthjet_n", &_b_truthjet_n, "truthjet_n/I");
0806   _tree->Branch("truthjet_parton_flavor", _b_truthjet_parton_flavor, "truthjet_parton_flavor[truthjet_n]/I");
0807   _tree->Branch("truthjet_hadron_flavor", _b_truthjet_hadron_flavor, "truthjet_hadron_flavor[truthjet_n]/I");
0808   _tree->Branch("truthjet_pt", _b_truthjet_pt, "truthjet_pt[truthjet_n]/F");
0809   _tree->Branch("truthjet_eta", _b_truthjet_eta, "truthjet_eta[truthjet_n]/F");
0810   _tree->Branch("truthjet_phi", _b_truthjet_phi, "truthjet_phi[truthjet_n]/F");
0811 
0812   _tree->Branch("recojet_valid", _b_recojet_valid, "recojet_valid[truthjet_n]/I");
0813   _tree->Branch("recojet_pt", _b_recojet_pt, "recojet_pt[truthjet_n]/F");
0814   _tree->Branch("recojet_eta", _b_recojet_eta, "recojet_eta[truthjet_n]/F");
0815   _tree->Branch("recojet_phi", _b_recojet_phi, "recojet_phi[truthjet_n]/F");
0816 
0817   _tree->Branch("particle_n", &_b_particle_n, "particle_n/I");
0818   _tree->Branch("particle_pt", _b_particle_pt, "particle_pt[particle_n]/F");
0819   _tree->Branch("particle_eta", _b_particle_eta, "particle_eta[particle_n]/F");
0820   _tree->Branch("particle_phi", _b_particle_phi, "particle_phi[particle_n]/F");
0821   _tree->Branch("particle_pid", _b_particle_pid, "particle_pid[particle_n]/I");
0822   _tree->Branch("particle_embed", _b_particle_embed, "particle_embed[particle_n]/i");
0823 
0824   _tree->Branch("particle_vertex_x", _b_particle_vertex_x, "particle_vertex_x[particle_n]/F");
0825   _tree->Branch("particle_vertex_y", _b_particle_vertex_y, "particle_vertex_y[particle_n]/F");
0826   _tree->Branch("particle_vertex_z", _b_particle_vertex_z, "particle_vertex_z[particle_n]/F");
0827   _tree->Branch("particle_dca_xy", _b_particle_dca_xy, "particle_dca_xy[particle_n]/F");
0828   _tree->Branch("particle_dca_z", _b_particle_dca_z, "particle_dca_z[particle_n]/F");
0829 
0830   _tree->Branch("track_n", &_b_track_n, "track_n/I");
0831   _tree->Branch("track_pt", _b_track_pt, "track_pt[track_n]/F");
0832   _tree->Branch("track_eta", _b_track_eta, "track_eta[track_n]/F");
0833   _tree->Branch("track_phi", _b_track_phi, "track_phi[track_n]/F");
0834 
0835   _tree->Branch("track_dca2d", _b_track_dca2d, "track_dca2d[track_n]/F");
0836   _tree->Branch("track_dca2d_error", _b_track_dca2d_error, "track_dca2d_error[track_n]/F");
0837 
0838   _tree->Branch("track_dca3d_xy", _b_track_dca3d_xy, "track_dca3d_xy[track_n]/F");
0839   _tree->Branch("track_dca3d_xy_error", _b_track_dca3d_xy_error, "track_dca3d_xy_error[track_n]/F");
0840 
0841   _tree->Branch("track_dca3d_z", _b_track_dca3d_z, "track_dca3d_z[track_n]/F");
0842   _tree->Branch("track_dca3d_z_error", _b_track_dca3d_z_error, "track_dca3d_z_error[track_n]/F");
0843 
0844   _tree->Branch("track_dca2d_calc", _b_track_dca2d_calc, "track_dca2d_calc[track_n]/F");
0845   _tree->Branch("track_dca2d_calc_truth", _b_track_dca2d_calc_truth, "track_dca2d_calc_truth[track_n]/F");
0846 
0847   _tree->Branch("track_dca3d_calc", _b_track_dca3d_calc, "track_dca3d_calc[track_n]/F");
0848   _tree->Branch("track_dca3d_calc_truth", _b_track_dca3d_calc_truth, "track_dca3d_calc_truth[track_n]/F");
0849 
0850   _tree->Branch("track_pca_phi", _b_track_pca_phi, "track_pca_phi[track_n]/F");
0851   _tree->Branch("track_pca_x", _b_track_pca_x, "track_pca_x[track_n]/F");
0852   _tree->Branch("track_pca_y", _b_track_pca_y, "track_pca_y[track_n]/F");
0853   _tree->Branch("track_pca_z", _b_track_pca_z, "track_pca_z[track_n]/F");
0854 
0855   _tree->Branch("track_quality", _b_track_quality, "track_quality[track_n]/F");
0856   _tree->Branch("track_chisq", _b_track_chisq, "track_chisq[track_n]/F");
0857   _tree->Branch("track_ndf", _b_track_ndf, "track_ndf[track_n]/I");
0858 
0859   _tree->Branch("track_nmaps", _b_track_nmaps, "track_nmaps[track_n]/I");
0860 
0861   _tree->Branch("track_nclusters", _b_track_nclusters, "track_nclusters[track_n]/i");
0862   _tree->Branch("track_nclusters_by_layer", _b_track_nclusters_by_layer, "track_nclusters_by_layer[track_n]/i");
0863 
0864   _tree->Branch("track_best_nclusters", _b_track_best_nclusters, "track_best_nclusters[track_n]/i");
0865   _tree->Branch("track_best_nclusters_by_layer", _b_track_best_nclusters_by_layer, "track_best_nclusters_by_layer[track_n]/i");
0866 
0867   _tree->Branch("track_best_embed", _b_track_best_embed, "track_best_embed[track_n]/i");
0868   _tree->Branch("track_best_primary", _b_track_best_primary, "track_best_primary[track_n]/O");
0869   _tree->Branch("track_best_pid", _b_track_best_pid, "track_best_pid[track_n]/I");
0870   _tree->Branch("track_best_pt", _b_track_best_pt, "track_best_pt[track_n]/F");
0871   _tree->Branch("track_best_in", _b_track_best_in, "track_best_in[track_n]/I");
0872   _tree->Branch("track_best_out", _b_track_best_out, "track_best_out[track_n]/I");
0873   _tree->Branch("track_best_parent_pid", _b_track_best_parent_pid, "track_best_parent_pid[track_n]/I");
0874 
0875   _tree->Branch("track_best_vertex_x", _b_track_best_vertex_x, "track_best_vertex_x[track_n]/F");
0876   _tree->Branch("track_best_vertex_y", _b_track_best_vertex_y, "track_best_vertex_y[track_n]/F");
0877   _tree->Branch("track_best_vertex_z", _b_track_best_vertex_z, "track_best_vertex_z[track_n]/F");
0878 
0879   _tree->Branch("track_best_dca_xy", _b_track_best_dca_xy, "track_best_dca_xy[track_n]/F");
0880   _tree->Branch("track_best_dca_z", _b_track_best_dca_z, "track_best_dca_z[track_n]/F");
0881 
0882 }
0883 int BJetModule::End(PHCompositeNode *topNode)
0884 {
0885   //_f->ls();
0886   //_tree->Write();
0887   _f->Write();
0888   _f->Close();
0889 
0890   return 0;
0891 }
[ 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!