RandomConeAna/src/RandomConeTree.cc

0001 #include "RandomConeTree.h"
0002
0003 // random cones
0004 #include <randomcones/RandomCone.h>
0005 #include <randomcones/RandomConev1.h>
0006
0007 // tower rho
0008 #include <towerrho/TowerRho.h>
0009 #include <towerrho/TowerRhov1.h>
0010
0011 // fun4all includes
0012 #include <fun4all/Fun4AllReturnCodes.h>
0013 #include <fun4all/PHTFileServer.h>
0014
0015 // phool includes
0016 #include <phool/PHCompositeNode.h>
0017 #include <phool/getClass.h>
0018
0019 // jet reco
0020 #include <jetbase/Jet.h>
0021 #include <jetbase/JetMap.h>
0022 #include <jetbase/Jetv1.h>
0023 #include <jetbase/JetMapv1.h>
0024 #include <jetbase/JetContainer.h>
0025 #include <jetbase/JetContainerv1.h>
0026 #include <jetbase/Jetv2.h>
0027
0028
0029 // centrality
0030 #include <centrality/CentralityInfo.h>
0031
0032 // global vertex
0033 #include <globalvertex/GlobalVertex.h>
0034 #include <globalvertex/GlobalVertexMap.h>
0035 #include <globalvertex/GlobalVertexMapv1.h>
0036
0037 // mbd
0038 #include <mbd/MbdOut.h>
0039
0040 // root includes
0041 #include <TTree.h>
0042
0043 // standard includes
0044 #include <cmath>
0045 #include <utility>
0046 #include <cstdlib>
0047 #include <iostream>
0048 #include <vector>
0049 #include <string>
0050 #include <array>
0051
0052
0053 RandomConeTree::RandomConeTree(const std::string& outputfilename)
0054  : SubsysReco("RandomConeTree")
0055   , m_outputfile_name(outputfilename)
0056 {
0057 }
0058
0059 RandomConeTree::~RandomConeTree()
0060 {
0061   // probably don't need to delete these but
0062   if(m_run_info_tree) delete m_run_info_tree;
0063   if(m_event_tree) delete m_event_tree;
0064 }
0065
0066 int RandomConeTree::Init(PHCompositeNode *topNode)
0067 {
0068   // create output file
0069   PHTFileServer::get().open(m_outputfile_name, "RECREATE");
0070
0071   // create trees
0072   m_run_info_tree = new TTree("run_info", "RunInfo");
0073   m_event_tree = new TTree("event_info", "EventInfo");
0074
0075   // run info tree
0076   if(m_do_gvtx_cut)
0077   {
0078     m_run_info_tree->Branch("zvtx_low", &m_gvtx_cut.first, "zvtx_low/F");
0079     m_run_info_tree->Branch("zvtx_high", &m_gvtx_cut.second, "zvtx_high/F");
0080   }
0081   if(m_MC_do_event_selection)
0082   {
0083     m_run_info_tree->Branch("MC_event_selection_jetpT_low", &m_MC_event_selection_jetpT_range.first, "MC_event_selection_jetpT_low/F");
0084     m_run_info_tree->Branch("MC_event_selection_jetpT_high", &m_MC_event_selection_jetpT_range.second, "MC_event_selection_jetpT_high/F");
0085   }
0086   if(m_MC_add_weight_to_ttree)
0087   {
0088     m_run_info_tree->Branch("MC_event_weight", &m_MC_event_weight, "MC_event_weight/D");
0089   }
0090
0091   // fill run info tree
0092   m_run_info_tree->Fill();
0093
0094   // event info tree
0095   m_event_tree->Branch("event_id", &m_event_id, "event_id/I");
0096   if(m_do_cent_info)
0097   {
0098     m_event_tree->Branch("centrality", &m_centrality, "centrality/I");
0099     if(!m_do_data) m_event_tree->Branch("impactparam", &m_impactparam, "impactparam/F");
0100     m_event_tree->Branch("mbd_NS", &m_mbd_NS, "mbd_NS/D");
0101   }
0102   m_event_tree->Branch("zvtx", &m_zvtx, "zvtx/F");
0103
0104   // add rho nodes
0105   for(unsigned int i = 0; i < m_rho_nodes.size(); ++i)
0106   {
0107     std::string name = m_rho_nodes.at(i)+"_rho";
0108     m_event_tree->Branch(name.c_str(), &m_rhos[i], (name + "/F").c_str());
0109     name = m_rho_nodes.at(i)+"_sigma";
0110     m_event_tree->Branch(name.c_str(), &m_sigma_rhos[i], (name + "/F").c_str());
0111   }
0112
0113   // add random cone nodes
0114   for(unsigned int i = 0; i < m_random_cone_nodes.size(); ++i)
0115   {
0116     std::string name = m_random_cone_nodes.at(i)+"_R";
0117     m_event_tree->Branch(name.c_str(), &m_cone_R[i], (name + "/F").c_str());
0118     name = m_random_cone_nodes.at(i)+"_eta";
0119     m_event_tree->Branch(name.c_str(), &m_cone_etas[i], (name + "/F").c_str());
0120     name = m_random_cone_nodes.at(i)+"_phi";
0121     m_event_tree->Branch(name.c_str(), &m_cone_phis[i], (name + "/F").c_str());
0122     name = m_random_cone_nodes.at(i)+"_pt";
0123     m_event_tree->Branch(name.c_str(), &m_cone_pts[i], (name + "/F").c_str());
0124     name = m_random_cone_nodes.at(i)+"_nclustered";
0125     m_event_tree->Branch(name.c_str(), &m_cone_nclustereds[i], (name + "/I").c_str());
0126   }
0127
0128
0129   m_event_id = -1;
0130
0131
0132   return Fun4AllReturnCodes::EVENT_OK;
0133
0134 }
0135
0136 int RandomConeTree::ResetEvent(PHCompositeNode *topNode)
0137 {
0138   // reset event variables
0139   m_centrality = -1;
0140   m_impactparam = NAN;
0141   m_zvtx = NAN;
0142   m_mbd_NS = NAN;
0143
0144   for(unsigned int i = 0; i < m_rho_nodes.size(); ++i)
0145   {
0146     m_rhos[i] = NAN;
0147     m_sigma_rhos[i] = NAN;
0148   }
0149
0150   for(unsigned int i = 0; i < m_random_cone_nodes.size(); ++i)
0151   {
0152     m_cone_R[i] = NAN;
0153     m_cone_etas[i] = NAN;
0154     m_cone_phis[i] = NAN;
0155     m_cone_pts[i] = NAN;
0156     m_cone_nclustereds[i] = -1;
0157   }
0158
0159   return Fun4AllReturnCodes::EVENT_OK;
0160 }
0161
0162 int RandomConeTree::process_event(PHCompositeNode *topNode)
0163 {
0164
0165   ++m_event_id;
0166
0167   // event selection
0168   if(!GoodEvent(topNode)) return Fun4AllReturnCodes::ABORTEVENT;
0169
0170   // Centrality info
0171   if(m_do_cent_info)  GetCentInfo(topNode);
0172
0173
0174   // rho info
0175   for(unsigned int i = 0; i < m_rho_nodes.size(); ++i)
0176   {
0177     TowerRho* towerrho = findNode::getClass<TowerRhov1>(topNode, m_rho_nodes.at(i));
0178     if (!towerrho)
0179     {
0180       std::cout << "RandomConeTree - Error can not find towerrho " << m_rho_nodes.at(i) << std::endl;
0181       exit(-1);
0182     }
0183     m_rhos[i] = towerrho->get_rho();
0184     m_sigma_rhos[i] = towerrho->get_sigma();
0185   }
0186
0187   // random cone info
0188   for(unsigned int i = 0; i < m_random_cone_nodes.size(); ++i)
0189   {
0190     RandomCone* rcone = findNode::getClass<RandomConev1>(topNode, m_random_cone_nodes.at(i));
0191     if (!rcone)
0192     {
0193       std::cout << "RandomConeTree - Error can not find random cone " << m_random_cone_nodes.at(i) << std::endl;
0194       exit(-1);
0195     }
0196     m_cone_R[i] = rcone->get_R();
0197     m_cone_etas[i] = rcone->get_eta();
0198     m_cone_phis[i] = rcone->get_phi();
0199     m_cone_pts[i] = rcone->get_pt();
0200     m_cone_nclustereds[i] = rcone->n_clustered();
0201   }
0202
0203   // fill event tree
0204   m_event_tree->Fill();
0205
0206   return Fun4AllReturnCodes::EVENT_OK;
0207
0208
0209
0210   // //==================================
0211   // //get towerrho
0212
0213   // TowerRho* towerrho_area = findNode::getClass<TowerRhov1>(topNode, "TowerRho_AREA");
0214   // if (!towerrho_area)
0215   // {
0216   //   std::cout << "RandomConeTree - Error can not find towerrho area" << std::endl;
0217   //   exit(-1);
0218   // }
0219
0220   // TowerRho* towerrho_mult = findNode::getClass<TowerRhov1>(topNode, "TowerRho_MULT");
0221   // if (!towerrho_mult)
0222   // {
0223   //   std::cout << "RandomConeTree - Error can not find towerrho mult" << std::endl;
0224   //   exit(-1);
0225   // }
0226
0227   // m_rho_area = towerrho_area->get_rho();
0228   // m_rho_area_sigma = towerrho_area->get_sigma();
0229   // m_rho_mult = towerrho_mult->get_rho();
0230   // m_rho_mult_sigma = towerrho_mult->get_sigma();
0231
0232
0233   // RandomCone* rcone_basic = findNode::getClass<RandomConev1>(topNode, "RandomCone_TOWER_basic_r04");
0234   // RandomCone* rcone_randomize_etaphi = findNode::getClass<RandomConev1>(topNode, "RandomCone_TOWER_randomize_etaphi_r04");
0235   // RandomCone* rcone_avoid_lead_jet = findNode::getClass<RandomConev1>(topNode, "RandomCone_TOWER_avoid_lead_jet_r04");
0236   // RandomCone* rcone_sub_basic = findNode::getClass<RandomConev1>(topNode, "RandomCone_TOWER_SUB1_basic_r04");
0237   // RandomCone* rcone_sub_randomize_etaphi = findNode::getClass<RandomConev1>(topNode, "RandomCone_TOWER_SUB1_randomize_etaphi_r04");
0238   // RandomCone* rcone_sub_avoid_lead_jet = findNode::getClass<RandomConev1>(topNode, "RandomCone_TOWER_SUB1_avoid_lead_jet_r04");
0239
0240   // if(!rcone_basic){
0241   //   std::cout << "RandomConeTree::process_event - Error can not find basic random cone" << std::endl;
0242   //   exit(-1);
0243   // }
0244   // if(!rcone_randomize_etaphi){
0245   //   std::cout << "RandomConeTree::process_event - Error can not find randomize etaphi random cone" << std::endl;
0246   //   exit(-1);
0247   // }
0248   // if(!rcone_avoid_lead_jet){
0249   //   std::cout << "RandomConeTree::process_event - Error can not find avoid leading jet random cone" << std::endl;
0250   //   exit(-1);
0251   // }
0252   // if(!rcone_sub_basic){
0253   //   std::cout << "RandomConeTree::process_event - Error can not find basic random cone" << std::endl;
0254   //   exit(-1);
0255   // }
0256   // if(!rcone_sub_randomize_etaphi){
0257   //   std::cout << "RandomConeTree::process_event - Error can not find randomize etaphi random cone" << std::endl;
0258   //   exit(-1);
0259   // }
0260   // if(!rcone_sub_avoid_lead_jet){
0261   //   std::cout << "RandomConeTree::process_event - Error can not find avoid leading jet random cone" << std::endl;
0262   //   exit(-1);
0263   // }
0264
0265   // m_cone_area = rcone_basic->get_area();
0266
0267   // std::vector<Jet *> inputs;  // owns memory
0268   // for (auto &_input : _inputs)
0269   // {
0270   //   std::vector<Jet *> parts = _input->get_input(topNode);
0271   //   for (auto &part : parts)
0272   //   {
0273   //     inputs.push_back(part);
0274   //     inputs.back()->set_id(inputs.size() - 1);  // unique ids ensured
0275   //   }
0276   // }
0277
0278   // std::vector<Jet *> iter_inputs;  // owns memory
0279   // for (auto &_iter_input : _iter_inputs)
0280   // {
0281   //   std::vector<Jet *> parts = _iter_input->get_input(topNode);
0282   //   for (auto &part : parts)
0283   //   {
0284   //     iter_inputs.push_back(part);
0285   //     iter_inputs.back()->set_id(iter_inputs.size() - 1);  // unique ids ensured
0286   //   }
0287   // }
0288
0289
0290
0291   // // get pseudojets
0292   // std::vector<fastjet::PseudoJet> pseudojets = jets_to_pseudojets(inputs, false);
0293   // std::vector<fastjet::PseudoJet> iter_pseudojets = jets_to_pseudojets(iter_inputs, false);
0294
0295   // // select random cone
0296   // m_cone_eta = rcone_basic->get_eta();
0297   // m_cone_phi = rcone_basic->get_phi();
0298   // m_cone_from_RC_pt = rcone_basic->get_pt();
0299   // m_cone_nclustered_from_RC = rcone_basic->n_clustered();
0300
0301   // // control cone
0302   // m_cone_nclustered = 0;
0303   // m_cone_pt = 0;
0304   // for (unsigned int ipart = 0; ipart < pseudojets.size(); ++ipart)
0305   // {
0306   //   // calculate the distance between the particle and the cone center
0307   //   float deta = m_cone_eta - pseudojets[ipart].eta();
0308   //   float dphi = m_cone_phi - pseudojets[ipart].phi_std();
0309   //   if (dphi > M_PI) dphi -= 2 * M_PI;
0310   //   if (dphi < -M_PI) dphi += 2 * M_PI;
0311   //   float dr = sqrt(deta * deta + dphi * dphi);
0312   //   if (dr < rcone_basic->get_R())
0313   //   {
0314   //     m_cone_pt += pseudojets[ipart].pt();
0315   //     m_cone_nclustered++;
0316   //   }
0317   // }
0318
0319   // // randomize cone
0320   // m_randomized_cone_eta = rcone_randomize_etaphi->get_eta();
0321   // m_randomized_cone_phi = rcone_randomize_etaphi->get_phi();
0322   // m_randomized_cone_from_RC_pt = rcone_randomize_etaphi->get_pt();
0323   // m_randomized_cone_nclustered_from_RC = rcone_randomize_etaphi->n_clustered();
0324
0325   // // avoid leading jet cone
0326   // m_avoid_leading_cone_eta = rcone_avoid_lead_jet->get_eta();
0327   // m_avoid_leading_cone_phi = rcone_avoid_lead_jet->get_phi();
0328   // m_avoid_leading_cone_from_RC_pt = rcone_avoid_lead_jet->get_pt();
0329   // m_avoid_leading_cone_nclustered_from_RC = rcone_avoid_lead_jet->n_clustered();
0330
0331   // // iter sub cone
0332   // m_cone_eta_iter = rcone_sub_basic->get_eta();
0333   // m_cone_phi_iter = rcone_sub_basic->get_phi();
0334   // m_cone_from_RC_pt_iter = rcone_sub_basic->get_pt();
0335   // m_cone_nclustered_from_RC_iter = rcone_sub_basic->n_clustered();
0336
0337
0338   // // iter sub cone
0339   // m_cone_pt_iter = 0;
0340   // for (unsigned int ipart = 0; ipart < iter_pseudojets.size(); ++ipart)
0341   // {
0342   //   // calculate the distance between the particle and the cone center
0343   //   float deta = m_cone_eta_iter - iter_pseudojets[ipart].eta();
0344   //   float dphi = m_cone_phi_iter - iter_pseudojets[ipart].phi_std();
0345   //   if (dphi > M_PI) dphi -= 2 * M_PI;
0346   //   if (dphi < -M_PI) dphi += 2 * M_PI;
0347   //   float dr = sqrt(deta * deta + dphi * dphi);
0348   //   if (dr < rcone_sub_basic->get_R())
0349   //   {
0350   //     m_cone_pt_iter += iter_pseudojets[ipart].pt();
0351   //     m_cone_nclustered_iter++;
0352   //   }
0353
0354   // }
0355
0356   // // randomize cone
0357   // m_randomized_cone_eta_iter = rcone_sub_randomize_etaphi->get_eta();
0358   // m_randomized_cone_phi_iter = rcone_sub_randomize_etaphi->get_phi();
0359   // m_randomized_cone_from_RC_pt_iter = rcone_sub_randomize_etaphi->get_pt();
0360   // m_randomized_cone_nclustered_from_RC_iter = rcone_sub_randomize_etaphi->n_clustered();
0361
0362   // // avoid leading jet cone
0363   // m_avoid_leading_cone_eta_iter = rcone_sub_avoid_lead_jet->get_eta();
0364   // m_avoid_leading_cone_phi_iter = rcone_sub_avoid_lead_jet->get_phi();
0365   // m_avoid_leading_cone_from_RC_pt_iter = rcone_sub_avoid_lead_jet->get_pt();
0366   // m_avoid_leading_cone_nclustered_from_RC_iter = rcone_sub_avoid_lead_jet->n_clustered();
0367
0368
0369   // //==================================
0370   // // Fill tree
0371   // //==================================
0372   // m_tree->Fill();
0373
0374   // return Fun4AllReturnCodes::EVENT_OK;
0375 }
0376
0377 int RandomConeTree::End(PHCompositeNode *topNode)
0378 {
0379   std::cout << "RandomConeTree::End - Output to " << m_outputfile_name << std::endl;
0380     // write tree to file
0381   PHTFileServer::get().cd(m_outputfile_name);
0382   // m_tree->Write();
0383   // write file to disk
0384   PHTFileServer::get().write(m_outputfile_name);
0385   std::cout << "RandomConeTree::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0386   return Fun4AllReturnCodes::EVENT_OK;
0387 }
0388
0389
0390 //===========================================================
0391 // Event selection
0392 bool RandomConeTree::GoodEvent(PHCompositeNode *topNode)
0393 {
0394   // MC event selection based on leading R = 0.4 truth jet pT
0395   if(m_MC_do_event_selection)
0396   {
0397     // get truth jet nodes
0398     JetMap *truthjets = findNode::getClass<JetMap>(topNode, "AntiKt_Truth_r04");
0399     if(!truthjets)
0400     {
0401       std::cout << "RandomConeTree::GoodEvent(PHCompositeNode *topNode) Could not find truth jet nodes" << std::endl;
0402       exit(-1); // this is a fatal error
0403     }
0404
0405     // get leading truth jet pT
0406     float leading_truth_pt = -1;
0407     for(JetMap::Iter iter = truthjets->begin(); iter != truthjets->end(); ++iter)
0408     {
0409
0410       Jet *jet = iter->second;
0411
0412       if(jet->get_pt() > leading_truth_pt) leading_truth_pt = jet->get_pt();
0413
0414     }
0415
0416     // check if event passes selection
0417     if( (leading_truth_pt < m_MC_event_selection_jetpT_range.first) ||
0418         (leading_truth_pt > m_MC_event_selection_jetpT_range.second)
0419     )
0420     {
0421
0422       if(Verbosity() > 0) std::cout << "RandomConeTree::GoodEvent(PHCompositeNode *topNode) Event failed MC event selection" << std::endl;
0423       return false;
0424
0425     }
0426   }
0427
0428   // zvtx cut
0429   if(m_do_gvtx_cut)
0430   {
0431     // get global vertex map node
0432     GlobalVertexMap *vtxMap = findNode::getClass<GlobalVertexMapv1>(topNode,"GlobalVertexMap");
0433     if (!vtxMap)
0434     {
0435       std::cout << "RandomConeTree::GoodEvent(PHCompositeNode *topNode) Could not find global vertex map node" << std::endl;
0436       exit(-1); // this is a fatal error
0437     }
0438
0439     if (!vtxMap->get(0))
0440     {
0441       if(Verbosity() > 0) std::cout << "RandomConeTree::GoodEvent(PHCompositeNode *topNode) No primary vertex" << std::endl;
0442       return false;
0443     }
0444
0445     // get zvtx
0446     m_zvtx = vtxMap->get(0)->get_z();
0447
0448     if( (m_zvtx < m_gvtx_cut.first) || (m_zvtx > m_gvtx_cut.second))
0449     {
0450       if(Verbosity() > 0) std::cout << "RandomConeTree::GoodEvent(PHCompositeNode *topNode) Vertex cut failed" << std::endl;
0451       return false;
0452     }
0453
0454   }
0455
0456   return true;
0457
0458 }
0459
0460 //===========================================================
0461 // Cent info
0462 void RandomConeTree::GetCentInfo(PHCompositeNode *topNode)
0463 {
0464   //==================================
0465   // Get centrality info
0466   //==================================
0467   CentralityInfo* cent_node = findNode::getClass<CentralityInfo>(topNode, "CentralityInfo");
0468   if (!cent_node)
0469   {
0470       std::cout << "RandomConeTree::process_event() ERROR: Can't find CentralityInfo" << std::endl;
0471       exit(-1); // this is a fatal error
0472   }
0473
0474   if (!m_do_data)
0475   {
0476       m_centrality = cent_node->get_centile(CentralityInfo::PROP::bimp);
0477       m_impactparam = cent_node->get_quantity(CentralityInfo::PROP::bimp);
0478       m_mbd_NS = cent_node->get_quantity(CentralityInfo::PROP::mbd_NS);
0479   }
0480   else
0481   {
0482       m_centrality = (int)(100*cent_node->get_centile(CentralityInfo::PROP::mbd_NS));
0483
0484       PHNodeIterator iter(topNode);
0485       PHCompositeNode *mbdNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "MBD"));
0486       if(!mbdNode)
0487         {
0488             std::cerr << Name() << "::" <<  __PRETTY_FUNCTION__ << "MBD Node missing, doing nothing." << std::endl;
0489             throw std::runtime_error("Failed to find MBD node in RandomConeTree::GetCentInfo");
0490         }
0491
0492       MbdOut *_data_MBD = findNode::getClass<MbdOut>(mbdNode, "MbdOut");
0493
0494       if(!_data_MBD)
0495         {
0496           std::cerr << Name() << "::" <<  __PRETTY_FUNCTION__ << "MbdOut Node missing, doing nothing." << std::endl;
0497           throw std::runtime_error("Failed to find MbdOut node in RandomConeTree::GetCentInfo");
0498         }
0499
0500       m_mbd_NS = _data_MBD->get_q(0) + _data_MBD->get_q(1);
0501     }
0502
0503     return ;
0504 }