FillClusMatchTree.cc

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file FillClusMatchTree.cc

#include "FillClusMatchTree.h"
#include "TrkrClusterIsMatcher.h"
#include "g4evalfn.h"

#include <TFile.h>
#include <TTree.h>
#include <fun4all/PHTFileServer.h>
#include <fun4all/Fun4AllReturnCodes.h>
#include <g4main/PHG4TruthInfoContainer.h>
#include <g4tracking/EmbRecoMatch.h>
#include <g4tracking/EmbRecoMatchContainer.h>
#include <g4tracking/TrkrTruthTrack.h>
#include <g4tracking/TrkrTruthTrackContainer.h>
#include <iostream>
#include <phool/PHCompositeNode.h>
#include <phool/PHDataNode.h>
#include <phool/PHIODataNode.h>
#include <phool/PHNode.h>
#include <phool/PHNodeIterator.h>
#include <phool/PHObject.h>  // for PHObject
#include <phool/getClass.h>
#include <phool/phool.h>  // for PHWHERE
#include <trackbase/ClusHitsVerbose.h>
#include <trackbase/TrkrDefs.h>
#include <trackbase_historic/SvtxPHG4ParticleMap_v1.h>
#include <trackbase_historic/SvtxTrack.h>
#include <trackbase_historic/SvtxTrackMap.h>
#include <TH2D.h>
#include <trackbase/TrkrCluster.h>

using std::cout;
using std::endl;
//____________________________________________________________________________..
FillClusMatchTree::FillClusMatchTree(
     TrkrClusterIsMatcher* _ismatcher
    , const std::string& _outfile_name
    , bool _fill_clusters
    , bool _fill_clusverbose
    , bool _fill_SvUnMatched
    )
    : m_ismatcher        { _ismatcher        }
    , m_outfile_name     { _outfile_name     }
    , m_fill_clusters    { _fill_clusters    }
    , m_fill_clusverbose { _fill_clusverbose }
    , m_fill_SvUnmatched  { _fill_SvUnMatched }
{
  m_TCEval.ismatcher = m_ismatcher;

  PHTFileServer::get().open(m_outfile_name, "RECREATE");

  h2_G4_nPixelsPhi = new TH2D("G4_nPixelsPhi", "PHG4 Emb Tracks; cluster pixel width Phi; layer",
      100, -0.5, 99.5, 56, -0.5, 55.5);
  h2_G4_nPixelsZ = new TH2D("G4_nPixelsZ", "PHG4 Emb Tracks; cluster pixel width Z; layer",
      100, -0.5, 99.5, 56, -0.5, 55.5);
  h2_Sv_nPixelsPhi = new TH2D("Sv_nPixelsPhi", "Svtx Reco Tracks; cluster pixel width Phi; layer",
      100, -0.5, 99.5, 56, -0.5, 55.5);
  h2_Sv_nPixelsZ = new TH2D("Sv_nPixelsZ", "Svtx Reco Tracks; cluster pixel width Z; layer",
      100, -0.5, 99.5, 56, -0.5, 55.5);

  m_ttree = new TTree("T", "Tracks (and sometimes clusters)");

  m_ttree->Branch("event",                  &nevent);
  m_ttree->Branch("nphg4_part",             &nphg4_part);
  m_ttree->Branch("centrality",             &centrality);
  m_ttree->Branch("ntrackmatches",          &ntrackmatches);
  m_ttree->Branch("nphg4",                  &nphg4);
  m_ttree->Branch("nsvtx",                  &nsvtx);

  m_ttree ->Branch("trackid"    , &b_trackid    );
  m_ttree  ->Branch("is_G4track" , &b_is_g4track );
  m_ttree  ->Branch("is_Svtrack" , &b_is_Svtrack );
  m_ttree  ->Branch("is_matched" , &b_is_matched );

  m_ttree ->Branch("trkpt"  , &b_trkpt  );
  m_ttree ->Branch("trketa" , &b_trketa );
  m_ttree ->Branch("trkphi" , &b_trkphi );

  m_ttree  ->Branch("nclus"         , &b_nclus         );
  m_ttree  ->Branch("nclustpc"      , &b_nclustpc      );
  m_ttree  ->Branch("nclusmvtx"     , &b_nclusmvtx     );
  m_ttree  ->Branch("nclusintt"     , &b_nclusintt     );
  m_ttree  ->Branch("matchrat"      , &b_matchrat      );
  m_ttree  ->Branch("matchrat_intt" , &b_matchrat_intt );
  m_ttree  ->Branch("matchrat_mvtx" , &b_matchrat_mvtx );
  m_ttree  ->Branch("matchrat_tpc"  , &b_matchrat_tpc  );

  if (m_fill_clusters) {
    m_ttree ->Branch("clus_match"   , &b_clusmatch     );
    m_ttree ->Branch("clus_x"       , &b_clus_x        );
    m_ttree ->Branch("clus_y"       , &b_clus_y        );
    m_ttree ->Branch("clus_z"       , &b_clus_z        );
    m_ttree ->Branch("clus_r"       , &b_clus_r        );
    m_ttree ->Branch("clus_layer"   , &b_clus_layer    );
    /* m_ttree ->Branch("nphibins"     , &b_clus_nphibins ); */
    /* m_ttree ->Branch("nzbins"       , &b_clus_ntbins   ); */

    m_ttree ->Branch("clus_lphi"     , &b_clus_lphi     );
    m_ttree ->Branch("clus_lphisize" , &b_clus_lphisize );
    m_ttree ->Branch("clus_lz"       , &b_clus_lz       );
    m_ttree ->Branch("clus_lzsize"   , &b_clus_lzsize   );
    
    if (m_fill_clusverbose) {
      m_ttree ->Branch("phibins"     , &b_phibins     );
      m_ttree ->Branch("phibins_cut" , &b_phibins_cut );
      m_ttree ->Branch("zbins"       , &b_zbins     );
      m_ttree ->Branch("zbins_cut"   , &b_zbins_cut );

      m_ttree ->Branch("phibinsE"     , &b_phibinsE     );
      m_ttree ->Branch("phibinsE_cut" , &b_phibinsE_cut );
      m_ttree ->Branch("zbinsE"       , &b_zbinsE     );
      m_ttree ->Branch("zbinsE_cut"   , &b_zbinsE_cut );
    }
  }
}

