Back to home page

sPhenix code displayed by LXR
 
 

    

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

0001 #include "TagAndProbe.h"
0002 
0003 #include <globalvertex/SvtxVertexMap.h>
0004 #include <trackbase_historic/SvtxTrackMap.h>
0005 
0006 #include <fun4all/Fun4AllReturnCodes.h>
0007 
0008 #include <phool/getClass.h>
0009 
0010 #include <TEntryList.h>
0011 #include <TFile.h>
0012 #include <TLeaf.h>
0013 #include <TSystem.h>
0014 #include <TTree.h>  // for getting the TTree from the file
0015 #include <TVector3.h>
0016 #include <TLorentzVector.h>
0017 
0018 //#include <KFPVertex.h>
0019 //#include <KFParticle.h>           // for KFParticle
0020 #include <fun4all/Fun4AllBase.h>  // for Fun4AllBase::VERBOSITY...
0021 #include <fun4all/SubsysReco.h>   // for SubsysReco
0022 
0023 #include <ffamodules/CDBInterface.h>  // for accessing the field map file from the CDB
0024 #include <cctype>                     // for toupper
0025 #include <cmath>                      // for sqrt
0026 #include <cstdlib>                    // for size_t, exit
0027 #include <filesystem>
0028 #include <iostream>  // for operator<<, endl, basi...
0029 #include <map>       // for map
0030 #include <tuple>     // for tie, tuple
0031 
0032 class PHCompositeNode;
0033 
0034 // Tag and Probe Constructor
0035 TagAndProbe::TagAndProbe()
0036   : SubsysReco("TAGANDPROBE")
0037   , m_save_output(true)
0038   , m_outfile_name("outputTAP.root")
0039   , m_outfile(nullptr)
0040   , m_run(0)
0041   , m_segment(0)
0042 {
0043 }
0044 
0045 TagAndProbe::TagAndProbe(const std::string &name, const float run, const float seg)
0046   : SubsysReco(name)
0047   , m_save_output(true)
0048   , m_outfile_name("outputTAP.root")
0049   , m_outfile(nullptr)
0050   , m_run(run)
0051   , m_segment(seg)
0052 {
0053 }
0054 
0055 int TagAndProbe::Init(PHCompositeNode *topNode)
0056 {
0057   assert(topNode);
0058   if (m_save_output && Verbosity() >= VERBOSITY_SOME)
0059   {
0060     std::cout << "Output nTuple: " << m_outfile_name << std::endl;
0061   }
0062   if (m_save_output)
0063   {
0064     m_outfile = TFile::Open(m_outfile_name.c_str(),"RECREATE");
0065   }
0066 
0067   initializeBranches();
0068 
0069   m_event = 0;
0070 
0071   return 0;
0072 }
0073 
0074 int TagAndProbe::InitRun(PHCompositeNode *topNode)
0075 {
0076   assert(topNode);
0077 
0078   m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0079   if (!m_tGeometry)
0080   {
0081     std::cout << PHWHERE << "Error, can't find acts tracking geometry" << std::endl;
0082     return Fun4AllReturnCodes::ABORTEVENT;
0083   }
0084 
0085   //getField();
0086 
0087   m_cutInfoTree->Fill();
0088  
0089   return 0;
0090 }
0091 
0092 int TagAndProbe::process_event(PHCompositeNode *topNode)
0093 { 
0094   SvtxTrackMap *m_trackmap = findNode::getClass<SvtxTrackMap>(topNode, m_trk_map_node_name);
0095   int multiplicity = m_trackmap->size();
0096   if (multiplicity == 0)
0097   {
0098     if (Verbosity() >= VERBOSITY_SOME)
0099     {
0100       std::cout << "TagAndProbe: Event skipped as there are no tracks" << std::endl;
0101     }
0102     return Fun4AllReturnCodes::ABORTEVENT;
0103   }
0104 
0105   m_vertexMap = findNode::getClass<SvtxVertexMap>(topNode, m_vtx_map_node_name);
0106   if (m_vertexMap->size() == 0)
0107   {
0108     if (Verbosity() >= VERBOSITY_SOME)
0109     {
0110       std::cout << "TagAndProbe: Event skipped as there are no vertices" << std::endl;
0111     }
0112     return Fun4AllReturnCodes::ABORTEVENT;
0113   }
0114 
0115   int t1 = 0;
0116   for (const auto &[key1, track1] : *m_trackmap)
0117   {
0118     if (!track1)
0119     {
0120       continue;
0121     }
0122 
0123     m_trackID_1 = track1->get_id(); 
0124     m_crossing_1 = track1->get_crossing();   
0125     m_px_1 = track1->get_px();
0126     m_py_1 = track1->get_py();
0127     m_pz_1 = track1->get_pz();
0128     m_pt_1 = track1->get_pt();
0129     m_eta_1 = track1->get_eta();
0130     m_phi_1 = track1->get_phi();
0131     m_charge_1 = track1->get_charge();
0132     m_quality_1 = track1->get_quality();
0133     m_chisq_1 = track1->get_chisq();
0134     m_ndf_1 = track1->get_ndf();
0135 
0136     TrackSeed* silseed1 = track1->get_silicon_seed();
0137     TrackSeed* tpcseed1 = track1->get_tpc_seed();   
0138     m_nmaps_1 = 0;
0139     m_nintt_1 = 0;
0140     m_ntpc_1 = 0;
0141     m_nmms_1 = 0;
0142     m_nhits_1 = 0; 
0143     // +1 just in case _nlayers is zero
0144     std::vector<int> maps(_nlayers_maps+1, 0);
0145     std::vector<int> intt(_nlayers_intt+1, 0);
0146     std::vector<int> tpc(_nlayers_tpc+1, 0);
0147     std::vector<int> mms(_nlayers_mms+1, 0);
0148 
0149     if (tpcseed1)
0150     {
0151       m_nhits_1 += tpcseed1->size_cluster_keys();
0152       for (TrackSeed::ConstClusterKeyIter local_iter = tpcseed1->begin_cluster_keys();
0153            local_iter != tpcseed1->end_cluster_keys();
0154            ++local_iter)
0155       {
0156         TrkrDefs::cluskey cluster_key = *local_iter;
0157         //TrkrCluster* cluster = clustermap->findCluster(cluster_key);
0158         unsigned int local_layer = TrkrDefs::getLayer(cluster_key);
0159         if (_nlayers_maps > 0 && local_layer < _nlayers_maps)
0160         {
0161           maps[local_layer] = 1;
0162           m_nmaps_1++;
0163         }
0164         if (_nlayers_intt > 0 && local_layer >= _nlayers_maps && local_layer < _nlayers_maps + _nlayers_intt)
0165         {
0166           intt[local_layer - _nlayers_maps] = 1;
0167           m_nintt_1++;
0168         }
0169         if (_nlayers_tpc > 0 &&
0170             local_layer >= (_nlayers_maps + _nlayers_intt) &&
0171             local_layer < (_nlayers_maps + _nlayers_intt + _nlayers_tpc))
0172         {
0173           tpc[local_layer - (_nlayers_maps + _nlayers_intt)] = 1;
0174           m_ntpc_1++; 
0175         }
0176         if (_nlayers_mms > 0 && local_layer >= _nlayers_maps + _nlayers_intt + _nlayers_tpc)
0177         {
0178           mms[local_layer - (_nlayers_maps + _nlayers_intt + _nlayers_tpc)] = 1;
0179           m_nmms_1++;
0180         }
0181       }
0182     }
0183     if (silseed1)
0184     {
0185       m_nhits_1 += silseed1->size_cluster_keys();
0186       for (TrackSeed::ConstClusterKeyIter local_iter = silseed1->begin_cluster_keys();
0187            local_iter != silseed1->end_cluster_keys();
0188            ++local_iter)
0189       {
0190         TrkrDefs::cluskey cluster_key = *local_iter;
0191         //TrkrCluster* cluster = clustermap->findCluster(cluster_key);
0192         unsigned int local_layer = TrkrDefs::getLayer(cluster_key);
0193         if (_nlayers_maps > 0 && local_layer < _nlayers_maps)
0194         {
0195           maps[local_layer] = 1;
0196           m_nmaps_1++;
0197         }
0198         if (_nlayers_intt > 0 && local_layer >= _nlayers_maps && local_layer < _nlayers_maps + _nlayers_intt)
0199         {
0200           intt[local_layer - _nlayers_maps] = 1;
0201           m_nintt_1++;
0202         }
0203         if (_nlayers_tpc > 0 &&
0204             local_layer >= (_nlayers_maps + _nlayers_intt) &&
0205             local_layer < (_nlayers_maps + _nlayers_intt + _nlayers_tpc))
0206         {
0207           tpc[local_layer - (_nlayers_maps + _nlayers_intt)] = 1;
0208           m_ntpc_1++; 
0209         }
0210         if (_nlayers_mms > 0 && local_layer >= _nlayers_maps + _nlayers_intt + _nlayers_tpc)
0211         {
0212           mms[local_layer - (_nlayers_maps + _nlayers_intt + _nlayers_tpc)] = 1;
0213           m_nmms_1++;
0214         }
0215       }
0216     }
0217     float nlintt = 0;
0218     float nlmaps = 0;
0219     float nltpc = 0;
0220     float nlmms = 0;
0221     if (_nlayers_maps > 0)
0222     {
0223       for (unsigned int i = 0; i < _nlayers_maps; i++)
0224       {
0225         nlmaps += maps[i];
0226       }
0227     }  
0228     if (_nlayers_intt > 0)
0229     {
0230       for (unsigned int i = 0; i < _nlayers_intt; i++)
0231       {
0232         nlintt += intt[i];
0233       }
0234     }
0235     if (_nlayers_tpc > 0)
0236     {
0237       for (unsigned int i = 0; i < _nlayers_tpc; i++)
0238       {
0239         nltpc += tpc[i];
0240       }
0241     }
0242     if (_nlayers_mms > 0)
0243     {
0244       for (unsigned int i = 0; i < _nlayers_mms; i++)
0245       {
0246         nlmms += mms[i];
0247       }
0248     }
0249     m_nlayers_1 = nlmaps + nlintt + nltpc + nlmms;
0250 
0251     m_mapsstates_1 = 0;
0252     m_inttstates_1 = 0;
0253     m_tpcstates_1 = 0;
0254     m_mmsstates_1 = 0;
0255     for (auto state_iter = track1->begin_states();
0256        state_iter != track1->end_states();
0257        ++state_iter)
0258     {
0259       SvtxTrackState* tstate = state_iter->second;
0260       if (tstate->get_pathlength() != 0) //The first track state is an extrapolation so has no cluster
0261       {
0262         auto stateckey = tstate->get_cluskey();
0263         uint8_t id = TrkrDefs::getTrkrId(stateckey);
0264     
0265         switch (id)
0266         { 
0267           case TrkrDefs::mvtxId:
0268             ++m_mapsstates_1;
0269             break;
0270           case TrkrDefs::inttId:
0271             ++m_inttstates_1;
0272             break;
0273           case TrkrDefs::tpcId:
0274             ++m_tpcstates_1;
0275             break;
0276           case TrkrDefs::micromegasId:
0277             ++m_mmsstates_1;
0278             break;
0279           default:
0280            std::cout << "Cluster key doesnt match a tracking system, this shouldn't happen" << std::endl;
0281            break; 
0282         }
0283       }
0284     }
0285 
0286     if (m_include_pv_info)
0287     {
0288       // this is the global vertex id
0289       int _vertexID = track1->get_vertex_id();
0290       auto _vertex = m_vertexMap->get(_vertexID);
0291       if (_vertex)
0292       {
0293         m_vx = _vertex->get_x();
0294         m_vy = _vertex->get_y();
0295         m_vz = _vertex->get_z();
0296       }
0297       else 
0298       {
0299         m_vx = NAN;
0300         m_vy = NAN;
0301         m_vz = NAN;
0302       }
0303     }
0304   
0305     if (m_ntpc_1 >= m_nTPC_tag && m_nmaps_1 >= m_nMVTX_tag && m_nintt_1 >= m_nINTT_tag && 
0306         m_nmms_1 >= m_nTPOT_tag && m_quality_1 <= m_quality_tag)
0307     {
0308       m_tagpass_1 = true; 
0309     }
0310     else
0311     {
0312       m_tagpass_1 = false;
0313     }
0314 
0315     if (m_ntpc_1 >= m_nTPC_passprobe && m_nmaps_1 >= m_nMVTX_passprobe && m_nintt_1 >= m_nINTT_passprobe && 
0316         m_nmms_1 >= m_nTPOT_passprobe && m_quality_1 <= m_quality_passprobe)
0317     {
0318       m_probepass_1 = true; 
0319     }
0320     else
0321     {
0322       m_probepass_1 = false;
0323     }
0324     
0325     int t2 = 0;
0326     for (const auto &[key2, track2] : *m_trackmap)
0327     {
0328       if (!track2 || t1 >= t2 || track1->get_charge() == track2->get_charge() || track1->get_crossing() != track2->get_crossing())
0329       {
0330         ++t2;
0331         continue;
0332       }
0333       m_trackID_2 = track2->get_id(); 
0334       m_crossing_2 = track2->get_crossing();   
0335       m_px_2 = track2->get_px();
0336       m_py_2 = track2->get_py();
0337       m_pz_2 = track2->get_pz();
0338       m_pt_2 = track2->get_pt();
0339       m_eta_2 = track2->get_eta();
0340       m_phi_2 = track2->get_phi();
0341       m_charge_2 = track2->get_charge();
0342       m_quality_2 = track2->get_quality();
0343       m_chisq_2 = track2->get_chisq();
0344       m_ndf_2 = track2->get_ndf();
0345 
0346       TrackSeed* silseed2 = track2->get_silicon_seed();
0347       TrackSeed* tpcseed2 = track2->get_tpc_seed();   
0348       m_nmaps_2 = 0;
0349       m_nintt_2 = 0;
0350       m_ntpc_2 = 0;
0351       m_nmms_2 = 0;
0352       m_nhits_2 = 0; 
0353       maps.clear();
0354       intt.clear();
0355       tpc.clear();
0356       mms.clear();
0357 
0358       if (tpcseed2)
0359       {
0360         m_nhits_2 += tpcseed2->size_cluster_keys();
0361         for (TrackSeed::ConstClusterKeyIter local_iter = tpcseed2->begin_cluster_keys();
0362              local_iter != tpcseed2->end_cluster_keys();
0363              ++local_iter)
0364         {
0365           TrkrDefs::cluskey cluster_key = *local_iter;
0366           //TrkrCluster* cluster = clustermap->findCluster(cluster_key);
0367           unsigned int local_layer = TrkrDefs::getLayer(cluster_key);
0368           if (_nlayers_maps > 0 && local_layer < _nlayers_maps)
0369           {
0370             maps[local_layer] = 1;
0371             m_nmaps_2++;
0372           }
0373           if (_nlayers_intt > 0 && local_layer >= _nlayers_maps && local_layer < _nlayers_maps + _nlayers_intt)
0374           {
0375             intt[local_layer - _nlayers_maps] = 1;
0376             m_nintt_2++;
0377           }
0378           if (_nlayers_tpc > 0 &&
0379               local_layer >= (_nlayers_maps + _nlayers_intt) &&
0380               local_layer < (_nlayers_maps + _nlayers_intt + _nlayers_tpc))
0381           {
0382             tpc[local_layer - (_nlayers_maps + _nlayers_intt)] = 1;
0383             m_ntpc_2++; 
0384           }
0385           if (_nlayers_mms > 0 && local_layer >= _nlayers_maps + _nlayers_intt + _nlayers_tpc)
0386           {
0387             mms[local_layer - (_nlayers_maps + _nlayers_intt + _nlayers_tpc)] = 1;
0388             m_nmms_2++;
0389           }
0390         }
0391       }
0392       if (silseed2)
0393       {
0394         m_nhits_2 += silseed2->size_cluster_keys();
0395         for (TrackSeed::ConstClusterKeyIter local_iter = silseed2->begin_cluster_keys();
0396              local_iter != silseed2->end_cluster_keys();
0397              ++local_iter)
0398         {
0399           TrkrDefs::cluskey cluster_key = *local_iter;
0400           //TrkrCluster* cluster = clustermap->findCluster(cluster_key);
0401           unsigned int local_layer = TrkrDefs::getLayer(cluster_key);
0402           if (_nlayers_maps > 0 && local_layer < _nlayers_maps)
0403           {
0404             maps[local_layer] = 1;
0405             m_nmaps_2++;
0406           }
0407           if (_nlayers_intt > 0 && local_layer >= _nlayers_maps && local_layer < _nlayers_maps + _nlayers_intt)
0408           {
0409             intt[local_layer - _nlayers_maps] = 1;
0410             m_nintt_2++;
0411           }
0412           if (_nlayers_tpc > 0 &&
0413               local_layer >= (_nlayers_maps + _nlayers_intt) &&
0414               local_layer < (_nlayers_maps + _nlayers_intt + _nlayers_tpc))
0415           {
0416             tpc[local_layer - (_nlayers_maps + _nlayers_intt)] = 1;
0417             m_ntpc_2++; 
0418           }
0419           if (_nlayers_mms > 0 && local_layer >= _nlayers_maps + _nlayers_intt + _nlayers_tpc)
0420           {
0421             mms[local_layer - (_nlayers_maps + _nlayers_intt + _nlayers_tpc)] = 1;
0422             m_nmms_2++;
0423           }
0424         }
0425       }
0426       nlintt = 0;
0427       nlmaps = 0;
0428       nltpc = 0;
0429       nlmms = 0;
0430       if (_nlayers_maps > 0)
0431       {
0432         for (unsigned int i = 0; i < _nlayers_maps; i++)
0433         {
0434           nlmaps += maps[i];
0435         }
0436       }  
0437       if (_nlayers_intt > 0)
0438       {
0439         for (unsigned int i = 0; i < _nlayers_intt; i++)
0440         {
0441           nlintt += intt[i];
0442         }
0443       }
0444       if (_nlayers_tpc > 0)
0445       {
0446         for (unsigned int i = 0; i < _nlayers_tpc; i++)
0447         {
0448           nltpc += tpc[i];
0449         }
0450       }
0451       if (_nlayers_mms > 0)
0452       {
0453         for (unsigned int i = 0; i < _nlayers_mms; i++)
0454         {
0455           nlmms += mms[i];
0456         }
0457       }
0458       m_nlayers_2 = nlmaps + nlintt + nltpc + nlmms;
0459 
0460       m_mapsstates_2 = 0;
0461       m_inttstates_2 = 0;
0462       m_tpcstates_2 = 0;
0463       m_mmsstates_2 = 0;
0464       for (auto state_iter = track2->begin_states();
0465          state_iter != track2->end_states();
0466          ++state_iter)
0467       {
0468         SvtxTrackState* tstate = state_iter->second;
0469         if (tstate->get_pathlength() != 0) //The first track state is an extrapolation so has no cluster
0470         {
0471           auto stateckey = tstate->get_cluskey();
0472           uint8_t id = TrkrDefs::getTrkrId(stateckey);
0473     
0474           switch (id)
0475           { 
0476             case TrkrDefs::mvtxId:
0477               ++m_mapsstates_2;
0478               break;
0479             case TrkrDefs::inttId:
0480               ++m_inttstates_2;
0481               break;
0482             case TrkrDefs::tpcId:
0483               ++m_tpcstates_2;
0484               break;
0485             case TrkrDefs::micromegasId:
0486               ++m_mmsstates_2;
0487               break;
0488             default:
0489              std::cout << "Cluster key doesnt match a tracking system, this shouldn't happen" << std::endl;
0490              break; 
0491           }
0492         }
0493       }
0494 
0495       if (m_ntpc_2 >= m_nTPC_tag && m_nmaps_2 >= m_nMVTX_tag && m_nintt_2 >= m_nINTT_tag && 
0496           m_nmms_2 >= m_nTPOT_tag && m_quality_2 <= m_quality_tag)
0497       {
0498         m_tagpass_2 = true; 
0499       }
0500       else
0501       {
0502         m_tagpass_2 = false;
0503       }
0504 
0505       if (m_ntpc_2 >= m_nTPC_passprobe && m_nmaps_2 >= m_nMVTX_passprobe && m_nintt_2 >= m_nINTT_passprobe && 
0506           m_nmms_2 >= m_nTPOT_passprobe && m_quality_2 <= m_quality_passprobe)
0507       {
0508         m_probepass_2 = true; 
0509       }
0510       else
0511       {
0512         m_probepass_2 = false;
0513       }
0514     
0515       double pair_dca;
0516       Acts::Vector3 pca_rel1;
0517       Acts::Vector3 pca_rel2;
0518       // assuming stright line tracks for a rough estimate
0519       Acts::Vector3 A_pos1(track1->get_x(), track1->get_y(), track1->get_z());
0520       Acts::Vector3 A_mom1(track1->get_px(), track1->get_py(), track1->get_pz());
0521       Acts::Vector3 A_pos2(track2->get_x(), track2->get_y(), track2->get_z());
0522       Acts::Vector3 A_mom2(track2->get_px(), track2->get_py(), track2->get_pz());
0523       findPcaTwoTracks(A_pos1, A_pos2, A_mom1, A_mom2, pca_rel1, pca_rel2, pair_dca);     
0524 
0525       Eigen::Vector3d projected_pos1;
0526       Eigen::Vector3d projected_mom1;
0527       Eigen::Vector3d projected_pos2;
0528       Eigen::Vector3d projected_mom2;
0529 
0530       bool ret1 = projectTrackToPoint(track1, pca_rel1, projected_pos1, projected_mom1);
0531       bool ret2 = projectTrackToPoint(track2, pca_rel2, projected_pos2, projected_mom2);
0532      
0533 
0534       // Invariant Mass Calculation Info
0535       double E1 = sqrt(pow(projected_mom1(0), 2) + pow(projected_mom1(1), 2) + pow(projected_mom1(2), 2) + pow(_pionMass, 2));
0536       double E2 = sqrt(pow(projected_mom2(0), 2) + pow(projected_mom2(1), 2) + pow(projected_mom2(2), 2) + pow(_pionMass, 2));
0537       TLorentzVector tra1(projected_mom1(0), projected_mom1(1), projected_mom1(2), E1);
0538       TLorentzVector tra2(projected_mom2(0), projected_mom2(1), projected_mom2(2), E2);
0539       TLorentzVector tsum;
0540       tsum = tra1 + tra2;
0541       m_invM = tsum.M();
0542 
0543       if (!ret1 || !ret2)
0544       {
0545         if (Verbosity() > 5)
0546         {
0547           std::cout << "Failed to make track params: track1 - " << ret1 << ", track2 - " << ret2 << std::endl;
0548         }
0549       }  
0550 
0551       m_px_proj_1 = projected_mom1(0);
0552       m_py_proj_1 = projected_mom1(1);
0553       m_pz_proj_1 = projected_mom1(2);
0554       m_px_proj_2 = projected_mom2(0);
0555       m_py_proj_2 = projected_mom2(1);
0556       m_pz_proj_2 = projected_mom2(2);
0557 
0558       double pair_dca_proj;
0559       Acts::Vector3 pca_rel1_proj;
0560       Acts::Vector3 pca_rel2_proj;
0561       findPcaTwoTracks(projected_pos1, projected_pos2, projected_mom1, projected_mom2, pca_rel1_proj, pca_rel2_proj, pair_dca_proj);
0562 
0563       m_x_proj_1 = pca_rel1_proj(0);
0564       m_y_proj_1 = pca_rel1_proj(1);
0565       m_z_proj_1 = pca_rel1_proj(2);
0566       m_x_proj_2 = pca_rel2_proj(0);
0567       m_y_proj_2 = pca_rel2_proj(1);
0568       m_z_proj_2 = pca_rel2_proj(2);
0569       m_dca = pair_dca_proj;
0570       
0571       //KFParticle tr1_kfp = makeParticle(track1);
0572       //KFParticle tr2_kfp = makeParticle(track2);
0573       //KFParticle mother;
0574       //mother.AddDaughter(tr1_kfp);
0575       //mother.AddDaughter(tr2_kfp);
0576       //tr1_kfp.SetProductionVertex(mother);
0577       //tr2_kfp.SetProductionVertex(mother);
0578       //float calculated_mass_err;
0579       //mother.GetMass(m_invM_kfp, calculated_mass_err); 
0580 
0581       m_TAPTree->Fill();
0582  
0583       ++t2;
0584     } 
0585     clearValues();
0586     ++t1;
0587   } 
0588   
0589   ++m_event;
0590 
0591   return Fun4AllReturnCodes::EVENT_OK;
0592 }
0593 
0594 int TagAndProbe::End(PHCompositeNode * /*topNode*/)
0595 {
0596   std::cout << "TagAndProbe identification finished" << std::endl;
0597 
0598   if (m_save_output)
0599   {
0600     m_outfile->cd();
0601     m_cutInfoTree->Write();
0602     m_TAPTree->Write();
0603     m_outfile->Close();
0604   }
0605 
0606   return 0;
0607 }
0608 
0609 void TagAndProbe::initializeBranches()
0610 {
0611   m_outfile = new TFile(m_outfile_name.c_str(), "RECREATE");  
0612   m_cutInfoTree = new TTree("CutInfoTree", "CutInfoTree");
0613   m_TAPTree = new TTree("TAPTree","TAPTree");
0614 
0615   m_cutInfoTree->OptimizeBaskets();
0616   m_cutInfoTree->SetAutoSave(-5e6);  // Save the output file every 5MB
0617   m_TAPTree->OptimizeBaskets();
0618   m_TAPTree->SetAutoSave(-5e6);  // Save the output file every 5MB
0619 
0620   // Tree containing the cut info used for the entries in this root file
0621   m_cutInfoTree->Branch("tag_nTPC_cut", &m_nTPC_tag, "tag_nTPC_cut/F");
0622   m_cutInfoTree->Branch("tag_nINTT_cut", &m_nINTT_tag, "tag_nINTT_cut/F");
0623   m_cutInfoTree->Branch("tag_nMVTX_cut", &m_nMVTX_tag, "tag_nMVTX_cut/F");
0624   m_cutInfoTree->Branch("tag_nTPOT_cut", &m_nTPOT_tag, "tag_nTPOT_cut/F");
0625   m_cutInfoTree->Branch("tag_quality_cut", &m_quality_tag, "tag_quality_cut/F");
0626   m_cutInfoTree->Branch("probe_nTPC_cut", &m_nTPC_passprobe, "probe_nTPC_cut/F");
0627   m_cutInfoTree->Branch("probe_nINTT_cut", &m_nINTT_passprobe, "probe_nINTT_cut/F");
0628   m_cutInfoTree->Branch("probe_nMVTX_cut", &m_nMVTX_passprobe, "probe_nMVTX_cut/F");
0629   m_cutInfoTree->Branch("probe_nTPOT_cut", &m_nTPOT_passprobe, "probe_nTPOT_cut/F");
0630   m_cutInfoTree->Branch("probe_quality_cut", &m_quality_passprobe, "probe_quality_cut/F");
0631 
0632   // Tree containing all information for pairs of tracks
0633   // Includes a tag identifying if they passed or failed the cut being investigated 
0634   m_TAPTree->Branch("run", &m_run, "run/F");
0635   m_TAPTree->Branch("segment", &m_segment, "segment/F");
0636   m_TAPTree->Branch("event", &m_event, "event/F");
0637   //track 1
0638   m_TAPTree->Branch("trackID_1", &m_trackID_1, "trackID_1/F");
0639   m_TAPTree->Branch("crossing_1", &m_crossing_1, "crossing_1/F");
0640   m_TAPTree->Branch("px_1", &m_px_1, "px_1/F");
0641   m_TAPTree->Branch("py_1", &m_py_1, "py_1/F");
0642   m_TAPTree->Branch("pz_1", &m_pz_1, "pz_1/F");
0643   m_TAPTree->Branch("pt_1", &m_pt_1, "pt_1/F");
0644   m_TAPTree->Branch("px_proj_1", &m_px_proj_1, "px_proj_1/F");
0645   m_TAPTree->Branch("py_proj_1", &m_py_proj_1, "py_proj_1/F");
0646   m_TAPTree->Branch("pz_proj_1", &m_pz_proj_1, "pz_proj_1/F");
0647   //m_TAPTree->Branch("px_proj_kfp_1", &m_px_proj_kfp_1, "px_proj_kfp_1/F");
0648   //m_TAPTree->Branch("py_proj_kfp_1", &m_py_proj_kfp_1, "py_proj_kfp_1/F");
0649   //m_TAPTree->Branch("pz_proj_kfp_1", &m_pz_proj_kfp_1, "pz_proj_kfp_1/F");
0650   m_TAPTree->Branch("x_proj_1", &m_x_proj_1, "x_proj_1/F");
0651   m_TAPTree->Branch("y_proj_1", &m_y_proj_1, "y_proj_1/F");
0652   m_TAPTree->Branch("z_proj_1", &m_z_proj_1, "z_proj_1/F");
0653   m_TAPTree->Branch("eta_1", &m_eta_1, "eta_1/F");
0654   m_TAPTree->Branch("phi_1", &m_phi_1, "phi_1/F");
0655   m_TAPTree->Branch("charge_1", &m_charge_1, "charge_1/F");
0656   m_TAPTree->Branch("quality_1", &m_quality_1, "quality_1/F");
0657   m_TAPTree->Branch("chisq_1", &m_chisq_1, "chisq_1/F");
0658   m_TAPTree->Branch("ndf_1", &m_ndf_1, "ndf_1/F");
0659   m_TAPTree->Branch("nhits_1", &m_nhits_1, "nhits_1/F");
0660   m_TAPTree->Branch("nlayers_1", &m_nlayers_1, "nlayers_1/F");
0661   m_TAPTree->Branch("nmaps_1", &m_nmaps_1, "nmaps_1/F");
0662   m_TAPTree->Branch("nintt_1", &m_nintt_1, "nintt_1/F");
0663   m_TAPTree->Branch("ntpc_1", &m_ntpc_1, "ntpc_1/F");
0664   m_TAPTree->Branch("nmms_1", &m_nmms_1, "nmms_1/F");
0665   m_TAPTree->Branch("mapsstates_1", &m_mapsstates_1, "mapsstates_1/F");
0666   m_TAPTree->Branch("inttstates_1", &m_inttstates_1, "inttstates_1/F");
0667   m_TAPTree->Branch("tpcstates_1", &m_tpcstates_1, "tpcstates_1/F");
0668   m_TAPTree->Branch("mmsstates_1", &m_mmsstates_1, "mmsstates_1/F");
0669   m_TAPTree->Branch("pca_x_1", &m_pca_x_1, "pca_x_1/F");
0670   m_TAPTree->Branch("pca_y_1", &m_pca_y_1, "pca_y_1/F");
0671   m_TAPTree->Branch("pca_z_1", &m_pca_z_1, "pca_z_1/F");
0672   m_TAPTree->Branch("tagpass_1", &m_tagpass_1, "tagpass_1/F");
0673   m_TAPTree->Branch("probepass_1", &m_probepass_1, "probepass_1/F");
0674   //track 2
0675   m_TAPTree->Branch("trackID_2", &m_trackID_2, "trackID_2/F");
0676   m_TAPTree->Branch("crossing_2", &m_crossing_2, "crossing_2/F");
0677   m_TAPTree->Branch("px_2", &m_px_2, "px_2/F");
0678   m_TAPTree->Branch("py_2", &m_py_2, "py_2/F");
0679   m_TAPTree->Branch("pz_2", &m_pz_2, "pz_2/F");
0680   m_TAPTree->Branch("pt_2", &m_pt_2, "pt_2/F");
0681   m_TAPTree->Branch("px_proj_2", &m_px_proj_2, "px_proj_2/F");
0682   m_TAPTree->Branch("py_proj_2", &m_py_proj_2, "py_proj_2/F");
0683   m_TAPTree->Branch("pz_proj_2", &m_pz_proj_2, "pz_proj_2/F");
0684   //m_TAPTree->Branch("px_proj_kfp_2", &m_px_proj_kfp_2, "px_proj_kfp_2/F");
0685   //m_TAPTree->Branch("py_proj_kfp_2", &m_py_proj_kfp_2, "py_proj_kfp_2/F");
0686   //m_TAPTree->Branch("pz_proj_kfp_2", &m_pz_proj_kfp_2, "pz_proj_kfp_2/F");
0687   m_TAPTree->Branch("x_proj_2", &m_x_proj_2, "x_proj_2/F");
0688   m_TAPTree->Branch("y_proj_2", &m_y_proj_2, "y_proj_2/F");
0689   m_TAPTree->Branch("z_proj_2", &m_z_proj_2, "z_proj_2/F");
0690   m_TAPTree->Branch("eta_2", &m_eta_2, "eta_2/F");
0691   m_TAPTree->Branch("phi_2", &m_phi_2, "phi_2/F");
0692   m_TAPTree->Branch("charge_2", &m_charge_2, "charge_2/F");
0693   m_TAPTree->Branch("quality_2", &m_quality_2, "quality_2/F");
0694   m_TAPTree->Branch("chisq_2", &m_chisq_2, "chisq_2/F");
0695   m_TAPTree->Branch("ndf_2", &m_ndf_2, "ndf_2/F");
0696   m_TAPTree->Branch("nhits_2", &m_nhits_2, "nhits_2/F");
0697   m_TAPTree->Branch("nlayers_2", &m_nlayers_2, "nlayers_2/F");
0698   m_TAPTree->Branch("nmaps_2", &m_nmaps_2, "nmaps_2/F");
0699   m_TAPTree->Branch("nintt_2", &m_nintt_2, "nintt_2/F");
0700   m_TAPTree->Branch("ntpc_2", &m_ntpc_2, "ntpc_2/F");
0701   m_TAPTree->Branch("nmms_2", &m_nmms_2, "nmms_2/F");
0702   m_TAPTree->Branch("mapsstates_2", &m_mapsstates_2, "mapsstates_2/F");
0703   m_TAPTree->Branch("inttstates_2", &m_inttstates_2, "inttstates_2/F");
0704   m_TAPTree->Branch("tpcstates_2", &m_tpcstates_2, "tpcstates_2/F");
0705   m_TAPTree->Branch("mmsstates_2", &m_mmsstates_2, "mmsstates_2/F");
0706   m_TAPTree->Branch("pca_x_2", &m_pca_x_2, "pca_x_2/F");
0707   m_TAPTree->Branch("pca_y_2", &m_pca_y_2, "pca_y_2/F");
0708   m_TAPTree->Branch("pca_z_2", &m_pca_z_2, "pca_z_2/F");
0709   m_TAPTree->Branch("tagpass_2", &m_tagpass_2, "tagpass_2/F");
0710   m_TAPTree->Branch("probepass_2", &m_probepass_2, "probepass_2/F");
0711   //both tracks
0712   m_TAPTree->Branch("invM", &m_invM, "invM/F");
0713   //m_TAPTree->Branch("invM_kfp", &m_invM_kfp, "invM_kfp/F");
0714   m_TAPTree->Branch("dca", &m_dca, "dca/F");
0715   if (m_include_pv_info)
0716   {
0717     m_TAPTree->Branch("vx", &m_vx, "vx/F");
0718     m_TAPTree->Branch("vy", &m_vy, "vy/F");
0719     m_TAPTree->Branch("vz", &m_vz, "vz/F");
0720   }
0721 }
0722 
0723 void TagAndProbe::clearValues()
0724 {
0725   m_trackID_1 = NAN; m_trackID_2 = NAN;
0726   m_crossing_1 = NAN; m_crossing_2 = NAN;
0727   m_px_1 = NAN; m_px_2 = NAN;
0728   m_py_1 = NAN; m_py_2 = NAN;
0729   m_pz_1 = NAN; m_pz_2 = NAN;
0730   m_px_proj_1 = NAN; m_px_proj_2 = NAN;
0731   m_py_proj_1 = NAN; m_py_proj_2 = NAN;
0732   m_pz_proj_1 = NAN; m_pz_proj_2 = NAN;
0733   //m_px_proj_kfp_1 = NAN; m_px_proj_kfp_2 = NAN;
0734   //m_py_proj_kfp_1 = NAN; m_py_proj_kfp_2 = NAN;
0735   //m_pz_proj_kfp_1 = NAN; m_pz_proj_kfp_2 = NAN;
0736   m_x_proj_1 = NAN; m_x_proj_2 = NAN;
0737   m_y_proj_1 = NAN; m_y_proj_2 = NAN;
0738   m_z_proj_1 = NAN; m_z_proj_2 = NAN;
0739   m_eta_1 = NAN; m_eta_2 = NAN;
0740   m_phi_1 = NAN; m_phi_2 = NAN;
0741   m_charge_1 = NAN; m_charge_2 = NAN;
0742   m_quality_1 = NAN; m_quality_2 = NAN;
0743   m_chisq_1 = NAN; m_chisq_2 = NAN;
0744   m_ndf_1 = NAN; m_ndf_2 = NAN;
0745   m_nhits_1 = NAN; m_nhits_2 = NAN;
0746   m_nlayers_1 = NAN; m_nlayers_2 = NAN;
0747   m_nmaps_1 = NAN; m_nmaps_2 = NAN;
0748   m_nintt_1 = NAN; m_nintt_2 = NAN;
0749   m_ntpc_1 = NAN; m_ntpc_2 = NAN;
0750   m_nmms_1 = NAN; m_nmms_2 = NAN;
0751   m_mapsstates_1 = NAN; m_mapsstates_2 = NAN;
0752   m_inttstates_1 = NAN; m_inttstates_2 = NAN;
0753   m_tpcstates_1 = NAN; m_tpcstates_2 = NAN;
0754   m_mmsstates_1 = NAN; m_mmsstates_2 = NAN;
0755   m_pca_x_1 = NAN; m_pca_x_2 = NAN;
0756   m_pca_y_1 = NAN; m_pca_y_2 = NAN;
0757   m_pca_z_1 = NAN; m_pca_z_2 = NAN;
0758   m_tagpass_1 = NAN; m_tagpass_2 = NAN;
0759   m_probepass_1 = NAN; m_probepass_2 = NAN;
0760 
0761   m_invM = NAN;
0762   //m_invM_kfp = NAN;
0763   m_dca = NAN;
0764   if (m_include_pv_info)
0765   {
0766     m_vx = NAN;
0767     m_vy = NAN;
0768     m_vz = NAN;
0769   } 
0770 }
0771 
0772 // borrowed from KshortReconstruction
0773 bool TagAndProbe::projectTrackToPoint(SvtxTrack* track, Eigen::Vector3d PCA, Eigen::Vector3d& pos, Eigen::Vector3d& mom)
0774 {
0775   bool ret = true;
0776 
0777   track->identify();
0778 
0779   /// create perigee surface
0780   ActsPropagator actsPropagator(m_tGeometry);
0781   auto perigee = actsPropagator.makeVertexSurface(PCA);  // PCA is in cm here
0782   auto params = actsPropagator.makeTrackParams(track, m_vertexMap);
0783   if (!params.ok())
0784   {
0785     return false;
0786   }
0787   auto result = actsPropagator.propagateTrack(params.value(), perigee);
0788 
0789   if (result.ok())
0790   {
0791     auto projectionPos = result.value().second.position(m_tGeometry->geometry().getGeoContext());
0792     const auto momentum = result.value().second.momentum();
0793     pos(0) = projectionPos.x() / Acts::UnitConstants::cm;
0794     pos(1) = projectionPos.y() / Acts::UnitConstants::cm;
0795     pos(2) = projectionPos.z() / Acts::UnitConstants::cm;
0796 
0797     if (Verbosity() > 2)
0798     {
0799       std::cout << "                 Input PCA " << PCA << "  projection out " << pos << std::endl;
0800     }
0801 
0802     mom(0) = momentum.x();
0803     mom(1) = momentum.y();
0804     mom(2) = momentum.z();
0805   }
0806   else
0807   {
0808     pos(0) = track->get_x();
0809     pos(1) = track->get_y();
0810     pos(2) = track->get_z();
0811 
0812     mom(0) = track->get_px();
0813     mom(1) = track->get_py();
0814     mom(2) = track->get_pz();
0815 
0816     if(Verbosity() > 0)
0817       {
0818     std::cout << result.error() << std::endl;
0819     std::cout << result.error().message() << std::endl;
0820     std::cout << " Failed projection of track with: " << std::endl;
0821     std::cout << " x,y,z = " << track->get_x() << "  " << track->get_y() << "  " << track->get_z() << std::endl;
0822     std::cout << " px,py,pz = " << track->get_px() << "  " << track->get_py() << "  " << track->get_pz() << std::endl;
0823     std::cout << " to point (x,y,z) = " << PCA(0) / Acts::UnitConstants::cm << "  " << PCA(1) / Acts::UnitConstants::cm << "  " << PCA(2) / Acts::UnitConstants::cm << std::endl;
0824       }
0825 
0826     //    ret = false;
0827   }
0828 
0829   return ret;
0830 }
0831 
0832 // borrowed from KshortReconstruction
0833 void TagAndProbe::findPcaTwoTracks(const Acts::Vector3& pos1, const Acts::Vector3& pos2, Acts::Vector3 mom1, Acts::Vector3 mom2, Acts::Vector3& pca1, Acts::Vector3& pca2, double& dca)
0834 {
0835   TLorentzVector v1;
0836   TLorentzVector v2;
0837 
0838   double px1 = mom1(0);
0839   double py1 = mom1(1);
0840   double pz1 = mom1(2);
0841   double px2 = mom2(0);
0842   double py2 = mom2(1);
0843   double pz2 = mom2(2);
0844 
0845   Float_t E1 = sqrt(pow(px1, 2) + pow(py1, 2) + pow(pz1, 2) + pow(_pionMass, 2));
0846   Float_t E2 = sqrt(pow(px2, 2) + pow(py2, 2) + pow(pz2, 2) + pow(_pionMass, 2));
0847 
0848   v1.SetPxPyPzE(px1, py1, pz1, E1);
0849   v2.SetPxPyPzE(px2, py2, pz2, E2);
0850 
0851   // calculate lorentz vector
0852   const Eigen::Vector3d& a1 = pos1;
0853   const Eigen::Vector3d& a2 = pos2;
0854 
0855   Eigen::Vector3d b1(v1.Px(), v1.Py(), v1.Pz());
0856   Eigen::Vector3d b2(v2.Px(), v2.Py(), v2.Pz());
0857 
0858   // The shortest distance between two skew lines described by
0859   //  a1 + c * b1
0860   //  a2 + d * b2
0861   // where a1, a2, are vectors representing points on the lines, b1, b2 are direction vectors, and c and d are scalars
0862   // dca = (b1 x b2) .(a2-a1) / |b1 x b2|
0863 
0864   // bcrossb/mag_bcrossb is a unit vector perpendicular to both direction vectors b1 and b2
0865   auto bcrossb = b1.cross(b2);
0866   auto mag_bcrossb = bcrossb.norm();
0867   // a2-a1 is the vector joining any arbitrary points on the two lines
0868   auto aminusa = a2 - a1;
0869 
0870   // The DCA of these two lines is the projection of a2-a1 along the direction of the perpendicular to both
0871   // remember that a2-a1 is longer than (or equal to) the dca by definition
0872   dca = 999;
0873   if (mag_bcrossb != 0)
0874   {
0875     dca = bcrossb.dot(aminusa) / mag_bcrossb;
0876   }
0877   else
0878   {
0879     return;  // same track, skip combination
0880   }
0881   
0882   // get the points at which the normal to the lines intersect the lines, where the lines are perpendicular
0883   double X = b1.dot(b2) - b1.dot(b1) * b2.dot(b2) / b2.dot(b1);
0884   double Y = (a2.dot(b2) - a1.dot(b2)) - (a2.dot(b1) - a1.dot(b1)) * b2.dot(b2) / b2.dot(b1);
0885   double c = Y / X;
0886 
0887   double F = b1.dot(b1) / b2.dot(b1);
0888   double G = -(a2.dot(b1) - a1.dot(b1)) / b2.dot(b1);
0889   double d = c * F + G;
0890 
0891   // then the points of closest approach are:
0892   pca1 = a1 + c * b1;
0893   pca2 = a2 + d * b2;
0894 
0895   return;
0896 }
0897 
0898 /*
0899 // borrowed from KFParticle_sPHENIX
0900 void TagAndProbe::getField()
0901 {
0902   //This sweeps the sPHENIX magnetic field map from some point radially then grabs the first event that passes the selection
0903   m_magField = std::filesystem::exists(m_magField) ? m_magField : CDBInterface::instance()->getUrl(m_magField);
0904 
0905   if (Verbosity() > 0)
0906   {
0907     std::cout << PHWHERE << ": using fieldmap : " << m_magField << std::endl;
0908   }
0909 
0910   TFile *fin = new TFile(m_magField.c_str());
0911   TTree *fieldmap = (TTree *) fin->Get("fieldmap");
0912 
0913   float Bz = 0.;
0914   unsigned int r = 0.;
0915   float z = 0.;
0916 
0917   double arc = M_PI/2;
0918   unsigned int n = 0;
0919 
0920   while (Bz == 0)
0921   {
0922     if (n == 4)
0923     {
0924       ++r;
0925     }
0926 
0927     if (r == 3) //Dont go too far out radially
0928     {
0929       ++z;
0930     }
0931 
0932     n = n & 0x3U; //Constrains n from 0 to 3
0933     r = r & 0x2U;
0934 
0935     double x = r*std::cos(n*arc);
0936     double y = r*std::sin(n*arc);
0937 
0938     std::string sweep = "x == " + std::to_string(x) + " && y == " + std::to_string(y) + " && z == " + std::to_string(z);
0939 
0940     fieldmap->Draw(">>elist", sweep.c_str(), "entrylist");
0941     TEntryList *elist = (TEntryList*)gDirectory->Get("elist");
0942     if (elist->GetEntry(0))
0943     {
0944       TLeaf *fieldValue = fieldmap->GetLeaf("bz");
0945       fieldValue->GetBranch()->GetEntry(elist->GetEntry(0));
0946       Bz = fieldValue->GetValue();
0947     }
0948 
0949     ++n;
0950 
0951     if (r == 0) // No point in rescanning (0,0)
0952     {
0953       ++r;
0954       n = 0;
0955     }
0956   }
0957   // The actual unit of KFParticle is in kilo Gauss (kG), which is equivalent to 0.1 T, instead of Tesla (T). The positive value indicates the B field is in the +z direction
0958   Bz *= 10;  // Factor of 10 to convert the B field unit from kG to T
0959   KFParticle::SetField((double) Bz);
0960 
0961   fieldmap->Delete();
0962   fin->Close();
0963 }
0964 
0965 
0966 // borrowed from KFParticle_sPHENIX
0967 KFParticle TagAndProbe::makeParticle(SvtxTrack* track_kfp)  /// Return a KFPTrack from track vector and covariance matrix. No mass or vertex constraints
0968 {
0969   KFParticle kfp_particle;
0970 
0971   float f_trackParameters[6] = {track_kfp->get_x(),
0972                                 track_kfp->get_y(),
0973                                 track_kfp->get_z(),
0974                                 track_kfp->get_px(),
0975                                 track_kfp->get_py(),
0976                                 track_kfp->get_pz()};
0977 
0978   float f_trackCovariance[21];
0979   unsigned int iterate = 0;
0980   for (unsigned int i = 0; i < 6; ++i)
0981   {
0982     for (unsigned int j = 0; j <= i; ++j)
0983     {
0984       f_trackCovariance[iterate] = track_kfp->get_error(i, j);
0985       ++iterate;
0986     }
0987   }
0988 
0989   kfp_particle.Create(f_trackParameters, f_trackCovariance, (Int_t) track_kfp->get_charge(), -1);
0990   kfp_particle.NDF() = track_kfp->get_ndf();
0991   kfp_particle.Chi2() = track_kfp->get_chisq();
0992   kfp_particle.SetId(track_kfp->get_id());
0993 
0994   return kfp_particle;
0995 }
0996 */
[ 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!