.. _program_listing_file_Root_FatJetContainer.cxx:

Program Listing for File FatJetContainer.cxx
============================================

|exhale_lsh| :ref:`Return to documentation for file <file_Root_FatJetContainer.cxx>` (``Root/FatJetContainer.cxx``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   #include "xAODAnaHelpers/FatJetContainer.h"
   #include <xAODAnaHelpers/HelperFunctions.h>
   #include <iostream>
   #include "xAODTruth/TruthEventContainer.h"
   
   using namespace xAH;
   
   FatJetContainer::FatJetContainer(const std::string& name, const std::string& detailStr, const std::string& subjetDetailStr, const std::string& suffix,
                    float units, bool mc)
     : ParticleContainer(name,detailStr,units,mc, true, true, suffix)
   {
     if (m_infoSwitch.m_scales) {
         m_JetConstitScaleMomentum_eta       = new std::vector<float>();
         m_JetConstitScaleMomentum_phi       = new std::vector<float>();
         m_JetConstitScaleMomentum_m       = new std::vector<float>();
         m_JetConstitScaleMomentum_pt        = new std::vector<float>();
   
         m_JetEMScaleMomentum_eta        = new std::vector<float>();
         m_JetEMScaleMomentum_phi        = new std::vector<float>();
         m_JetEMScaleMomentum_m        = new std::vector<float>();
         m_JetEMScaleMomentum_pt       = new std::vector<float>();
     }
   
     if (m_infoSwitch.m_area) {
       m_GhostArea = new std::vector<float>();
       m_ActiveArea = new std::vector<float>();
       m_VoronoiArea = new std::vector<float>();
   
       m_ActiveArea4vec_pt = new std::vector<float>();
       m_ActiveArea4vec_eta = new std::vector<float>();
       m_ActiveArea4vec_phi = new std::vector<float>();
       m_ActiveArea4vec_m = new std::vector<float>();
     }
   
     if ( m_infoSwitch.m_substructure ) {
       m_Split12           = new std::vector<float>();
       m_Split23           = new std::vector<float>();
       m_Split34           = new std::vector<float>();
       m_tau1_wta          = new std::vector<float>();
       m_tau2_wta          = new std::vector<float>();
       m_tau3_wta          = new std::vector<float>();
       m_tau21_wta         = new std::vector<float>();
       m_tau32_wta         = new std::vector<float>();
       m_ECF1              = new std::vector<float>();
       m_ECF2              = new std::vector<float>();
       m_ECF3              = new std::vector<float>();
       m_C2                = new std::vector<float>();
       m_D2                = new std::vector<float>();
       m_NTrimSubjets      = new std::vector<float>();
       m_NClusters         = new std::vector<int>  ();
       m_nTracks           = new std::vector<int>  ();
       m_ungrtrk500    = new std::vector<int>  ();
       m_EMFrac        = new std::vector<float>();
       m_nChargedParticles = new std::vector<int>();
     }
   
     if ( m_infoSwitch.m_ntrimsubjets && !m_infoSwitch.m_substructure ) {
       m_NTrimSubjets = new std::vector<float>();
     }
   
     if ( m_infoSwitch.m_constituent) {
       m_numConstituents    = new std::vector< int >();
     }
   
     if ( m_infoSwitch.m_constituentAll) {
       m_constituentWeights  = new std::vector< std::vector<float> >();
       m_constituent_pt      = new std::vector< std::vector<float> >();
       m_constituent_eta     = new std::vector< std::vector<float> >();
       m_constituent_phi     = new std::vector< std::vector<float> >();
       m_constituent_e       = new std::vector< std::vector<float> >();
     }
   
     
     if ( m_infoSwitch.m_truth && m_mc ) {
       m_truth_m  =new std::vector<float>;
       m_truth_pt =new std::vector<float>;
       m_truth_phi=new std::vector<float>;
       m_truth_eta=new std::vector<float>;
     }
   
     if ( m_infoSwitch.m_bosonCount && m_mc) {
       m_nTQuarks  = new std::vector< int > ();
       m_nHBosons  = new std::vector< int > ();
       m_nWBosons  = new std::vector< int > ();
       m_nZBosons  = new std::vector< int > ();
     }
   
     if (m_infoSwitch.m_muonCorrection) {
       m_muonCorrected_pt  = new std::vector<float>();
       m_muonCorrected_eta = new std::vector<float>();
       m_muonCorrected_phi = new std::vector<float>();
       m_muonCorrected_m   = new std::vector<float>();
     }
   
     for(const auto& trackJetName : m_infoSwitch.m_trackJetNames)
       {
         std::string trkJetName = name;
         if( !suffix.empty() ){ trkJetName += "_"+suffix; }
         trkJetName += "_"+trackJetName;
         m_trkJets[trackJetName] = new xAH::JetContainer(trkJetName, subjetDetailStr, m_units, m_mc);
   
         m_trkJetsIdx[trackJetName] = new std::vector<std::vector<unsigned int> > ();
       }
   
   }
   
   FatJetContainer::~FatJetContainer()
   {
     if(m_debug) std::cout << " Deleting FatJetContainer "  << std::endl;
   
     if ( m_infoSwitch.m_scales ) {
         delete m_JetConstitScaleMomentum_eta ;
         delete m_JetConstitScaleMomentum_phi ;
         delete m_JetConstitScaleMomentum_m   ;
         delete m_JetConstitScaleMomentum_pt  ;
   
         delete m_JetEMScaleMomentum_eta      ;
         delete m_JetEMScaleMomentum_phi      ;
         delete m_JetEMScaleMomentum_m        ;
         delete m_JetEMScaleMomentum_pt       ;
     }
   
     if ( m_infoSwitch.m_area ) {
       delete m_GhostArea;
       delete m_ActiveArea;
       delete m_VoronoiArea;
   
       delete m_ActiveArea4vec_pt;
       delete m_ActiveArea4vec_eta;
       delete m_ActiveArea4vec_phi;
       delete m_ActiveArea4vec_m;
     }
   
     if ( m_infoSwitch.m_substructure ) {
       delete m_Split12     ;
       delete m_Split23     ;
       delete m_Split34     ;
       delete m_tau1_wta    ;
       delete m_tau2_wta    ;
       delete m_tau3_wta    ;
       delete m_tau21_wta   ;
       delete m_tau32_wta   ;
       delete m_ECF1        ;
       delete m_ECF2        ;
       delete m_ECF3        ;
       delete m_C2          ;
       delete m_D2          ;
       delete m_NTrimSubjets;
       delete m_NClusters   ;
       delete m_nTracks   ;
       delete m_ungrtrk500     ;
       delete m_EMFrac     ;
       delete m_nChargedParticles  ;
     }
   
     if ( m_infoSwitch.m_ntrimsubjets && !m_infoSwitch.m_substructure ){
       delete m_NTrimSubjets;
     }
   
     if ( m_infoSwitch.m_constituent) {
       delete m_numConstituents;
     }
   
     if ( m_infoSwitch.m_constituentAll) {
       delete m_constituentWeights;
       delete m_constituent_pt    ;
       delete m_constituent_eta   ;
       delete m_constituent_phi   ;
       delete m_constituent_e     ;
     }
   
     if ( m_infoSwitch.m_truth && m_mc ) {
       delete m_truth_m;
       delete m_truth_pt;
       delete m_truth_phi;
       delete m_truth_eta;
     }
     
     if ( m_infoSwitch.m_bosonCount && m_mc) {
       delete m_nTQuarks;
       delete m_nHBosons;
       delete m_nWBosons;
       delete m_nZBosons;
     }
     
     if ( m_infoSwitch.m_muonCorrection) {
       delete m_muonCorrected_pt;
       delete m_muonCorrected_eta;
       delete m_muonCorrected_phi;
       delete m_muonCorrected_m;
     }
   
     if( !m_infoSwitch.m_trackJetNames.empty() ){
       for(const auto& kv : m_trkJets)
         {
       delete m_trkJets   [kv.first];
       delete m_trkJetsIdx[kv.first];
         }
       m_trkJets   .clear();
       m_trkJetsIdx.clear();
     }
   
   }
   
   void FatJetContainer::setTree(TTree *tree)
   {
     //
     // Connect branches
     ParticleContainer::setTree(tree);
   
     if( m_infoSwitch.m_scales ) {
       connectBranch<float>(tree, "JetConstitScaleMomentum_eta", &m_JetConstitScaleMomentum_eta);
       connectBranch<float>(tree, "JetConstitScaleMomentum_phi", &m_JetConstitScaleMomentum_phi);
       connectBranch<float>(tree, "JetConstitScaleMomentum_m", &m_JetConstitScaleMomentum_m);
       connectBranch<float>(tree, "JetConstitScaleMomentum_pt", &m_JetConstitScaleMomentum_pt);
   
       connectBranch<float>(tree, "JetEMScaleMomentum_eta", &m_JetEMScaleMomentum_eta);
       connectBranch<float>(tree, "JetEMScaleMomentum_phi", &m_JetEMScaleMomentum_phi);
       connectBranch<float>(tree, "JetEMScaleMomentum_m", &m_JetEMScaleMomentum_m);
       connectBranch<float>(tree, "JetEMScaleMomentum_pt", &m_JetEMScaleMomentum_pt);
     }
   
     if ( m_infoSwitch.m_area ) {
       connectBranch<float>(tree, "m_GhostArea", &m_GhostArea);
       connectBranch<float>(tree, "m_ActiveArea", &m_ActiveArea);
       connectBranch<float>(tree, "m_VoronoiArea", &m_VoronoiArea);
   
       connectBranch<float>(tree, "m_ActiveArea4vec_pt", &m_ActiveArea4vec_pt);
       connectBranch<float>(tree, "m_ActiveArea4vec_eta", &m_ActiveArea4vec_eta);
       connectBranch<float>(tree, "m_ActiveArea4vec_phi", &m_ActiveArea4vec_phi);
       connectBranch<float>(tree, "m_ActiveArea4vec_m", &m_ActiveArea4vec_m);
     }
   
     if ( m_infoSwitch.m_substructure ) {
       connectBranch<float>(tree, "Split12",      &m_Split12);
       connectBranch<float>(tree, "Split23",      &m_Split23);
       connectBranch<float>(tree, "Split34",      &m_Split34);
       connectBranch<float>(tree, "tau1_wta",     &m_tau1_wta);
       connectBranch<float>(tree, "tau2_wta",     &m_tau2_wta);
       connectBranch<float>(tree, "tau3_wta",     &m_tau3_wta);
       connectBranch<float>(tree, "tau21_wta",    &m_tau21_wta);
       connectBranch<float>(tree, "tau32_wta",    &m_tau32_wta);
       connectBranch<float>(tree, "ECF1",         &m_ECF1);
       connectBranch<float>(tree, "ECF2",         &m_ECF2);
       connectBranch<float>(tree, "ECF3",         &m_ECF3);
       connectBranch<float>(tree, "C2",           &m_C2);
       connectBranch<float>(tree, "D2",           &m_D2);
       connectBranch<float>(tree, "NTrimSubjets", &m_NTrimSubjets);
       connectBranch<int>  (tree, "Nclusters",    &m_NClusters);
       connectBranch<int>  (tree, "nTracks",      &m_nTracks);
       connectBranch<int>  (tree, "ungrtrk500",        &m_ungrtrk500);
       connectBranch<float>(tree, "EMFrac",        &m_EMFrac);
       connectBranch<int>  (tree, "nChargedParticles", &m_nChargedParticles);
     }
   
     if ( m_infoSwitch.m_ntrimsubjets && !m_infoSwitch.m_substructure ){
       connectBranch<float>(tree, "NTrimSubjets", &m_NTrimSubjets);
     }
   
     if ( m_infoSwitch.m_constituent) {
       connectBranch<int>  (tree, "numConstituents",    &m_numConstituents);
     }
   
     if ( m_infoSwitch.m_constituentAll) {
       connectBranch< std::vector<float> >(tree, "constituentWeights",  &m_constituentWeights);
       connectBranch< std::vector<float> >(tree, "constituent_pt",      &m_constituent_pt);
       connectBranch< std::vector<float> >(tree, "constituent_eta",     &m_constituent_eta);
       connectBranch< std::vector<float> >(tree, "constituent_phi",     &m_constituent_phi);
       connectBranch< std::vector<float> >(tree, "constituent_e",       &m_constituent_e);
     }
   
     if(m_infoSwitch.m_truth)
       {
         connectBranch<float>(tree,"truth_m",   &m_truth_m);
         connectBranch<float>(tree,"truth_pt",  &m_truth_pt);
         connectBranch<float>(tree,"truth_phi", &m_truth_phi);
         connectBranch<float>(tree,"truth_eta", &m_truth_eta);
       }
   
     if ( m_infoSwitch.m_bosonCount) {
       connectBranch< int >(tree, "nTQuarks",  &m_nTQuarks);
       connectBranch< int >(tree, "nHBosons",  &m_nHBosons);
       connectBranch< int >(tree, "nWBosons",  &m_nWBosons);
       connectBranch< int >(tree, "nZBosons",  &m_nZBosons);
     }
   
     if (m_infoSwitch.m_muonCorrection) {
       connectBranch< float >(tree, "muonCorrected_pt" , &m_muonCorrected_pt );
       connectBranch< float >(tree, "muonCorrected_eta", &m_muonCorrected_eta);
       connectBranch< float >(tree, "muonCorrected_phi", &m_muonCorrected_phi);
       connectBranch< float >(tree, "muonCorrected_m"  , &m_muonCorrected_m  );
       
     } 
   
     for(const std::pair< std::string, std::vector<std::vector<unsigned int>>* >& kv : m_trkJetsIdx)
       {
         m_trkJets[kv.first]->JetContainer::setTree(tree);
         if(tree->GetBranch(branchName("trkJetsIdx").c_str()))
       connectBranch< std::vector<unsigned int> >(tree, "trkJetsIdx", &m_trkJetsIdx[kv.first]);
         else
       connectBranch< std::vector<unsigned int> >(tree, "trkJetsIdx_"+kv.first, &m_trkJetsIdx[kv.first]);
       }
   }
   
   void FatJetContainer::updateParticle(uint idx, FatJet& fatjet)
   {
     if(m_debug) std::cout << "in FatJetContainer::updateParticle " << std::endl;
     ParticleContainer::updateParticle(idx,fatjet);
   
     if ( m_infoSwitch.m_scales ) {
         fatjet.JetConstitScaleMomentum_eta = m_JetConstitScaleMomentum_eta ->at(idx);
         fatjet.JetConstitScaleMomentum_phi = m_JetConstitScaleMomentum_phi ->at(idx);
         fatjet.JetConstitScaleMomentum_m = m_JetConstitScaleMomentum_m ->at(idx);
         fatjet.JetConstitScaleMomentum_pt = m_JetConstitScaleMomentum_pt ->at(idx);
   
         fatjet.JetEMScaleMomentum_eta = m_JetEMScaleMomentum_eta ->at(idx);
         fatjet.JetEMScaleMomentum_phi = m_JetEMScaleMomentum_phi ->at(idx);
         fatjet.JetEMScaleMomentum_m = m_JetEMScaleMomentum_m ->at(idx);
         fatjet.JetEMScaleMomentum_pt = m_JetEMScaleMomentum_pt ->at(idx);
     }
   
     if ( m_infoSwitch.m_area ) {
       fatjet.GhostArea = m_GhostArea->at(idx);
       fatjet.ActiveArea = m_ActiveArea->at(idx);
       fatjet.VoronoiArea = m_VoronoiArea->at(idx);
       fatjet.ActiveArea4vec_pt = m_ActiveArea4vec_pt->at(idx);
       fatjet.ActiveArea4vec_eta = m_ActiveArea4vec_eta->at(idx);
       fatjet.ActiveArea4vec_phi = m_ActiveArea4vec_phi->at(idx);
       fatjet.ActiveArea4vec_m = m_ActiveArea4vec_m->at(idx);
     }
   
     if ( m_infoSwitch.m_substructure ) {
       fatjet.Split12      = m_Split12     ->at(idx);
       fatjet.Split23      = m_Split23     ->at(idx);
       fatjet.Split34      = m_Split34     ->at(idx);
       fatjet.tau1_wta     = m_tau1_wta    ->at(idx);
       fatjet.tau2_wta     = m_tau2_wta    ->at(idx);
       fatjet.tau3_wta     = m_tau3_wta    ->at(idx);
       fatjet.tau21_wta    = m_tau21_wta   ->at(idx);
       fatjet.tau32_wta    = m_tau32_wta   ->at(idx);
       fatjet.ECF1         = m_ECF1        ->at(idx);
       fatjet.ECF2         = m_ECF2        ->at(idx);
       fatjet.ECF3         = m_ECF3        ->at(idx);
       fatjet.C2           = m_C2          ->at(idx);
       fatjet.D2           = m_D2          ->at(idx);
       fatjet.NTrimSubjets = m_NTrimSubjets->at(idx);
       fatjet.NClusters    = m_NClusters   ->at(idx);
       fatjet.nTracks      = m_nTracks     ->at(idx);
       fatjet.ungrtrk500       = m_ungrtrk500      ->at(idx);
       fatjet.EMFrac       = m_EMFrac      ->at(idx);
       fatjet.nChargedParticles    = m_nChargedParticles   ->at(idx);
     }
   
     if ( m_infoSwitch.m_ntrimsubjets && !m_infoSwitch.m_substructure ){
       fatjet.NTrimSubjets = m_NTrimSubjets->at(idx);
     }
   
     if ( m_infoSwitch.m_constituent) {
       fatjet.numConstituents    = m_numConstituents   ->at(idx);
     }
   
     if ( m_infoSwitch.m_constituentAll) {
       fatjet.constituentWeights = m_constituentWeights  ->at(idx);
       fatjet.constituent_pt     = m_constituent_pt      ->at(idx);
       fatjet.constituent_eta    = m_constituent_eta     ->at(idx);
       fatjet.constituent_phi    = m_constituent_phi     ->at(idx);
       fatjet.constituent_e      = m_constituent_e       ->at(idx);
     }
   
     if(m_infoSwitch.m_truth)
       {
         fatjet.truth_p4.SetPtEtaPhiE(m_truth_pt ->at(idx),
                      m_truth_eta->at(idx),
                      m_truth_phi->at(idx),
                      m_truth_m  ->at(idx));
       }
     
     if (m_infoSwitch.m_bosonCount) {
       fatjet.nTQuarks = m_nTQuarks->at(idx);
       fatjet.nHBosons = m_nHBosons->at(idx);
       fatjet.nWBosons = m_nWBosons->at(idx);
       fatjet.nZBosons = m_nZBosons->at(idx);
     }
     
     if (m_infoSwitch.m_muonCorrection) {
       fatjet.muonCorrected_pt  = m_muonCorrected_pt ->at(idx);
       fatjet.muonCorrected_eta = m_muonCorrected_eta->at(idx);
       fatjet.muonCorrected_phi = m_muonCorrected_phi->at(idx);
       fatjet.muonCorrected_m   = m_muonCorrected_m  ->at(idx);
     }
   
     for(const auto& kv : m_trkJets)
       {
         fatjet.trkJets[kv.first].clear();
         std::vector<std::vector<unsigned int>> *trkJetsIdx=m_trkJetsIdx[kv.first];
         for(unsigned int iTrkJet : trkJetsIdx->at(idx))
       {
         Jet thisTrkJet;
         kv.second->updateParticle(iTrkJet, thisTrkJet);
         fatjet.trkJets[kv.first].push_back(thisTrkJet);
       }
       }
   
     if(m_debug) std::cout << "leave FatJetContainer::updateParticle " << std::endl;
     return;
   }
   
   
   
   void FatJetContainer::setBranches(TTree *tree)
   {
     ParticleContainer::setBranches(tree);
   
     if ( m_infoSwitch.m_scales ) {
         setBranch<float>(tree, "JetConstitScaleMomentum_eta", m_JetConstitScaleMomentum_eta);
         setBranch<float>(tree, "JetConstitScaleMomentum_phi", m_JetConstitScaleMomentum_phi);
         setBranch<float>(tree, "JetConstitScaleMomentum_m", m_JetConstitScaleMomentum_m);
         setBranch<float>(tree, "JetConstitScaleMomentum_pt", m_JetConstitScaleMomentum_pt);
   
         setBranch<float>(tree, "JetEMScaleMomentum_eta", m_JetEMScaleMomentum_eta);
         setBranch<float>(tree, "JetEMScaleMomentum_phi", m_JetEMScaleMomentum_phi);
         setBranch<float>(tree, "JetEMScaleMomentum_m", m_JetEMScaleMomentum_m);
         setBranch<float>(tree, "JetEMScaleMomentum_pt", m_JetEMScaleMomentum_pt);
     }
   
     if ( m_infoSwitch.m_area ) {
       setBranch<float>(tree, "GhostArea", m_GhostArea);
       setBranch<float>(tree, "ActiveArea", m_ActiveArea);
       setBranch<float>(tree, "VoronoiArea", m_VoronoiArea);
       setBranch<float>(tree, "ActiveArea4vec_pt", m_ActiveArea4vec_pt);
       setBranch<float>(tree, "ActiveArea4vec_eta", m_ActiveArea4vec_eta);
       setBranch<float>(tree, "ActiveArea4vec_phi", m_ActiveArea4vec_phi);
       setBranch<float>(tree, "ActiveArea4vec_m", m_ActiveArea4vec_m);
   
     }
   
     if ( m_infoSwitch.m_substructure ) {
       setBranch<float>(tree, "Split12",      m_Split12);
       setBranch<float>(tree, "Split23",      m_Split23);
       setBranch<float>(tree, "Split34",      m_Split34);
       setBranch<float>(tree, "tau1_wta",     m_tau1_wta);
       setBranch<float>(tree, "tau2_wta",     m_tau2_wta);
       setBranch<float>(tree, "tau3_wta",     m_tau3_wta);
       setBranch<float>(tree, "tau21_wta",    m_tau21_wta);
       setBranch<float>(tree, "tau32_wta",    m_tau32_wta);
       setBranch<float>(tree, "ECF1",         m_ECF1);
       setBranch<float>(tree, "ECF2",         m_ECF2);
       setBranch<float>(tree, "ECF3",         m_ECF3);
       setBranch<float>(tree, "C2",           m_C2);
       setBranch<float>(tree, "D2",           m_D2);
       setBranch<float>(tree, "NTrimSubjets", m_NTrimSubjets);
       setBranch<int>  (tree, "Nclusters",    m_NClusters);
       setBranch<int>  (tree, "nTracks",      m_nTracks);
       setBranch<int>  (tree, "ungrtrk500",    m_ungrtrk500);
       setBranch<float>(tree, "EMFrac",        m_EMFrac);
       setBranch<int>  (tree, "nChargedParticles", m_nChargedParticles);
     }
   
     if ( m_infoSwitch.m_ntrimsubjets && !m_infoSwitch.m_substructure){
       setBranch<float>(tree, "NTrimSubjets", m_NTrimSubjets);
     }
   
     if ( m_infoSwitch.m_constituent) {
       setBranch<int>  (tree, "numConstituents",    m_numConstituents);
     }
   
     if ( m_infoSwitch.m_constituentAll) {
       setBranch< std::vector<float> >(tree, "constituentWeights",  m_constituentWeights);
       setBranch< std::vector<float> >(tree, "constituent_pt",      m_constituent_pt);
       setBranch< std::vector<float> >(tree, "constituent_eta",     m_constituent_eta);
       setBranch< std::vector<float> >(tree, "constituent_phi",     m_constituent_phi);
       setBranch< std::vector<float> >(tree, "constituent_e",       m_constituent_e);
     }
   
     if ( m_infoSwitch.m_truth && m_mc ) {
       setBranch<float>(tree, "truth_m"  , m_truth_m  );
       setBranch<float>(tree, "truth_pt" , m_truth_pt );
       setBranch<float>(tree, "truth_phi", m_truth_phi);
       setBranch<float>(tree, "truth_eta", m_truth_eta);
     }
       
     if ( m_infoSwitch.m_bosonCount && m_mc ) {
       setBranch< int >(tree, "nTQuarks",       m_nTQuarks);
       setBranch< int >(tree, "nHBosons",       m_nHBosons);
       setBranch< int >(tree, "nWBosons",       m_nWBosons);
       setBranch< int >(tree, "nZBosons",       m_nZBosons);
     }
     if (m_infoSwitch.m_muonCorrection) {
       setBranch<float> (tree, "muonCorrected_pt" , m_muonCorrected_pt );
       setBranch<float> (tree, "muonCorrected_eta", m_muonCorrected_eta);
       setBranch<float> (tree, "muonCorrected_phi", m_muonCorrected_phi);
       setBranch<float> (tree, "muonCorrected_m"  , m_muonCorrected_m  );
     }
   
     for(const auto& kv : m_trkJets)
       {
         kv.second->setBranches(tree);
         setBranch< std::vector<unsigned int> >(tree, "trkJetsIdx_"+kv.first, m_trkJetsIdx[kv.first]);
       }
   
     return;
   }
   
   
   void FatJetContainer::clear()
   {
   
     ParticleContainer::clear();
   
     if ( m_infoSwitch.m_scales ) {
       m_JetConstitScaleMomentum_eta   ->clear();
       m_JetConstitScaleMomentum_phi   ->clear();
       m_JetConstitScaleMomentum_m     ->clear();
       m_JetConstitScaleMomentum_pt    ->clear();
   
       m_JetEMScaleMomentum_eta        ->clear();
       m_JetEMScaleMomentum_phi        ->clear();
       m_JetEMScaleMomentum_m          ->clear();
       m_JetEMScaleMomentum_pt         ->clear();
     }
   
     if ( m_infoSwitch.m_area ) {
       m_GhostArea->clear();
       m_ActiveArea->clear();
       m_VoronoiArea->clear();
       m_ActiveArea4vec_pt->clear();
       m_ActiveArea4vec_eta->clear();
       m_ActiveArea4vec_phi->clear();
       m_ActiveArea4vec_m->clear();
     }
   
     if ( m_infoSwitch.m_substructure ) {
       m_Split12     ->clear();
       m_Split23     ->clear();
       m_Split34     ->clear();
       m_tau1_wta    ->clear();
       m_tau2_wta    ->clear();
       m_tau3_wta    ->clear();
       m_tau21_wta   ->clear();
       m_tau32_wta   ->clear();
       m_ECF1        ->clear();
       m_ECF2        ->clear();
       m_ECF3        ->clear();
       m_C2          ->clear();
       m_D2          ->clear();
       m_NTrimSubjets->clear();
       m_NClusters   ->clear();
       m_nTracks     ->clear();
       m_ungrtrk500    ->clear();
       m_EMFrac        ->clear();
       m_nChargedParticles ->clear();
     }
   
     if ( m_infoSwitch.m_ntrimsubjets && !m_infoSwitch.m_substructure ){
       m_NTrimSubjets->clear();
     }
   
     if ( m_infoSwitch.m_constituent) {
       m_numConstituents->clear();
     }
   
     if ( m_infoSwitch.m_constituentAll) {
       m_constituentWeights->clear();
       m_constituent_pt    ->clear();
       m_constituent_eta   ->clear();
       m_constituent_phi   ->clear();
       m_constituent_e     ->clear();
     }
   
     if ( m_infoSwitch.m_truth && m_mc ) {
       m_truth_m  ->clear();
       m_truth_pt ->clear();
       m_truth_phi->clear();
       m_truth_eta->clear();
     }
     
     if ( m_infoSwitch.m_bosonCount && m_mc) {
       m_nTQuarks->clear();
       m_nHBosons->clear();
       m_nWBosons->clear();
       m_nZBosons->clear();
     }
   
     if ( m_infoSwitch.m_muonCorrection) {
       m_muonCorrected_pt ->clear();
       m_muonCorrected_eta->clear();
       m_muonCorrected_phi->clear();
       m_muonCorrected_m  ->clear();
     }
     
     for(const std::pair< std::string, std::vector<std::vector<unsigned int>>* >& kv : m_trkJetsIdx)
       {
         m_trkJets   [kv.first]->clear();
         m_trkJetsIdx[kv.first]->clear();
       }
   
     return;
   }
   
   void FatJetContainer::FillFatJet( const xAOD::Jet* jet, int pvLocation ){
     return FillFatJet(static_cast<const xAOD::IParticle*>(jet), pvLocation);
   }
   
   void FatJetContainer::FillFatJet( const xAOD::IParticle* particle, int pvLocation ){
   
     ParticleContainer::FillParticle(particle);
   
     const xAOD::Jet* fatjet=dynamic_cast<const xAOD::Jet*>(particle);
   
     if( m_infoSwitch.m_scales ){
       static SG::AuxElement::ConstAccessor<float> acc_JetConstitScaleMomentum_eta("JetConstitScaleMomentum_eta");
       safeFill<float, float, xAOD::Jet>(fatjet, acc_JetConstitScaleMomentum_eta, m_JetConstitScaleMomentum_eta, -999);
       static SG::AuxElement::ConstAccessor<float> acc_JetConstitScaleMomentum_phi("JetConstitScaleMomentum_phi");
       safeFill<float, float, xAOD::Jet>(fatjet, acc_JetConstitScaleMomentum_phi, m_JetConstitScaleMomentum_phi, -999);
       static SG::AuxElement::ConstAccessor<float> acc_JetConstitScaleMomentum_m("JetConstitScaleMomentum_m");
       safeFill<float, float, xAOD::Jet>(fatjet, acc_JetConstitScaleMomentum_m, m_JetConstitScaleMomentum_m, -999, m_units);
       static SG::AuxElement::ConstAccessor<float> acc_JetConstitScaleMomentum_pt("JetConstitScaleMomentum_pt");
       safeFill<float, float, xAOD::Jet>(fatjet, acc_JetConstitScaleMomentum_pt, m_JetConstitScaleMomentum_pt, -999, m_units);
   
       static SG::AuxElement::ConstAccessor<float> acc_JetEMScaleMomentum_eta("JetEMScaleMomentum_eta");
       safeFill<float, float, xAOD::Jet>(fatjet, acc_JetEMScaleMomentum_eta, m_JetEMScaleMomentum_eta, -999);
       static SG::AuxElement::ConstAccessor<float> acc_JetEMScaleMomentum_phi("JetEMScaleMomentum_phi");
       safeFill<float, float, xAOD::Jet>(fatjet, acc_JetEMScaleMomentum_phi, m_JetEMScaleMomentum_phi, -999);
       static SG::AuxElement::ConstAccessor<float> acc_JetEMScaleMomentum_m("JetEMScaleMomentum_m");
       safeFill<float, float, xAOD::Jet>(fatjet, acc_JetEMScaleMomentum_m, m_JetEMScaleMomentum_m, -999, m_units);
       static SG::AuxElement::ConstAccessor<float> acc_JetEMScaleMomentum_pt("JetEMScaleMomentum_pt");
       safeFill<float, float, xAOD::Jet>(fatjet, acc_JetEMScaleMomentum_pt, m_JetEMScaleMomentum_pt, -999, m_units);
     }
   
     if ( m_infoSwitch.m_area ) {
       static SG::AuxElement::ConstAccessor<float> acc_GhostArea("GhostArea");
       safeFill<float, float, xAOD::Jet>(fatjet, acc_GhostArea, m_GhostArea, -999);
       static SG::AuxElement::ConstAccessor<float> acc_ActiveArea("ActiveArea");
       safeFill<float, float, xAOD::Jet>(fatjet, acc_ActiveArea, m_ActiveArea, -999);
       static SG::AuxElement::ConstAccessor<float> acc_VoronoiArea("VoronoiArea");
       safeFill<float, float, xAOD::Jet>(fatjet, acc_VoronoiArea, m_VoronoiArea, -999);
       static SG::AuxElement::ConstAccessor<float> acc_ActiveArea4vec_pt("ActiveArea4vec_pt");
       safeFill<float, float, xAOD::Jet>(fatjet, acc_ActiveArea4vec_pt, m_ActiveArea4vec_pt, -999, m_units);
       static SG::AuxElement::ConstAccessor<float> acc_ActiveArea4vec_eta("ActiveArea4vec_eta");
       safeFill<float, float, xAOD::Jet>(fatjet, acc_ActiveArea4vec_eta, m_ActiveArea4vec_eta, -999);
       static SG::AuxElement::ConstAccessor<float> acc_ActiveArea4vec_phi("ActiveArea4vec_phi");
       safeFill<float, float, xAOD::Jet>(fatjet, acc_ActiveArea4vec_phi, m_ActiveArea4vec_phi, -999);
       static SG::AuxElement::ConstAccessor<float> acc_ActiveArea4vec_m("ActiveArea4vec_m");
       safeFill<float, float, xAOD::Jet>(fatjet, acc_ActiveArea4vec_m, m_ActiveArea4vec_m, -999, m_units);
     }
   
     if( m_infoSwitch.m_substructure ){
       static SG::AuxElement::ConstAccessor<float> acc_Split12("Split12");
       safeFill<float, float, xAOD::Jet>(fatjet, acc_Split12, m_Split12, -999, m_units);
   
       static SG::AuxElement::ConstAccessor<float> acc_Split23("Split23");
       safeFill<float, float, xAOD::Jet>(fatjet, acc_Split23, m_Split23, -999, m_units);
   
       static SG::AuxElement::ConstAccessor<float> acc_Split34("Split34");
       safeFill<float, float, xAOD::Jet>(fatjet, acc_Split34, m_Split34, -999, m_units);
   
       static SG::AuxElement::ConstAccessor<float> acc_tau1_wta ("Tau1_wta");
       safeFill<float, float, xAOD::Jet>(fatjet, acc_tau1_wta, m_tau1_wta, -999);
   
       static SG::AuxElement::ConstAccessor<float> acc_tau2_wta ("Tau2_wta");
       safeFill<float, float, xAOD::Jet>(fatjet, acc_tau2_wta, m_tau2_wta, -999);
   
       static SG::AuxElement::ConstAccessor<float> acc_tau3_wta ("Tau3_wta");
       safeFill<float, float, xAOD::Jet>(fatjet, acc_tau3_wta, m_tau3_wta, -999);
   
       static SG::AuxElement::ConstAccessor<float> acc_tau21_wta ("Tau21_wta");
       if(acc_tau21_wta.isAvailable( *fatjet )){
         m_tau21_wta->push_back( acc_tau21_wta( *fatjet ) );
       } else if ( acc_tau1_wta.isAvailable( *fatjet ) && acc_tau2_wta.isAvailable( *fatjet ) ) {
         m_tau21_wta->push_back( acc_tau2_wta( *fatjet ) / acc_tau1_wta( *fatjet ) );
       } else { m_tau21_wta->push_back( -999 ); }
   
       static SG::AuxElement::ConstAccessor<float> acc_tau32_wta ("Tau32_wta");
       if(acc_tau32_wta.isAvailable( *fatjet )){
         m_tau32_wta->push_back( acc_tau32_wta( *fatjet ) );
       } else if ( acc_tau2_wta.isAvailable( *fatjet ) && acc_tau3_wta.isAvailable( *fatjet ) ) {
         m_tau32_wta->push_back( acc_tau3_wta( *fatjet ) / acc_tau2_wta( *fatjet ) );
       } else { m_tau32_wta->push_back( -999 ); }
   
       static SG::AuxElement::ConstAccessor<int> acc_NClusters ("MyNClusters");
       safeFill<int, int, xAOD::Jet>(fatjet, acc_NClusters, m_NClusters, -999);
   
       static SG::AuxElement::ConstAccessor<float> acc_ECF1 ("ECF1");
       safeFill<float, float, xAOD::Jet>(fatjet, acc_ECF1, m_ECF1, -999, m_units);
   
       static SG::AuxElement::ConstAccessor<float> acc_ECF2("ECF2");
       safeFill<float, float, xAOD::Jet>(fatjet, acc_ECF2, m_ECF2, -999, m_units);
   
       static SG::AuxElement::ConstAccessor<float> acc_ECF3 ("ECF3");
       safeFill<float, float, xAOD::Jet>(fatjet, acc_ECF3, m_ECF3, -999, m_units);
   
       static SG::AuxElement::ConstAccessor<int> NTrimSubjets("NTrimSubjets");
       safeFill<int, float, xAOD::Jet>(fatjet, NTrimSubjets, m_NTrimSubjets, -999);
   
       static SG::AuxElement::ConstAccessor<float> acc_D2 ("D2");
       if( acc_D2.isAvailable( *fatjet ) ) {
         m_D2->push_back( acc_D2( *fatjet ));
       } else if (acc_ECF1.isAvailable( *fatjet ) && acc_ECF2.isAvailable( *fatjet ) && acc_ECF3.isAvailable( *fatjet )){
         float e2=(acc_ECF2( *fatjet )/(acc_ECF1( *fatjet )*acc_ECF1( *fatjet )));
         float e3=(acc_ECF3( *fatjet )/(acc_ECF1( *fatjet )*acc_ECF1( *fatjet )*acc_ECF1( *fatjet )));
         m_D2->push_back( e3/(e2*e2*e2) );
       } else{ m_D2->push_back(-999); }
   
       static SG::AuxElement::ConstAccessor<float> acc_C2("C2");
       if(acc_C2.isAvailable(*fatjet)){
         m_C2->push_back(acc_C2(*fatjet));
       } else if( acc_ECF1.isAvailable(*fatjet) && acc_ECF2.isAvailable(*fatjet) && acc_ECF3.isAvailable(*fatjet)){
         m_C2->push_back( acc_ECF3(*fatjet)*acc_ECF1(*fatjet)/pow(acc_ECF2(*fatjet),2.0));
       } else{ m_C2->push_back(-999); }
   
       static SG::AuxElement::ConstAccessor<int> acc_GhostTrackCount("GhostTrackCount");
       if( acc_GhostTrackCount.isAvailable( *fatjet ) ) {
         m_nTracks->push_back( acc_GhostTrackCount( *fatjet ));
       } else { m_nTracks->push_back(-999); }
   
       static SG::AuxElement::ConstAccessor<ElementLink<xAOD::JetContainer>> acc_parent("Parent");
       if (acc_parent.isAvailable(*fatjet)) {
         ElementLink<xAOD::JetContainer> fatjetParentLink = acc_parent(*fatjet);
         if (fatjetParentLink.isValid()) {
           const xAOD::Jet* fatjetParent {*fatjetParentLink};
           static SG::AuxElement::ConstAccessor< std::vector<int> > acc_NumTrkPt500("NumTrkPt500");
           if ( acc_NumTrkPt500.isAvailable( *fatjetParent ) ) {
             m_ungrtrk500->push_back( acc_NumTrkPt500( *fatjetParent )[pvLocation] );
           } else { 
         //Perhaps the case if we are dealing with reclustered jets
         int sumUngrtrk500 = 0;
         const xAOD::Jet* subjet(nullptr);
         for(auto constit: fatjet->getConstituents()){
               subjet = static_cast<const xAOD::Jet*>(constit->rawConstituent());
               if (subjet->type() != xAOD::Type::Jet) continue;
           if ( acc_NumTrkPt500.isAvailable( *subjet ) ) sumUngrtrk500+=acc_NumTrkPt500( *subjet )[pvLocation];
         }
             m_ungrtrk500->push_back( sumUngrtrk500 );
           }
         } else {
           m_ungrtrk500->push_back(-999);
         }
       } else {
         m_ungrtrk500->push_back(-999);
       }
   
       m_EMFrac->push_back(GetEMFrac (*fatjet));
       
       static SG::AuxElement::ConstAccessor<int> acc_nChargedParticles("nChargedParticles");
       if( acc_nChargedParticles.isAvailable( *fatjet ) ) {
         m_nChargedParticles->push_back( acc_nChargedParticles( *fatjet ));
       } else { m_nChargedParticles->push_back(-999); }
   
     }
   
     if ( m_infoSwitch.m_ntrimsubjets && !m_infoSwitch.m_substructure ){
       static SG::AuxElement::ConstAccessor<int> NTrimSubjets("NTrimSubjets");
       safeFill<int, float, xAOD::Jet>(fatjet, NTrimSubjets, m_NTrimSubjets, -999);
     }
   
     if( m_infoSwitch.m_constituent ){
       m_numConstituents->push_back( fatjet->numConstituents() );
     }
   
   
     if( m_infoSwitch.m_constituentAll ){
       m_constituentWeights->push_back( fatjet->getAttribute< std::vector<float> >( "constituentWeights" ) );
       std::vector<float> pt;
       std::vector<float> eta;
       std::vector<float> phi;
       std::vector<float> e;
       xAOD::JetConstituentVector consVec = fatjet->getConstituents();
   
       if( consVec.isValid() ) {
   
         // use the example provided in
         // http://acode-browser.usatlas.bnl.gov/lxr/source/atlas/Event/xAOD/xAODJet/xAODJet/JetConstituentVector.h
         xAOD::JetConstituentVector::iterator constit = consVec.begin();
         xAOD::JetConstituentVector::iterator constitE = consVec.end();
         for( ; constit != constitE; constit++){
       pt. push_back( constit->pt() / m_units );
       eta.push_back( constit->eta() );
       phi.push_back( constit->phi() );
       e.  push_back( constit->e() / m_units  );
         }
       }
       m_constituent_pt ->push_back( pt  );
       m_constituent_eta->push_back( eta );
       m_constituent_phi->push_back( phi );
       m_constituent_e  ->push_back( e   );
     }
   
     if ( m_infoSwitch.m_truth && m_mc ) {
       const xAOD::Jet* truthJet = HelperFunctions::getLink<xAOD::Jet>( fatjet, "GhostTruthAssociationLink" );
       if(truthJet) {
         m_truth_pt->push_back ( truthJet->pt() / m_units );
         m_truth_eta->push_back( truthJet->eta() );
         m_truth_phi->push_back( truthJet->phi() );
         m_truth_m->push_back  ( truthJet->m() / m_units );
       } else {
         m_truth_pt->push_back ( -999 );
         m_truth_eta->push_back( -999 );
         m_truth_phi->push_back( -999 );
         m_truth_m->push_back  ( -999 );
       }
   
     }
   
     if(m_infoSwitch.m_bosonCount && m_mc){
   
       const xAOD::Jet* fatjet_parent = fatjet; // Trimmed jet area will be used for leading calo-jet if parent link fails    
   
       try
         {
       auto el = fatjet->auxdata<ElementLink<xAOD::JetContainer> >("Parent");
       if(el.isValid())
         fatjet_parent = (*el);
       else
         Warning("execute()", "Invalid link to \"Parent\" from leading calo-jet");
         }
       catch(...)
         {
       Warning("execute()", "Unable to fetch \"Parent\" link from leading calo-jet");
         }
   
       static SG::AuxElement::ConstAccessor< int > truthfatjet_TQuarks("GhostTQuarksFinalCount");
       safeFill<int, int, xAOD::Jet>(fatjet_parent, truthfatjet_TQuarks, m_nTQuarks, -999);
   
       static SG::AuxElement::ConstAccessor< int > truthfatjet_WBosons("GhostWBosonsCount");
       safeFill<int, int, xAOD::Jet>(fatjet_parent, truthfatjet_WBosons, m_nWBosons, -999);
   
       static SG::AuxElement::ConstAccessor< int > truthfatjet_ZBosons("GhostZBosonsCount");
       safeFill<int, int, xAOD::Jet>(fatjet_parent, truthfatjet_ZBosons, m_nZBosons, -999);
   
       static SG::AuxElement::ConstAccessor< int > truthfatjet_HBosons("GhostHBosonsCount");
       safeFill<int, int, xAOD::Jet>(fatjet_parent, truthfatjet_HBosons, m_nHBosons, -999);
     }
   
     if (m_infoSwitch.m_muonCorrection) {
       static const SG::AuxElement::ConstAccessor<TLorentzVector> acc_correctedFatJets_tlv("correctedFatJets_tlv");
         m_muonCorrected_pt ->push_back(acc_correctedFatJets_tlv(*fatjet).Pt() / m_units);
         m_muonCorrected_eta->push_back(acc_correctedFatJets_tlv(*fatjet).Eta());
         m_muonCorrected_phi->push_back(acc_correctedFatJets_tlv(*fatjet).Phi());
         m_muonCorrected_m  ->push_back(acc_correctedFatJets_tlv(*fatjet).M()  / m_units);
     }
   
   
     //
     // Associated track jets
     //
     if( !m_infoSwitch.m_trackJetNames.empty() ){
   
       // Find the fat jet parent
       const xAOD::Jet* fatjet_parent = fatjet; // Trimmed jet area will be used for leading calo-jet if parent link fails
   
       try{
         auto el = fatjet->auxdata<ElementLink<xAOD::JetContainer> >("Parent");
         if(el.isValid())
       fatjet_parent = (*el);
         else
       Warning("execute()", "Invalid link to \"Parent\" from leading calo-jet");
       }
       catch (...){
         Warning("execute()", "Unable to fetch \"Parent\" link from leading calo-jet");
       }
   
       // Associate the different track jet collections
       std::vector<const xAOD::Jet*> assotrkjets;
       for(const auto& trackJetName : m_infoSwitch.m_trackJetNames)
         {
       try{
         assotrkjets = fatjet_parent->getAssociatedObjects<xAOD::Jet>(trackJetName);
         std::sort( assotrkjets.begin(), assotrkjets.end(), HelperFunctions::sort_pt );
       }
       catch (...){
         Warning("execute()", "Unable to fetch \"%s\" link from leading calo-jet", trackJetName.data());
       }
   
       std::vector<unsigned int> trkJetsIdx;
       for(auto TrackJet : assotrkjets){
         if(!SelectTrackJet(TrackJet)) continue;
         trkJetsIdx.push_back(m_trkJets[trackJetName]->m_n);
         m_trkJets[trackJetName]->FillJet(TrackJet, 0 , 0);
       }
       m_trkJetsIdx[trackJetName]->push_back(trkJetsIdx);
         }
     }
   
     return;
   }
   
   
   bool FatJetContainer::SelectTrackJet(const xAOD::Jet* TrackJet)
   {
     if( TrackJet->pt() < m_trackJetPtCut )            return false;
     if( fabs(TrackJet->eta()) > m_trackJetEtaCut )    return false;
     if( TrackJet->numConstituents() < 2 )             return false;
   
     return true;
   }
   
   float FatJetContainer::GetEMFrac(const xAOD::Jet& jet) {
       float eInSample = 0.; 
       float eInSampleFull = 0.; 
       float emfrac = 0.; 
       const xAOD::JetConstituentVector constituents = jet.getConstituents();
       if (!constituents.isValid()){
         //ATH_MSG_WARNING("Unable to retrieve valid constituents from parent of large R jet");
         return -1.;
       }
       for (const auto& constituent : constituents) {
         if(!constituent) continue;
         if(constituent->rawConstituent()->type()!=xAOD::Type::CaloCluster) {
           //ATH_MSG_WARNING("Tried to call fillEperSamplingCluster with a jet constituent that is not a cluster!");
           continue;
         }
         const xAOD::CaloCluster* constit = static_cast<const xAOD::CaloCluster*>(constituent->rawConstituent());  
         if(!constit) continue;
         for (int s=0;s<CaloSampling::Unknown; s++){
             eInSampleFull += constit->eSample(CaloSampling::CaloSample(s));
         }
         eInSample += constit->eSample(CaloSampling::EMB1);
         eInSample += constit->eSample(CaloSampling::EMB2);
         eInSample += constit->eSample(CaloSampling::EMB3);
         eInSample += constit->eSample(CaloSampling::EME1);
         eInSample += constit->eSample(CaloSampling::EME2);
         eInSample += constit->eSample(CaloSampling::EME3);
         eInSample += constit->eSample(CaloSampling::FCAL1);
       }
       emfrac  = eInSample/eInSampleFull;
       if ( emfrac > 1.0 ) emfrac = 1.;
       else if ( emfrac < 0.0 ) emfrac = 0.;
   
       return emfrac;
   }