//____________________________________________________________________________..
FillClusMatchTree::~FillClusMatchTree()
{
}

//____________________________________________________________________________..
int FillClusMatchTree::Init(PHCompositeNode *topNode)
{
  if (Verbosity()>1) {
    std::cout << " Beginning FillClusMatchTree " << std::endl;
    topNode->print();
  }

  return Fun4AllReturnCodes::EVENT_OK;
}

//____________________________________________________________________________..
int FillClusMatchTree::InitRun(PHCompositeNode *topNode)
{
  /* auto init_status = m_cluster_comp.init(topNode); */
  auto init_status = m_ismatcher->init(topNode);
  if (init_status == Fun4AllReturnCodes::ABORTRUN) return init_status;

  if (createNodes(topNode) != Fun4AllReturnCodes::EVENT_OK)
  {
    return Fun4AllReturnCodes::ABORTEVENT;
  }

  return Fun4AllReturnCodes::EVENT_OK;
}

int FillClusMatchTree::createNodes(PHCompositeNode *topNode)
{
  PHNodeIterator iter(topNode);

  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));

  if (!dstNode)
  {
    std::cout << PHWHERE << " DST node is missing, quitting" << std::endl;
    std::cerr << PHWHERE << " DST node is missing, quitting" << std::endl;
    throw std::runtime_error("Failed to find DST node in FillClusMatchTree::createNodes");
  }
  topNode->print();

  m_PHG4ClusHitVerb = findNode::getClass<ClusHitsVerbose>(topNode,"Trkr_TruthClusHitsVerbose");
  if (m_PHG4ClusHitVerb) {
    std::cout << " Found truth cluster hits verbose " << std::endl;
  } else {
    std::cout << " Did not find truth cluster hits verbose " << std::endl;
  }

  m_SvtxClusHitVerb = findNode::getClass<ClusHitsVerbose>(topNode,"Trkr_SvtxClusHitsVerbose");
  if (m_SvtxClusHitVerb) {
    std::cout << " Found Svtx cluster hits verbose " << std::endl;
  } else {
    std::cout << " Did not find Svtx cluster hits verbose " << std::endl;
  }

  m_EmbRecoMatchContainer = findNode::getClass<EmbRecoMatchContainer>(topNode,"TRKR_EMBRECOMATCHCONTAINER");
  if (!m_EmbRecoMatchContainer) {
    std::cout << PHWHERE << " Cannot find node TRKR_EMBRECOMATCHCONTAINER on node tree; quitting " << std::endl;
    std::cerr << PHWHERE << " Cannot find node TRKR_EMBRECOMATCHCONTAINER on node tree; quitting " << std::endl;
    throw std::runtime_error(" Cannot find node TRKR_EMBRECOMATCHCONTAINER on node tree; quitting");
  }

  PHCompositeNode *svtxNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "SVTX"));
  if (!svtxNode)
  {
    svtxNode = new PHCompositeNode("SVTX");
    dstNode->addNode(svtxNode);
  }

  m_PHG4TruthInfoContainer = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
  if (!m_PHG4TruthInfoContainer)
  {
    std::cout << "Could not locate G4TruthInfo node when running " 
      << "\"TruthRecoTrackMatching\" module." << std::endl;
    return Fun4AllReturnCodes::ABORTEVENT;
  }

  m_SvtxTrackMap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
  if (!m_SvtxTrackMap)
  {
    std::cout << "Could not locate SvtxTrackMap node when running "
      <<"\"TruthRecoTrackMatching\" module." << std::endl;
    return Fun4AllReturnCodes::ABORTEVENT;
  }

  m_TrkrTruthTrackContainer = findNode::getClass<TrkrTruthTrackContainer>(topNode, 
      "TRKR_TRUTHTRACKCONTAINER");
  if (!m_TrkrTruthTrackContainer)
  {
    std::cout << "Could not locate TRKR_TRUTHTRACKCONTAINER node when running "
      << "\"TruthRecoTrackMatching\" module." << std::endl;
    return Fun4AllReturnCodes::ABORTEVENT;
  }

  return Fun4AllReturnCodes::EVENT_OK;
}

