FastJetMedianBkg/src_JetMedianTree/JetRhoMedian.cc

0001 #include "JetRhoMedian.h"
0002 #include "fastjet/AreaDefinition.hh"
0003 #include "fastjet/ClusterSequenceArea.hh"
0004 #include "fastjet/Selector.hh"
0005 #include "fastjet/tools/BackgroundEstimatorBase.hh"
0006 #include "fastjet/tools/JetMedianBackgroundEstimator.hh"
0007
0008
0009 #include <TFile.h>
0010 #include <TH1F.h>
0011 #include <TH2F.h>
0012 #include <TRandom3.h>
0013 #include <TString.h>
0014 #include <TTree.h>
0015 #include <TH2D.h>
0016 #include <TVector3.h>
0017 #include <algorithm>
0018 #include <cassert>
0019 #include <centrality/CentralityInfo.h>
0020 #include <cmath>
0021 #include <fastjet/ClusterSequence.hh>
0022 #include <fastjet/JetDefinition.hh>
0023 #include <fastjet/PseudoJet.hh>
0024 #include <fun4all/Fun4AllReturnCodes.h>
0025 #include <fun4all/Fun4AllServer.h>
0026 #include <fun4all/PHTFileServer.h>
0027 #include <g4eval/JetEvalStack.h>
0028 #include <jetbase/JetInput.h>
0029 #include <jetbase/JetMapv1.h>
0030 #include <iostream>
0031 #include <limits>
0032 #include <phool/PHCompositeNode.h>
0033 #include <phool/getClass.h>
0034 #include <stdexcept>
0035 #include <trackbase_historic/SvtxTrackMap.h>
0036 #include <globalvertex/GlobalVertexMap.h>
0037
0038 // using namespace fastjet;
0039
0040 JetRhoMedian::JetRhoMedian
0041      ( const std::string& _outputfilename
0042      , const float        _jet_R
0043      , const std::string& truthjetname
0044      , const std::string& sub1jetname
0045      , const float min_lead_truth_pt
0046      )
0047   : SubsysReco("JetRhoMedian")
0048   , m_outputfilename { _outputfilename }
0049   , jet_R            { _jet_R          }
0050   , m_truthJetName   { truthjetname    }
0051   , m_sub1JetName    { sub1jetname     }
0052   , m_min_lead_truth_pt { min_lead_truth_pt }
0053   , m_T            { nullptr         }
0054   , m_inputs         {}
0055 { }
0056
0057 JetRhoMedian::~JetRhoMedian()
0058 {
0059   for (unsigned int i = 0; i < m_inputs.size(); ++i) delete m_inputs[i];
0060   m_inputs.clear();
0061 }
0062
0063 int JetRhoMedian::Init(PHCompositeNode* topNode)
0064 {
0065   if (Verbosity() >= JetRhoMedian::VERBOSITY_SOME)
0066     std::cout << "JetRhoMedian::Init - Outoput to " << m_outputfilename << std::endl;
0067
0068   std::cout << " File out: " << m_outputfilename << std::endl;
0069   PHTFileServer::get().open(m_outputfilename, "RECREATE");
0070
0071   //Tree
0072   m_T = new TTree("T", "JetRhoMedian Tree");
0073
0074   //      int m_event;
0075   /* m_T->Branch("id",                &m_id); */
0076   m_T->Branch("rho",               &m_rho);
0077   m_T->Branch("rho_sigma",         &m_rho_sigma);
0078   m_T->Branch("cent",          &m_cent);
0079   m_T->Branch("cent_mdb",          &m_cent_mdb);
0080   m_T->Branch("cent_epd",          &m_cent_epd);
0081   m_T->Branch("impactparam",       &m_impactparam);
0082
0083
0084   //Truth Jets
0085   m_T->Branch("TruthJetEta",  &m_TruthJetEta);
0086   m_T->Branch("TruthJetPhi",  &m_TruthJetPhi);
0087   m_T->Branch("TruthJetPt",   &m_TruthJetPt);
0088
0089   //Sub1 Jets
0090   m_T->Branch("sub1JetEta",  &m_Sub1JetEta);
0091   m_T->Branch("sub1JetPhi",  &m_Sub1JetPhi);
0092   m_T->Branch("sub1JetPt",   &m_Sub1JetPt);
0093
0094   m_T->Branch("rhoAJetEta",        &m_CaloJetEta);
0095   m_T->Branch("rhoAJetPhi",        &m_CaloJetPhi);
0096   m_T->Branch("rhoAJetPt",          &m_CaloJetPt);
0097   m_T->Branch("rhoAJetPtLessRhoA", &m_CaloJetPtLessRhoA);
0098   m_T->Branch("rhoAJetArea",       &m_CaloJetArea);
0099
0100   return Fun4AllReturnCodes::EVENT_OK;
0101 }
0102
0103 int JetRhoMedian::End(PHCompositeNode* topNode)
0104 {
0105   std::cout << "JetRhoMedian::End - Output to " << m_outputfilename << std::endl;
0106   PHTFileServer::get().cd(m_outputfilename);
0107
0108   m_T->Write();
0109   return Fun4AllReturnCodes::EVENT_OK;
0110 }
0111
0112 int JetRhoMedian::InitRun(PHCompositeNode* topNode)
0113 {
0114   topNode->print();
0115   std::cout << " Input Selections:" << std::endl;
0116   for (unsigned int i = 0; i < m_inputs.size(); ++i) m_inputs[i]->identify();
0117   return Fun4AllReturnCodes::EVENT_OK;
0118 }
0119
0120 int JetRhoMedian::process_event(PHCompositeNode* topNode)
0121 {
0122   clear_vectors();
0123   ++nevent;
0124
0125   // cut for only event with good vertices
0126   GlobalVertexMap *mapVtx = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0127   if (mapVtx->size() == 0) return Fun4AllReturnCodes::EVENT_OK;
0128
0129   // centrality
0130   CentralityInfo* cent_node = findNode::getClass<CentralityInfo>(topNode, "CentralityInfo");
0131   if (cent_node)
0132   {
0133     m_cent_epd  =  cent_node->get_centile(CentralityInfo::PROP::epd_NS);
0134     m_cent_mdb  =  cent_node->get_centile(CentralityInfo::PROP::mbd_NS);
0135     m_impactparam =  cent_node->get_quantity(CentralityInfo::PROP::bimp);
0136     m_cent =  cent_node->get_centile(CentralityInfo::PROP::bimp);
0137   }
0138
0139   const double ghost_max_rap = 1.1;
0140   const double ghost_R = 0.01;
0141   const double max_jet_rap = 1.1; // - jet_R;
0142
0143
0144   /* fastjet::JetDefinition jet_def(fastjet::kt_algorithm, 0.4);     //  JET DEFINITION */
0145   /* const double ghost_max_rap { 1.0 }; */
0146
0147   /* fastjet::AreaDefinition area_def_bkgd( */
0148   /*     fastjet::active_area_explicit_ghosts, fastjet::GhostedAreaSpec(ghost_max_rap, 1, ghost_R)); */
0149   /* fastjet::Selector selector_rm2 = jetrap * (!fastjet::SelectorNHardest(2)); // <-- */
0150
0151   //interface to truth jets
0152   JetMap* jetsMC = findNode::getClass<JetMap>(topNode, m_truthJetName);
0153   if (!jetsMC )
0154   {
0155     std::cout
0156       << "MyJetAnalysis::process_event - Error can not find DST Truth JetMap node "
0157       << m_truthJetName << std::endl;
0158     exit(-1);
0159   }
0160
0161   bool isfirst = true;
0162   int fixmek = 0;
0163   auto vjets = jetsMC->vec();
0164   std::sort(vjets.begin(), vjets.end(), [](Jet* a, Jet* b){return a->get_pt() > b->get_pt();});
0165   if (Verbosity() > 100) {
0166     if (m_cent_mdb<10) {
0167       std::cout << " NEW EVENT[" << (outevent++) << "]" << std::endl;
0168       std::cout << " ---  Truth Jets [set(" << outevent << ")]" <<  std::endl;
0169     }
0170   }
0171   for (auto& jet : vjets) {
0172     float pt  = jet->get_pt();
0173     if  (pt < 5.) continue;
0174     float eta = jet->get_eta();
0175     if (eta < -max_jet_rap || eta > max_jet_rap) continue;
0176     if (Verbosity()>100 && m_cent_mdb<10) {
0177       std::cout << Form("k[%i] pt:phi:eta[%5.2f:%5.2f:%5.2f]", fixmek++, jet->get_pt(), jet->get_phi(), jet->get_eta()) << std::endl;
0178     }
0179
0180     if (isfirst) {
0181       isfirst = false;
0182       if (pt < m_min_lead_truth_pt) return Fun4AllReturnCodes::EVENT_OK;
0183     }
0184     m_TruthJetPt .push_back(jet->get_pt());
0185     m_TruthJetEta.push_back(jet->get_eta());
0186     m_TruthJetPhi.push_back(jet->get_phi());
0187   }
0188
0189   //get the SUB1 jets
0190     if (Verbosity() > 100 && m_cent_mdb<10) {
0191       std::cout << "  --- Sub1 Jets [set(" << outevent << ")]" <<  std::endl;
0192        fixmek = 0;
0193     }
0194   JetMap* jets_sub1 = findNode::getClass<JetMap>(topNode, m_sub1JetName);
0195   if (jets_sub1 != nullptr) {
0196     vjets = jets_sub1->vec();
0197     std::sort(vjets.begin(), vjets.end(), [](Jet* a, Jet* b){return a->get_pt() > b->get_pt();});
0198     fixmek = 0;
0199     if (Verbosity()>100 && m_cent_mdb<10) std::cout << "SUB1 jets" << std::endl;
0200     for (auto& jet : jets_sub1->vec()) {
0201       // will return jets in order of descending pT
0202       float pt  = jet->get_pt();
0203       float eta = jet->get_eta();
0204       if (eta < -max_jet_rap || eta > max_jet_rap) continue;
0205       if  (pt < 5.) continue;
0206       if (Verbosity()>100 && m_cent_mdb<10) std::cout << Form("k[%i] pt:phi:eta[%5.2f:%5.2f:%5.2f]", fixmek++, jet->get_pt(), jet->get_phi(), jet->get_eta()) << std::endl;
0207       m_Sub1JetPt .push_back(jet->get_pt());
0208       m_Sub1JetEta.push_back(jet->get_eta());
0209       m_Sub1JetPhi.push_back(jet->get_phi());
0210     }
0211   }
0212
0213   const int TOW_PRINT_INT = 8000;
0214
0215   int n_inputs {0};
0216   std::vector<fastjet::PseudoJet> calo_pseudojets;
0217   for (auto inp : m_inputs) {
0218     if (Verbosity()>TOW_PRINT_INT) std::cout << " Input Tower Set("<<n_inputs++<<"): ";
0219     float sum_pt {0.};
0220     float n_towers {0};
0221     std::vector<Jet *> particles = inp->get_input(topNode);
0222     for (unsigned int i = 0; i < particles.size(); ++i)
0223     {
0224       auto& part = particles[i];
0225       float this_e = part->get_e();
0226       if (this_e == 0.) continue;
0227       float this_px = particles[i]->get_px();
0228       float this_py = particles[i]->get_py();
0229       float this_pz = particles[i]->get_pz();
0230       if (this_e < 0)
0231       {
0232         // make energy = +1 MeV for purposes of clustering
0233         float e_ratio = 0.001 / this_e;
0234         this_e  = this_e * e_ratio;
0235         this_px = this_px * e_ratio;
0236         this_py = this_py * e_ratio;
0237         this_pz = this_pz * e_ratio;
0238       }
0239       fastjet::PseudoJet pseudojet(this_px, this_py, this_pz, this_e);
0240       /* if (pseudojet.pt() < 0.2) continue; */
0241       calo_pseudojets    .push_back(pseudojet);
0242       if (Verbosity() > TOW_PRINT_INT) {
0243         sum_pt += pseudojet.perp();
0244         ++n_towers;
0245       }
0246     }
0247     if (Verbosity()>TOW_PRINT_INT) std::cout << " sumpt(" << sum_pt <<") from ntowers("<<n_towers <<")" << std::endl;
0248     for (auto &p : particles) delete p;
0249   }
0250
0251
0252   fastjet::Selector jetrap         = fastjet::SelectorAbsEtaMax(max_jet_rap);
0253   fastjet::Selector not_pure_ghost = !fastjet::SelectorIsPureGhost();
0254
0255   fastjet::Selector selection      = jetrap && not_pure_ghost;
0256   fastjet::AreaDefinition area_def(
0257       fastjet::active_area_explicit_ghosts, fastjet::GhostedAreaSpec(ghost_max_rap, 1, ghost_R));
0258
0259   fastjet::JetDefinition jet_def_antikt(fastjet::antikt_algorithm, jet_R);
0260
0261
0262   fastjet::JetDefinition jet_def_bkgd(fastjet::kt_algorithm, jet_R); // <--
0263   fastjet::Selector selector_rm2 = jetrap * (!fastjet::SelectorNHardest(2)); // <--
0264                                                                                                               //
0265   fastjet::JetMedianBackgroundEstimator bge_rm2 {selector_rm2, jet_def_bkgd, area_def};
0266   bge_rm2.set_particles(calo_pseudojets);
0267
0268   m_rho = bge_rm2.rho();
0269   m_rho_sigma = bge_rm2.sigma();
0270
0271   mBGE_mean_area    = bge_rm2 .mean_area();
0272   mBGE_empty_area   = bge_rm2 .empty_area();
0273   mBGE_n_empty_jets = bge_rm2 .n_empty_jets();
0274   mBGE_n_jets_used  = bge_rm2 .n_jets_used();
0275
0276
0277   fastjet::ClusterSequenceArea clustSeq(calo_pseudojets, jet_def_antikt, area_def);
0278   std::vector<fastjet::PseudoJet> jets
0279     = sorted_by_pt( selection( clustSeq.inclusive_jets(5.) ));
0280
0281   if (Verbosity()>100 && m_cent_mdb<10) {
0282     fixmek = 0;
0283     std::cout << "  --- rhoA jets [set("<<outevent<<")]" << std::endl;
0284     for (auto jet : jets) {
0285       if (jet.perp()-m_rho*jet.area() > 5.) {
0286         auto phi = jet.phi();
0287         while (phi > M_PI) phi -= 2*M_PI;
0288          std::cout << Form("k[%i] pt:phi:eta[%5.2f:%5.2f:%5.2f]", fixmek++, (jet.perp()-m_rho*jet.area()), phi, jet.eta()) << std::endl;
0289       }
0290     }
0291   }
0292
0293   for (auto jet : jets) {
0294     if (false && Verbosity() >= JetRhoMedian::VERBOSITY_SOME) {
0295         auto cst = fastjet::sorted_by_pt(jet.constituents());
0296         int i =0;
0297         for (auto p : cst)  {
0298           if (p.is_pure_ghost()) continue;
0299           std::cout << Form(" %3i|%4.2f",i++,p.perp());
0300           if (i%7==0) std::cout << std::endl;
0301         }
0302         std::cout << std::endl;
0303     }
0304
0305     m_CaloJetEta        .push_back( jet.eta()     );
0306     m_CaloJetPhi        .push_back( jet.phi_std() );
0307     m_CaloJetPt         .push_back( jet.pt()      );
0308     m_CaloJetPtLessRhoA .push_back( jet.pt()      - m_rho * jet.area());
0309     m_CaloJetArea       .push_back( jet.area());
0310   }
0311
0312
0313   m_T->Fill();
0314   return Fun4AllReturnCodes::EVENT_OK;
0315 }
0316
0317 void JetRhoMedian::clear_vectors() {
0318   /* if (nevent % 1 == 0) cout << " z0 " << m_CaloJetArea.size(); */
0319   m_CaloJetEta        .clear();
0320   m_CaloJetPhi        .clear();
0321   m_CaloJetPt          .clear();
0322   m_CaloJetPtLessRhoA .clear();
0323   m_CaloJetArea       .clear();
0324
0325   /* if (nevent % 1 == 0) cout << " " << m_TruthJetPt.size(); */
0326   m_TruthJetEta       .clear();
0327   m_TruthJetPhi       .clear();
0328   m_TruthJetPt        .clear();
0329
0330   m_Sub1JetEta       .clear();
0331   m_Sub1JetPhi       .clear();
0332   m_Sub1JetPt        .clear();
0333
0334 }