g4simulation/g4eval/FillClusMatchTree.cc

0001 #include "FillClusMatchTree.h"
0002 #include "TrkrClusterIsMatcher.h"
0003 #include "g4evalfn.h"
0004
0005 #include <g4tracking/EmbRecoMatch.h>
0006 #include <g4tracking/EmbRecoMatchContainer.h>
0007 #include <g4tracking/TrkrTruthTrack.h>
0008 #include <g4tracking/TrkrTruthTrackContainer.h>
0009
0010 #include <trackbase/ClusHitsVerbose.h>
0011 #include <trackbase/TrkrCluster.h>
0012 #include <trackbase/TrkrClusterContainer.h>
0013 #include <trackbase/TrkrDefs.h>
0014 #include <trackbase_historic/SvtxPHG4ParticleMap_v1.h>
0015 #include <trackbase_historic/SvtxTrack.h>
0016 #include <trackbase_historic/SvtxTrackMap.h>
0017
0018 #include <g4main/PHG4TruthInfoContainer.h>
0019
0020 #include <fun4all/Fun4AllReturnCodes.h>
0021 #include <fun4all/PHTFileServer.h>
0022
0023 #include <phool/PHCompositeNode.h>
0024 #include <phool/PHDataNode.h>
0025 #include <phool/PHIODataNode.h>
0026 #include <phool/PHNode.h>
0027 #include <phool/PHNodeIterator.h>
0028 #include <phool/PHObject.h>  // for PHObject
0029 #include <phool/getClass.h>
0030 #include <phool/phool.h>  // for PHWHERE
0031
0032 #include <TFile.h>
0033 #include <TH2.h>
0034 #include <TTree.h>
0035
0036 #include <boost/format.hpp>
0037
0038 #include <iostream>
0039
0040 //____________________________________________________________________________..
0041 FillClusMatchTree::FillClusMatchTree(
0042     TrkrClusterIsMatcher* _ismatcher, const std::string& _outfile_name, bool _fill_clusters, bool _fill_clusverbose, bool _fill_SvUnMatched)
0043   : m_ismatcher{_ismatcher}
0044   , m_outfile_name{_outfile_name}
0045   , m_fill_clusters{_fill_clusters}
0046   , m_fill_clusverbose{_fill_clusverbose}
0047   , m_fill_SvUnmatched{_fill_SvUnMatched}
0048 {
0049   m_TCEval.ismatcher = m_ismatcher;
0050
0051   PHTFileServer::get().open(m_outfile_name, "RECREATE");
0052
0053   h2_G4_nPixelsPhi = new TH2D("G4_nPixelsPhi", "PHG4 Emb Tracks; cluster pixel width Phi; layer",
0054                               100, -0.5, 99.5, 56, -0.5, 55.5);
0055   h2_G4_nPixelsZ = new TH2D("G4_nPixelsZ", "PHG4 Emb Tracks; cluster pixel width Z; layer",
0056                             100, -0.5, 99.5, 56, -0.5, 55.5);
0057   h2_Sv_nPixelsPhi = new TH2D("Sv_nPixelsPhi", "Svtx Reco Tracks; cluster pixel width Phi; layer",
0058                               100, -0.5, 99.5, 56, -0.5, 55.5);
0059   h2_Sv_nPixelsZ = new TH2D("Sv_nPixelsZ", "Svtx Reco Tracks; cluster pixel width Z; layer",
0060                             100, -0.5, 99.5, 56, -0.5, 55.5);
0061
0062   m_ttree = new TTree("T", "Tracks (and sometimes clusters)");
0063
0064   m_ttree->Branch("event", &nevent);
0065   m_ttree->Branch("nphg4_part", &nphg4_part);
0066   m_ttree->Branch("centrality", &centrality);
0067   m_ttree->Branch("ntrackmatches", &ntrackmatches);
0068   m_ttree->Branch("nphg4", &nphg4);
0069   m_ttree->Branch("nsvtx", &nsvtx);
0070
0071   m_ttree->Branch("trackid", &b_trackid);
0072   m_ttree->Branch("is_G4track", &b_is_g4track);
0073   m_ttree->Branch("is_Svtrack", &b_is_Svtrack);
0074   m_ttree->Branch("is_matched", &b_is_matched);
0075
0076   m_ttree->Branch("trkpt", &b_trkpt);
0077   m_ttree->Branch("trketa", &b_trketa);
0078   m_ttree->Branch("trkphi", &b_trkphi);
0079
0080   m_ttree->Branch("nclus", &b_nclus);
0081   m_ttree->Branch("nclustpc", &b_nclustpc);
0082   m_ttree->Branch("nclusmvtx", &b_nclusmvtx);
0083   m_ttree->Branch("nclusintt", &b_nclusintt);
0084   m_ttree->Branch("matchrat", &b_matchrat);
0085   m_ttree->Branch("matchrat_intt", &b_matchrat_intt);
0086   m_ttree->Branch("matchrat_mvtx", &b_matchrat_mvtx);
0087   m_ttree->Branch("matchrat_tpc", &b_matchrat_tpc);
0088
0089   if (m_fill_clusters)
0090   {
0091     m_ttree->Branch("clus_match", &b_clusmatch);
0092     m_ttree->Branch("clus_x", &b_clus_x);
0093     m_ttree->Branch("clus_y", &b_clus_y);
0094     m_ttree->Branch("clus_z", &b_clus_z);
0095     m_ttree->Branch("clus_r", &b_clus_r);
0096     m_ttree->Branch("clus_layer", &b_clus_layer);
0097     /* m_ttree ->Branch("nphibins"     , &b_clus_nphibins ); */
0098     /* m_ttree ->Branch("nzbins"       , &b_clus_ntbins   ); */
0099
0100     m_ttree->Branch("clus_lphi", &b_clus_lphi);
0101     m_ttree->Branch("clus_lphisize", &b_clus_lphisize);
0102     m_ttree->Branch("clus_lz", &b_clus_lz);
0103     m_ttree->Branch("clus_lzsize", &b_clus_lzsize);
0104
0105     if (m_fill_clusverbose)
0106     {
0107       m_ttree->Branch("phibins", &b_phibins);
0108       m_ttree->Branch("phibins_cut", &b_phibins_cut);
0109       m_ttree->Branch("zbins", &b_zbins);
0110       m_ttree->Branch("zbins_cut", &b_zbins_cut);
0111
0112       m_ttree->Branch("phibinsE", &b_phibinsE);
0113       m_ttree->Branch("phibinsE_cut", &b_phibinsE_cut);
0114       m_ttree->Branch("zbinsE", &b_zbinsE);
0115       m_ttree->Branch("zbinsE_cut", &b_zbinsE_cut);
0116     }
0117   }
0118 }
0119
0120 //____________________________________________________________________________..
0121 int FillClusMatchTree::Init(PHCompositeNode* topNode)
0122 {
0123   if (Verbosity() > 1)
0124   {
0125     std::cout << " Beginning FillClusMatchTree " << std::endl;
0126     topNode->print();
0127   }
0128
0129   return Fun4AllReturnCodes::EVENT_OK;
0130 }
0131
0132 //____________________________________________________________________________..
0133 int FillClusMatchTree::InitRun(PHCompositeNode* topNode)
0134 {
0135   /* auto init_status = m_cluster_comp.init(topNode); */
0136   auto init_status = m_ismatcher->init(topNode);
0137   if (init_status == Fun4AllReturnCodes::ABORTRUN)
0138   {
0139     return init_status;
0140   }
0141
0142   if (createNodes(topNode) != Fun4AllReturnCodes::EVENT_OK)
0143   {
0144     return Fun4AllReturnCodes::ABORTEVENT;
0145   }
0146
0147   return Fun4AllReturnCodes::EVENT_OK;
0148 }
0149
0150 int FillClusMatchTree::createNodes(PHCompositeNode* topNode)
0151 {
0152   PHNodeIterator iter(topNode);
0153
0154   PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0155
0156   if (!dstNode)
0157   {
0158     std::cout << PHWHERE << " DST node is missing, quitting" << std::endl;
0159     std::cout << PHWHERE << " DST node is missing, quitting" << std::endl;
0160     throw std::runtime_error("Failed to find DST node in FillClusMatchTree::createNodes");
0161   }
0162   topNode->print();
0163
0164   m_PHG4ClusHitVerb = findNode::getClass<ClusHitsVerbose>(topNode, "Trkr_TruthClusHitsVerbose");
0165   if (m_PHG4ClusHitVerb)
0166   {
0167     std::cout << " Found truth cluster hits verbose " << std::endl;
0168   }
0169   else
0170   {
0171     std::cout << " Did not find truth cluster hits verbose " << std::endl;
0172   }
0173
0174   m_SvtxClusHitVerb = findNode::getClass<ClusHitsVerbose>(topNode, "Trkr_SvtxClusHitsVerbose");
0175   if (m_SvtxClusHitVerb)
0176   {
0177     std::cout << " Found Svtx cluster hits verbose " << std::endl;
0178   }
0179   else
0180   {
0181     std::cout << " Did not find Svtx cluster hits verbose " << std::endl;
0182   }
0183
0184   m_EmbRecoMatchContainer = findNode::getClass<EmbRecoMatchContainer>(topNode, "TRKR_EMBRECOMATCHCONTAINER");
0185   if (!m_EmbRecoMatchContainer)
0186   {
0187     std::cout << PHWHERE << " Cannot find node TRKR_EMBRECOMATCHCONTAINER on node tree; quitting " << std::endl;
0188     std::cout << PHWHERE << " Cannot find node TRKR_EMBRECOMATCHCONTAINER on node tree; quitting " << std::endl;
0189     throw std::runtime_error(" Cannot find node TRKR_EMBRECOMATCHCONTAINER on node tree; quitting");
0190   }
0191
0192   PHCompositeNode* svtxNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "SVTX"));
0193   if (!svtxNode)
0194   {
0195     svtxNode = new PHCompositeNode("SVTX");
0196     dstNode->addNode(svtxNode);
0197   }
0198
0199   m_PHG4TruthInfoContainer = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0200   if (!m_PHG4TruthInfoContainer)
0201   {
0202     std::cout << "Could not locate G4TruthInfo node when running "
0203               << "\"TruthRecoTrackMatching\" module." << std::endl;
0204     return Fun4AllReturnCodes::ABORTEVENT;
0205   }
0206
0207   m_SvtxTrackMap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0208   if (!m_SvtxTrackMap)
0209   {
0210     std::cout << "Could not locate SvtxTrackMap node when running "
0211               << "\"TruthRecoTrackMatching\" module." << std::endl;
0212     return Fun4AllReturnCodes::ABORTEVENT;
0213   }
0214
0215   m_TrkrTruthTrackContainer = findNode::getClass<TrkrTruthTrackContainer>(topNode,
0216                                                                           "TRKR_TRUTHTRACKCONTAINER");
0217   if (!m_TrkrTruthTrackContainer)
0218   {
0219     std::cout << "Could not locate TRKR_TRUTHTRACKCONTAINER node when running "
0220               << "\"TruthRecoTrackMatching\" module." << std::endl;
0221     return Fun4AllReturnCodes::ABORTEVENT;
0222   }
0223
0224   return Fun4AllReturnCodes::EVENT_OK;
0225 }
0226
0227 int FillClusMatchTree::process_event(PHCompositeNode* /*topNode*/)
0228 {
0229   // We know that the data is in there and exists... why isn't it being transferred?
0230
0231   if (Verbosity() > 5)
0232   {
0233     std::cout << " FillClusMatchTree::process_event() " << std::endl;
0234   }
0235
0236   // fill in the event data
0237   ++nevent;
0238   nphg4 = m_TrkrTruthTrackContainer->getMap().size();
0239   nsvtx = m_SvtxTrackMap->size();
0240   ntrackmatches = m_EmbRecoMatchContainer->getMatches().size();
0241   // get centrality later...
0242
0243   // fill in pixel widths on truth tracks
0244   for (auto hitsetkey : m_TCEval.get_PHG4_clusters()->getHitSetKeys())
0245   {
0246     float layer = (float) TrkrDefs::getLayer(hitsetkey);
0247     auto range = m_TCEval.get_PHG4_clusters()->getClusters(hitsetkey);
0248     for (auto iter = range.first; iter != range.second; ++iter)
0249     {
0250       auto& cluster = iter->second;
0251       h2_G4_nPixelsPhi->Fill((float) cluster->getPhiSize(), layer);
0252       h2_G4_nPixelsZ->Fill((float) cluster->getZSize(), layer);
0253     }
0254   }
0255   // fill in pixel widths on reco tracks
0256   for (auto hitsetkey : m_TCEval.get_SVTX_clusters()->getHitSetKeys())
0257   {
0258     float layer = (float) TrkrDefs::getLayer(hitsetkey);
0259     auto range = m_TCEval.get_SVTX_clusters()->getClusters(hitsetkey);
0260     for (auto iter = range.first; iter != range.second; ++iter)
0261     {
0262       auto& cluster = iter->second;
0263       h2_Sv_nPixelsPhi->Fill((float) cluster->getPhiSize(), layer);
0264       h2_Sv_nPixelsZ->Fill((float) cluster->getZSize(), layer);
0265     }
0266   }
0267
0268   nphg4_part = 0;
0269   const auto range = m_PHG4TruthInfoContainer->GetPrimaryParticleRange();
0270   for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
0271        iter != range.second; ++iter)
0272   {
0273     nphg4_part++;
0274   }
0275
0276   // unmatched tracks are only entered once
0277   // matches can repeat a given svtx or phg4 track, depending on the
0278   // parameters in teh matching in FillClusMatchTree
0279   //
0280   // (1) fill unmatched phg4
0281   // (2) fill unmatched svtx
0282   // (3) fill matched phg4 and svtx
0283   clear_clusvecs(" nothing ");
0284
0285   if (Verbosity() > 2)
0286   {
0287     std::cout << " getting " << (int) m_EmbRecoMatchContainer->getMatches().size() << std::endl;
0288   }
0289   for (auto& match : m_EmbRecoMatchContainer->getMatches())
0290   {
0291     unsigned int g4_trkid = match->idTruthTrack();
0292     int sv_trkid = match->idRecoTrack();
0293
0294     auto g4trk = m_TrkrTruthTrackContainer->getTruthTrack(g4_trkid);
0295     auto svtrk = m_SvtxTrackMap->get(sv_trkid);
0296
0297     m_TCEval.addClusKeys(g4trk);
0298     m_TCEval.addClusKeys(svtrk);
0299     m_TCEval.find_matches();
0300
0301     // <- <- <- <- G4 Matched Tracks
0302     b_is_matched = true;
0303     b_is_g4track = true;
0304     b_is_Svtrack = false;
0305
0306     b_trackid = g4_trkid;
0307     b_trkpt = g4trk->getPt();
0308     b_trketa = g4trk->getPseudoRapidity();
0309     b_trkphi = g4trk->getPhi();
0310
0311     auto cnt = m_TCEval.phg4_cntclus();
0312     auto cnt_match = m_TCEval.phg4_cnt_matchedclus();
0313
0314     b_nclus = cnt[4];
0315     b_nclusmvtx = cnt[0];
0316     b_nclusintt = cnt[1];
0317     b_nclustpc = cnt[2];
0318
0319     b_matchrat = (float) cnt_match[4] / cnt[4];
0320     b_matchrat_mvtx = (float) cnt_match[0] / cnt[0];
0321     b_matchrat_intt = (float) cnt_match[1] / cnt[1];
0322     b_matchrat_tpc = (float) cnt_match[2] / cnt[2];
0323
0324     if (m_fill_clusters)
0325     {
0326       // - - - - UNMATCHED PHG4 hits in matched track
0327       auto clusters = m_TCEval.phg4_clusloc_unmatched();
0328       for (auto& clus : clusters)
0329       {
0330         b_clusmatch.push_back(false);
0331         fill_cluster_branches(clus, true);
0332       }
0333
0334       // - - - - MATCHED PHG4 hits in matched track
0335       clusters = m_TCEval.clusloc_matched();
0336       for (auto& clus : clusters)
0337       {
0338         b_clusmatch.push_back(true);
0339         fill_cluster_branches(clus, true);
0340       }
0341     }
0342     m_ttree->Fill();
0343     clear_clusvecs();
0344     /* clear_clusvecs("apple0 g4_matched"); */
0345
0346     // <- <- <- <- Svtx Matched Tracks
0347     b_is_g4track = false;
0348     b_is_Svtrack = true;
0349     b_trackid = sv_trkid;
0350     b_trkpt = svtrk->get_pt();
0351     b_trketa = svtrk->get_eta();
0352     b_trkphi = svtrk->get_phi();
0353
0354     cnt = m_TCEval.svtx_cntclus();
0355     b_nclus = cnt[4];
0356     b_nclusmvtx = cnt[0];
0357     b_nclusintt = cnt[1];
0358     b_nclustpc = cnt[2];
0359
0360     b_matchrat = (float) cnt_match[4] / cnt[4];
0361     b_matchrat_mvtx = (float) cnt_match[0] / cnt[0];
0362     b_matchrat_intt = (float) cnt_match[1] / cnt[1];
0363     b_matchrat_tpc = (float) cnt_match[2] / cnt[2];
0364
0365     /* int _ = 0; */
0366     if (m_fill_clusters)
0367     {
0368       auto clusters = m_TCEval.svtx_clusloc_unmatched();
0369       for (auto& clus : clusters)
0370       {
0371         b_clusmatch.push_back(false);
0372         fill_cluster_branches(clus, false);
0373       }
0374
0375       clusters = m_TCEval.clusloc_matched();
0376       for (auto& clus : clusters)
0377       {
0378         b_clusmatch.push_back(true);
0379         fill_cluster_branches(clus, false);
0380       }
0381     }
0382     m_ttree->Fill();
0383     clear_clusvecs();
0384   }
0385
0386   // <- <- <- <- G4 un-matched Tracks
0387   b_is_matched = false;
0388   b_is_g4track = true;
0389   b_is_Svtrack = false;
0390   for (auto& g4_trkid : m_EmbRecoMatchContainer->ids_TruthUnmatched())
0391   {
0392     auto g4trk = m_TrkrTruthTrackContainer->getTruthTrack(g4_trkid);
0393     m_TCEval.addClusKeys(g4trk);
0394
0395     b_trackid = g4_trkid;
0396     b_trkpt = g4trk->getPt();
0397     b_trketa = g4trk->getPseudoRapidity();
0398     b_trkphi = g4trk->getPhi();
0399
0400     auto cnt = m_TCEval.phg4_cntclus();
0401     b_nclus = cnt[4];
0402     b_nclusmvtx = cnt[0];
0403     b_nclusintt = cnt[1];
0404     b_nclustpc = cnt[2];
0405
0406     if (m_fill_clusters)
0407     {
0408       for (auto& clus : m_TCEval.phg4_clusloc_all())
0409       {
0410         b_clusmatch.push_back(false);
0411         fill_cluster_branches(clus, true);
0412       }
0413       // this is an unmatched track, so there are no matched clusters
0414     }
0415     m_ttree->Fill();
0416     clear_clusvecs();
0417   }
0418
0419   // <- <- <- <- Svtx un-matched Tracks
0420   b_is_matched = false;
0421   b_is_matched = false;
0422   b_is_g4track = false;
0423   b_is_Svtrack = true;
0424
0425   // put in all unmatched svtx tracks, too
0426   if (m_fill_SvUnmatched)
0427   {
0428     for (auto sv_trkid : g4evalfn::unmatchedSvtxTrkIds(m_EmbRecoMatchContainer, m_SvtxTrackMap))
0429     {
0430       auto svtrk = m_SvtxTrackMap->get(sv_trkid);
0431       m_TCEval.addClusKeys(svtrk);
0432       b_trackid = sv_trkid;
0433       b_trkpt = svtrk->get_pt();
0434       b_trketa = svtrk->get_eta();
0435       b_trkphi = svtrk->get_phi();
0436
0437       auto cnt = m_TCEval.svtx_cntclus();
0438       b_nclus = cnt[4];
0439       b_nclusmvtx = cnt[0];
0440       b_nclusintt = cnt[1];
0441       b_nclustpc = cnt[2];
0442
0443       if (m_fill_clusters)
0444       {
0445         for (auto& clus : m_TCEval.svtx_clusloc_all())
0446         {
0447           fill_cluster_branches(clus, false);
0448         }
0449       }
0450       m_ttree->Fill();
0451       clear_clusvecs();
0452     }
0453   }
0454
0455   if (Verbosity() > 100)
0456   {
0457     print_mvtx_diagnostics();
0458   }
0459   return Fun4AllReturnCodes::EVENT_OK;
0460 }
0461
0462 void FillClusMatchTree::print_mvtx_diagnostics()
0463 {
0464   std::cout << "To do: "
0465             << " (1)  number of truth tracks and total number of mvtx and ratio " << std::endl
0466             << " (2)  ditto for reco tracks " << std::endl;
0467
0468   double n_PHG4_tracks = m_TrkrTruthTrackContainer->getMap().size();
0469   // count how many mvtx clusters in the phg4 tracks
0470   double n_in_PHG4_tracks{0.};
0471   for (auto& pair : m_TrkrTruthTrackContainer->getMap())
0472   {
0473     m_TCEval.addClusKeys(pair.second);
0474     n_in_PHG4_tracks += m_TCEval.phg4_cntclus()[0];
0475   }
0476   // count how mant mvtx clusters in truth container (should be the same)
0477   double n_in_PHG4_clusters{0.};
0478   for (auto hitsetkey : m_TCEval.get_PHG4_clusters()->getHitSetKeys())
0479   {
0480     if (TrkrDefs::getLayer(hitsetkey) > 2)
0481     {
0482       continue;
0483     }
0484     auto range = m_TCEval.get_PHG4_clusters()->getClusters(hitsetkey);
0485     for (auto r = range.first; r != range.second; ++r)
0486     {
0487       n_in_PHG4_clusters += 1.;
0488     }
0489   }
0490
0491   // count how many svtx tracks
0492   double n_SVTX_tracks = m_SvtxTrackMap->size();
0493   // count how many mvtx clusters in svtx tracks
0494   double n_in_SVTX_tracks{0.};
0495   for (auto& entry : *m_SvtxTrackMap)
0496   {
0497     m_TCEval.addClusKeys(entry.second);
0498     n_in_SVTX_tracks += m_TCEval.svtx_cntclus()[0];
0499   }
0500   // count how many mvtx are total in the container
0501   double n_in_SVTX_clusters{0.};
0502   for (auto hitsetkey : m_TCEval.get_SVTX_clusters()->getHitSetKeys())
0503   {
0504     if (TrkrDefs::getLayer(hitsetkey) > 2)
0505     {
0506       continue;
0507     }
0508     auto range = m_TCEval.get_SVTX_clusters()->getClusters(hitsetkey);
0509     for (auto r = range.first; r != range.second; ++r)
0510     {
0511       n_in_SVTX_clusters += 1.;
0512     }
0513   }
0514
0515   std::cout << (boost::format("MVTX"
0516                               "\nPHG4:  Tracks(%.0f)   Clusters In tracks(%.0f)   Total (%.0f)"
0517                               "\n       ave. per track: %6.3f   ratio in all tracks: %6.2f") %
0518                 n_PHG4_tracks % n_in_PHG4_tracks % n_in_PHG4_clusters % (n_in_PHG4_tracks / n_PHG4_tracks) % (n_in_PHG4_tracks / n_in_PHG4_clusters))
0519                    .str()
0520             << std::endl;
0521   std::cout << (boost::format(
0522                     "\nSVTX:  Tracks(%.0f)   Clusters In tracks(%.0f)   Total (%.0f)"
0523                     "\n       ave. per track: %6.3f   ratio in all tracks: %6.2f") %
0524                 n_SVTX_tracks % n_in_SVTX_tracks % n_in_SVTX_clusters % (n_in_SVTX_tracks / n_SVTX_tracks) % (n_in_SVTX_tracks / n_in_SVTX_clusters))
0525                    .str()
0526             << std::endl;
0527 }
0528
0529 int FillClusMatchTree::End(PHCompositeNode* /*unused*/)
0530 {
0531   if (Verbosity() > 2)
0532   {
0533     std::cout << PHWHERE << ": ending FillClusMatchTree" << std::endl;
0534   }
0535   PHTFileServer::get().cd(m_outfile_name);
0536
0537   h2_G4_nPixelsPhi->Write();
0538   h2_G4_nPixelsZ->Write();
0539   h2_Sv_nPixelsPhi->Write();
0540   h2_Sv_nPixelsZ->Write();
0541
0542   m_ttree->Write();
0543   return Fun4AllReturnCodes::EVENT_OK;
0544 }
0545
0546 void FillClusMatchTree::clear_clusvecs(const std::string& tag)
0547 {
0548   if (tag != "" && b_clus_x.size() > 0)
0549   {
0550     for (auto x : b_clus_x)
0551     {
0552       std::cout << x << " ";
0553     }
0554     std::cout << std::endl;
0555   }
0556
0557   // Tracks and clustes
0558   if (m_fill_clusters)
0559   {
0560     b_clusmatch.clear();
0561     b_clus_x.clear();
0562     b_clus_y.clear();
0563     b_clus_z.clear();
0564     b_clus_r.clear();
0565     b_clus_layer.clear();
0566
0567     b_clus_lphi.clear();
0568     b_clus_lphisize.clear();
0569     b_clus_lz.clear();
0570     b_clus_lzsize.clear();
0571
0572     if (m_fill_clusverbose)
0573     {
0574       b_phibins.clear();
0575       b_phibins_cut.clear();
0576       b_zbins.clear();
0577       b_zbins_cut.clear();
0578
0579       b_phibinsE.clear();
0580       b_phibinsE_cut.clear();
0581       b_zbinsE.clear();
0582       b_zbinsE_cut.clear();
0583     }
0584   }
0585 }
0586
0587 void FillClusMatchTree::fill_cluster_branches(
0588     TrkrClusLoc& clus, bool isPHG4)
0589 {
0590   auto x = clus.gloc[0];  // std::get<1>(loc)[0];
0591   auto y = clus.gloc[1];
0592   auto r = pow(x * x + y * y, 0.5);
0593   b_clus_layer.push_back(clus.layer);
0594   b_clus_x.push_back(x);
0595   b_clus_y.push_back(y);
0596   b_clus_z.push_back(clus.gloc[2]);
0597   b_clus_r.push_back(r);
0598   b_clus_lz.push_back(clus.z);
0599   b_clus_lzsize.push_back(clus.zsize);
0600   b_clus_lphi.push_back(clus.phi);
0601   b_clus_lphisize.push_back(clus.phisize);
0602   if (!m_fill_clusverbose)
0603   {
0604     return;
0605   }
0606
0607   ClusHitsVerbose* data = (isPHG4) ? m_PHG4ClusHitVerb : m_SvtxClusHitVerb;
0608
0609   auto& key = clus.ckey;
0610   if (data->hasClusKey(key))
0611   {
0612     auto phidat = data->phiBins_pvecIE(key);
0613     b_phibins.push_back(phidat.first);
0614     b_phibinsE.push_back(phidat.second);
0615
0616     auto zdat = data->zBins_pvecIE(key);
0617     b_zbins.push_back(zdat.first);
0618     b_zbinsE.push_back(zdat.second);
0619
0620     auto phicut = data->phiCutBins_pvecIE(key);
0621     b_phibins_cut.push_back(phicut.first);
0622     b_phibinsE_cut.push_back(phicut.second);
0623
0624     auto zcut = data->zCutBins_pvecIE(key);
0625     b_zbins_cut.push_back(zcut.first);
0626     b_zbinsE_cut.push_back(zcut.second);
0627   }
0628   else
0629   {
0630     b_phibins.push_back({});
0631     b_phibinsE.push_back({});
0632
0633     b_zbins.push_back({});
0634     b_zbinsE.push_back({});
0635
0636     b_phibins_cut.push_back({});
0637     b_phibinsE_cut.push_back({});
0638
0639     b_zbins_cut.push_back({});
0640     b_zbinsE_cut.push_back({});
0641   }
0642 }