int FillClusMatchTree::process_event(PHCompositeNode * /*topNode*/)
{
  // We know that the data is in there and exists... why isn't it being transferred?

  if (Verbosity()>5) cout << " FillClusMatchTree::process_event() " << endl;

  // fill in the event data
  ++nevent; 
  nphg4 = m_TrkrTruthTrackContainer ->getMap().size();
  nsvtx = m_SvtxTrackMap            ->size();
  ntrackmatches = m_EmbRecoMatchContainer->getMatches().size();
  // get centrality later...

  // fill in pixel widths on truth tracks
  for (auto hitsetkey : m_TCEval.get_PHG4_clusters()->getHitSetKeys()) {
    float layer = (float) TrkrDefs::getLayer(hitsetkey);
    auto range = m_TCEval.get_PHG4_clusters()->getClusters(hitsetkey);
    for (auto iter = range.first; iter != range.second; ++iter) {
      auto& cluster = iter->second;
      h2_G4_nPixelsPhi ->Fill( (float)cluster->getPhiSize(), layer );
      h2_G4_nPixelsZ   ->Fill( (float)cluster->getZSize(),   layer );
    }
  }
  // fill in pixel widths on reco tracks
  for (auto hitsetkey : m_TCEval.get_SVTX_clusters()->getHitSetKeys()) {
    float layer = (float) TrkrDefs::getLayer(hitsetkey);
    auto range = m_TCEval.get_SVTX_clusters()->getClusters(hitsetkey);
    for (auto iter = range.first; iter != range.second; ++iter) {
      auto& cluster = iter->second;
      h2_Sv_nPixelsPhi ->Fill( (float)cluster->getPhiSize(), layer );
      h2_Sv_nPixelsZ   ->Fill( (float)cluster->getZSize(), layer );
    }
  }

  nphg4_part = 0;
  const auto range = m_PHG4TruthInfoContainer->GetPrimaryParticleRange();
  for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
      iter != range.second; ++iter)
  { nphg4_part++; }

  // unmatched tracks are only entered once
  // matches can repeat a given svtx or phg4 track, depending on the 
  // parameters in teh matching in FillClusMatchTree
  //
  // (1) fill unmatched phg4
  // (2) fill unmatched svtx
  // (3) fill matched phg4 and svtx
  clear_clusvecs(" nothing "); 

  if (Verbosity() > 2) std::cout << " getting " << (int) m_EmbRecoMatchContainer->getMatches().size() << std::endl;
  for (auto& match : m_EmbRecoMatchContainer->getMatches()) {
    unsigned int g4_trkid = match->idTruthTrack();
    int sv_trkid = match->idRecoTrack();

    auto g4trk = m_TrkrTruthTrackContainer->getTruthTrack(g4_trkid);
    auto svtrk = m_SvtxTrackMap->get(sv_trkid);

    m_TCEval.addClusKeys( g4trk );
    m_TCEval.addClusKeys( svtrk );
    m_TCEval.find_matches();


    // <- <- <- <- G4 Matched Tracks
    b_is_matched = true;
    b_is_g4track = true;
    b_is_Svtrack = false;

    b_trackid    = g4_trkid;
    b_trkpt      = g4trk->getPt();
    b_trketa     = g4trk->getPseudoRapidity();
    b_trkphi     = g4trk->getPhi();

    auto cnt = m_TCEval.phg4_cntclus();
    auto cnt_match = m_TCEval.phg4_cnt_matchedclus();

    b_nclus     = cnt[4] ;
    b_nclusmvtx = cnt[0] ;
    b_nclusintt = cnt[1] ;
    b_nclustpc  = cnt[2] ;

    b_matchrat      = (float)cnt_match[4] / cnt[4] ;
    b_matchrat_mvtx = (float)cnt_match[0] / cnt[0] ;
    b_matchrat_intt = (float)cnt_match[1] / cnt[1] ;
    b_matchrat_tpc  = (float)cnt_match[2] / cnt[2] ;

    if (m_fill_clusters) {
      // - - - - UNMATCHED PHG4 hits in matched track
      auto clusters = m_TCEval.phg4_clusloc_unmatched();
      for (auto& clus : clusters) {
          b_clusmatch     .push_back(false);
          fill_cluster_branches(clus, true);
      }

      // - - - - MATCHED PHG4 hits in matched track
      clusters = m_TCEval.clusloc_matched();
      for (auto& clus : clusters) {
        b_clusmatch     .push_back(true);
        fill_cluster_branches(clus, true);
      }
    }
    m_ttree->Fill();
    clear_clusvecs();
    /* clear_clusvecs("apple0 g4_matched"); */

    // <- <- <- <- Svtx Matched Tracks
    b_is_g4track = false;
    b_is_Svtrack = true;
    b_trackid    = sv_trkid;
    b_trkpt      = svtrk->get_pt();
    b_trketa     = svtrk->get_eta();
    b_trkphi     = svtrk->get_phi();

    cnt = m_TCEval.svtx_cntclus();
    b_nclus     = cnt[4] ;
    b_nclusmvtx = cnt[0] ;
    b_nclusintt = cnt[1] ;
    b_nclustpc  = cnt[2] ;

    b_matchrat      = (float)cnt_match[4] / cnt[4] ;
    b_matchrat_mvtx = (float)cnt_match[0] / cnt[0] ;
    b_matchrat_intt = (float)cnt_match[1] / cnt[1] ;
    b_matchrat_tpc  = (float)cnt_match[2] / cnt[2] ;

    /* int _ = 0; */
    if (m_fill_clusters) {
      auto clusters = m_TCEval.svtx_clusloc_unmatched();
      for (auto& clus : clusters) {
        b_clusmatch     .push_back(false);
        fill_cluster_branches(clus, false);
      }

      clusters = m_TCEval.clusloc_matched();
      for (auto& clus : clusters) {
        b_clusmatch     .push_back(true);
        fill_cluster_branches(clus, false);
      }
    }
    m_ttree->Fill();
    clear_clusvecs();
  }

  // <- <- <- <- G4 un-matched Tracks
  b_is_matched = false;
  b_is_g4track = true;
  b_is_Svtrack = false;
  for (auto& g4_trkid : m_EmbRecoMatchContainer->ids_TruthUnmatched()) {
    auto g4trk = m_TrkrTruthTrackContainer->getTruthTrack(g4_trkid);
    m_TCEval.addClusKeys( g4trk );

    b_trackid  = g4_trkid;
    b_trkpt    = g4trk->getPt();
    b_trketa   = g4trk->getPseudoRapidity();
    b_trkphi   = g4trk->getPhi();

    auto cnt = m_TCEval.phg4_cntclus();
    b_nclus     = cnt[4];
    b_nclusmvtx = cnt[0];
    b_nclusintt = cnt[1];
    b_nclustpc  = cnt[2];

    if (m_fill_clusters) {
      for (auto& clus : m_TCEval.phg4_clusloc_all()) {
        b_clusmatch     .push_back(false);
        fill_cluster_branches(clus, true);
      }
      //this is an unmatched track, so there are no matched clusters
    }
    m_ttree->Fill();
    clear_clusvecs();
  }

  // <- <- <- <- Svtx un-matched Tracks
  b_is_matched = false;
  b_is_matched = false;
  b_is_g4track = false;
  b_is_Svtrack = true;

  // put in all unmatched svtx tracks, too
  if (m_fill_SvUnmatched) {
    for (auto sv_trkid : g4evalfn::unmatchedSvtxTrkIds(m_EmbRecoMatchContainer, m_SvtxTrackMap)) {
      auto svtrk = m_SvtxTrackMap->get(sv_trkid);
      m_TCEval.addClusKeys( svtrk );
      b_trackid  = sv_trkid;
      b_trkpt    = svtrk->get_pt();
      b_trketa   = svtrk->get_eta();
      b_trkphi   = svtrk->get_phi();

      auto cnt = m_TCEval.svtx_cntclus();
      b_nclus     = cnt[4] ;
      b_nclusmvtx = cnt[0] ;
      b_nclusintt = cnt[1] ;
      b_nclustpc  = cnt[2] ;

      if (m_fill_clusters) {
        for (auto &clus : m_TCEval.svtx_clusloc_all()) {
            fill_cluster_branches(clus, false);
          }
      }
      m_ttree->Fill();
      clear_clusvecs();
    }
  }

  if (Verbosity()>100) print_mvtx_diagnostics();
  return Fun4AllReturnCodes::EVENT_OK;
}

void FillClusMatchTree::print_mvtx_diagnostics() {
  std::cout << "To do: "
    << " (1)  number of truth tracks and total number of mvtx and ratio " << std::endl
    << " (2)  ditto for reco tracks " << std::endl;

  double n_PHG4_tracks = m_TrkrTruthTrackContainer->getMap().size();
  // count how many mvtx clusters in the phg4 tracks
  double n_in_PHG4_tracks {0.};
  for (auto& pair : m_TrkrTruthTrackContainer->getMap()) {
    m_TCEval.addClusKeys( pair.second );
    n_in_PHG4_tracks += m_TCEval.phg4_cntclus()[0];
  }
  // count how mant mvtx clusters in truth container (should be the same)
  double n_in_PHG4_clusters {0.};
  for (auto hitsetkey : m_TCEval.get_PHG4_clusters()->getHitSetKeys()) {
    if (TrkrDefs::getLayer(hitsetkey) > 2) continue;
    auto range = m_TCEval.get_PHG4_clusters()->getClusters(hitsetkey);
    for (auto r = range.first; r != range.second; ++r)  {
      n_in_PHG4_clusters += 1.;
    }
  }

  // count how many svtx tracks
  double n_SVTX_tracks = m_SvtxTrackMap->size();
  // count how many mvtx clusters in svtx tracks
  double n_in_SVTX_tracks {0.};
  for (auto entry = m_SvtxTrackMap->begin(); entry != m_SvtxTrackMap->end(); ++entry) {
    m_TCEval.addClusKeys(entry->second);
    n_in_SVTX_tracks += m_TCEval.svtx_cntclus()[0];
  }
  // count how many mvtx are total in the container
  double n_in_SVTX_clusters {0.};
  for (auto hitsetkey : m_TCEval.get_SVTX_clusters()->getHitSetKeys()) {
    if (TrkrDefs::getLayer(hitsetkey) > 2) continue;
    auto range = m_TCEval.get_SVTX_clusters()->getClusters(hitsetkey);
    for (auto r = range.first; r != range.second; ++r)  {
      n_in_SVTX_clusters += 1.;
    }
  }

  std::cout << Form("MVTX"
      "\nPHG4:  Tracks(%.0f)   Clusters In tracks(%.0f)   Total (%.0f)"
      "\n       ave. per track: %6.3f   ratio in all tracks: %6.2f"
      , 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))
    << std::endl;
  std::cout << Form(
      "\nSVTX:  Tracks(%.0f)   Clusters In tracks(%.0f)   Total (%.0f)"
      "\n       ave. per track: %6.3f   ratio in all tracks: %6.2f"
      , 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))
    << std::endl;

}

  int FillClusMatchTree::End(PHCompositeNode *)
  {
    if (Verbosity()>2) std::cout << PHWHERE << ": ending FillClusMatchTree" << std::endl;
    PHTFileServer::get().cd(m_outfile_name);

    h2_G4_nPixelsPhi ->Write();
    h2_G4_nPixelsZ   ->Write();
    h2_Sv_nPixelsPhi ->Write();
    h2_Sv_nPixelsZ   ->Write();

    m_ttree->Write();
    return Fun4AllReturnCodes::EVENT_OK;
  }

void FillClusMatchTree::clear_clusvecs(const std::string& tag) {
  if (tag != "" && b_clus_x.size()>0) {
    for (auto x : b_clus_x) cout << x <<" ";
    cout << endl;
  }

  // Tracks and clustes
  if (m_fill_clusters) {
    b_clusmatch     .clear();
    b_clus_x        .clear();
    b_clus_y        .clear();
    b_clus_z        .clear();
    b_clus_r        .clear();
    b_clus_layer    .clear();

    b_clus_lphi      .clear();
    b_clus_lphisize  .clear();
    b_clus_lz        .clear();
    b_clus_lzsize    .clear();

    if (m_fill_clusverbose) {
      b_phibins     .clear();
      b_phibins_cut .clear();
      b_zbins       .clear();
      b_zbins_cut   .clear();

      b_phibinsE     .clear();
      b_phibinsE_cut .clear();
      b_zbinsE       .clear();
      b_zbinsE_cut   .clear();
    }
  }
}

void FillClusMatchTree::fill_cluster_branches(
    TrkrClusLoc& clus, bool isPHG4) 
{
  auto x = clus.gloc[0];//std::get<1>(loc)[0];
  auto y = clus.gloc[1];
  auto r = pow(x*x+y*y,0.5);
  b_clus_layer    .push_back(clus.layer);
  b_clus_x        .push_back(x);
  b_clus_y        .push_back(y);
  b_clus_z        .push_back(clus.gloc[2]);
  b_clus_r        .push_back(r);
  b_clus_lz       .push_back(clus.z);
  b_clus_lzsize   .push_back(clus.zsize);
  b_clus_lphi     .push_back(clus.phi);
  b_clus_lphisize .push_back(clus.phisize);
  if (!m_fill_clusverbose) return;

  ClusHitsVerbose* data = (isPHG4) ? m_PHG4ClusHitVerb : m_SvtxClusHitVerb;

  auto& key = clus.ckey;
  if (data->hasClusKey(key)) {
    auto phidat = data->phiBins_pvecIE(key);
    b_phibins  .push_back(phidat.first);
    b_phibinsE .push_back(phidat.second);

    auto zdat = data->zBins_pvecIE(key);
    b_zbins  .push_back(zdat.first);
    b_zbinsE .push_back(zdat.second);

    auto phicut = data->phiCutBins_pvecIE(key);
    b_phibins_cut  .push_back(phicut.first);
    b_phibinsE_cut .push_back(phicut.second);

    auto zcut = data->zCutBins_pvecIE(key);
    b_zbins_cut  .push_back(zcut.first);
    b_zbinsE_cut .push_back(zcut.second);
  } else {
    b_phibins  .push_back({});
    b_phibinsE .push_back({});

    b_zbins  .push_back({});
    b_zbinsE .push_back({});

    b_phibins_cut  .push_back({});
    b_phibinsE_cut .push_back({});

    b_zbins_cut  .push_back({});
    b_zbinsE_cut .push_back({});
  }
}