Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:14:13 

0001 #include "RhoFluct.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 RhoFluct::RhoFluct
0041      ( const std::string& _outputfilename
0042      , const float        _jet_R   
0043      )
0044   : SubsysReco("RhoFluct")
0045   , m_outputfilename { _outputfilename }
0046   , jet_R            { _jet_R          }
0047 { }
0048 
0049 RhoFluct::~RhoFluct()
0050 { 
0051   for (unsigned int i = 0; i < m_input_Sub1.size(); ++i) delete m_input_Sub1[i];
0052   m_input_Sub1.clear();
0053 
0054   for (unsigned int i = 0; i < m_input_rhoA.size(); ++i) delete m_input_rhoA[i];
0055   m_input_rhoA.clear();
0056 }
0057 
0058 int RhoFluct::Init(PHCompositeNode* topNode)
0059 {
0060   if (Verbosity() >= RhoFluct::VERBOSITY_SOME)
0061     std::cout << "RhoFluct::Init - Outoput to " << m_outputfilename << std::endl;
0062 
0063   std::cout << " File out: " << m_outputfilename << std::endl;
0064   PHTFileServer::get().open(m_outputfilename, "RECREATE");
0065 
0066   //Tree
0067   m_T = new TTree("T", "RhoFluct Tree");
0068 
0069   //      int m_event;
0070   /* m_T->Branch("id",                &m_id); */
0071   m_T->Branch("rho",         &m_rho);
0072   m_T->Branch("rho_sigma",   &m_rho_sigma);
0073   m_T->Branch("cent",        &m_cent);
0074   m_T->Branch("cent_mdb",    &m_cent_mdb);
0075   m_T->Branch("cent_epd",    &m_cent_epd);
0076   m_T->Branch("impactparam", &m_impactparam);
0077 
0078 
0079   // embedded objects
0080   m_T->Branch("emb_1TeV_phi", &m_1TeV_phi);
0081   m_T->Branch("emb_1TeV_eta", &m_1TeV_eta);
0082 
0083   //Sub1 Jets
0084   m_T->Branch("sub1_ismatched", &m_Sub1_ismatched);
0085   m_T->Branch("sub1JetPhi",     &m_Sub1JetPhi);
0086   m_T->Branch("sub1JetEta",     &m_Sub1JetEta);
0087   m_T->Branch("sub1JetPt",      &m_Sub1JetPt);
0088   m_T->Branch("sub1Jet_delPt",  &m_Sub1Jet_delpt);
0089 
0090   //rhoA jets
0091   m_T->Branch("rhoA_ismatched", &m_rhoA_ismatched);
0092   m_T->Branch("rhoAJetPhi",     &m_rhoAJetPhi);
0093   m_T->Branch("rhoAJetEta",     &m_rhoAJetEta);
0094   m_T->Branch("rhoAJetPt",      &m_rhoAJetPt);
0095   m_T->Branch("rhoAJetArea",    &m_rhoAJetArea);
0096   m_T->Branch("rhoAJetPtLessRhoA",    &m_rhoAJetPtLessRhoA);
0097   m_T->Branch("rhoAJet_delPt",  &m_rhoAJet_delpt);
0098 
0099   return Fun4AllReturnCodes::EVENT_OK;
0100 }
0101 
0102 int RhoFluct::End(PHCompositeNode* topNode)
0103 {
0104   std::cout << "RhoFluct::End - Output to " << m_outputfilename << std::endl;
0105   PHTFileServer::get().cd(m_outputfilename);
0106 
0107   m_T->Write();
0108   return Fun4AllReturnCodes::EVENT_OK;
0109 }
0110 
0111 int RhoFluct::InitRun(PHCompositeNode* topNode)
0112 {
0113   topNode->print();
0114   std::cout << " Input Selections:" << std::endl;
0115   for (unsigned int i = 0; i < m_input_rhoA.size(); ++i) m_input_rhoA[i]->identify();
0116   for (unsigned int i = 0; i < m_input_Sub1.size(); ++i) m_input_Sub1[i]->identify();
0117   return Fun4AllReturnCodes::EVENT_OK;
0118 }
0119 
0120 int RhoFluct::process_event(PHCompositeNode* topNode)
0121 {
0122   ++nevent;
0123 
0124   // cut for only event with good vertices
0125   GlobalVertexMap *mapVtx = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0126   if (mapVtx->size() == 0) return Fun4AllReturnCodes::EVENT_OK;
0127 
0128   // centrality
0129   CentralityInfo* cent_node = findNode::getClass<CentralityInfo>(topNode, "CentralityInfo");
0130   if (cent_node)
0131   {
0132     m_cent_epd  =  cent_node->get_centile(CentralityInfo::PROP::epd_NS);
0133     m_cent_mdb  =  cent_node->get_centile(CentralityInfo::PROP::mbd_NS);
0134     m_impactparam =  cent_node->get_quantity(CentralityInfo::PROP::bimp);
0135     m_cent =  cent_node->get_centile(CentralityInfo::PROP::bimp);
0136   }
0137 
0138   const double ghost_max_rap = 1.1;
0139   const double ghost_R = 0.01;
0140   const double max_jet_rap = 1.1; // - jet_R;
0141 
0142   //get the SUB1 jets
0143   m_1TeV_phi =  rand3.Uniform(-M_PI, M_PI);
0144   m_1TeV_eta =  rand3.Uniform(-0.7,  0.7);
0145 
0146   const int TOW_PRINT_INT = 8000;
0147 
0148   // build the towers for rhoA
0149   int n_inputs {0};
0150   std::vector<fastjet::PseudoJet> calo_pseudojets;
0151   for (auto inp : m_input_rhoA) {
0152     if (Verbosity()>TOW_PRINT_INT) std::cout << " Input Tower Set("<<n_inputs++<<"): ";
0153     float sum_pt {0.};
0154     float n_towers {0};
0155     std::vector<Jet *> particles = inp->get_input(topNode);
0156     for (unsigned int i = 0; i < particles.size(); ++i)
0157     {
0158       auto& part = particles[i];
0159       float this_e = part->get_e();
0160       if (this_e == 0.) continue;
0161       float this_px = particles[i]->get_px();
0162       float this_py = particles[i]->get_py();
0163       float this_pz = particles[i]->get_pz();
0164       if (this_e < 0)
0165       {
0166         // make energy = +1 MeV for purposes of clustering
0167         float e_ratio = 0.001 / this_e;
0168         this_e  = this_e * e_ratio;
0169         this_px = this_px * e_ratio;
0170         this_py = this_py * e_ratio;
0171         this_pz = this_pz * e_ratio;
0172       }
0173       fastjet::PseudoJet pseudojet(this_px, this_py, this_pz, this_e);
0174       /* if (pseudojet.pt() < 0.2) continue; */
0175       calo_pseudojets    .push_back(pseudojet);
0176       if (Verbosity() > TOW_PRINT_INT) {
0177         sum_pt += pseudojet.perp();
0178         ++n_towers;
0179       }
0180     }
0181     if (Verbosity()>TOW_PRINT_INT) std::cout << " sumpt(" << sum_pt <<") from ntowers("<<n_towers <<")" << std::endl;
0182     for (auto &p : particles) delete p;
0183   }
0184   // add the TeV particle
0185   auto embjet = fastjet::PseudoJet();
0186   embjet.reset_PtYPhiM(30., m_1TeV_eta, m_1TeV_phi);
0187   calo_pseudojets.push_back(embjet);
0188 
0189   fastjet::Selector jetrap         = fastjet::SelectorAbsEtaMax(max_jet_rap);
0190   fastjet::Selector not_pure_ghost = !fastjet::SelectorIsPureGhost();
0191 
0192   fastjet::Selector selection      = jetrap && not_pure_ghost;
0193   fastjet::AreaDefinition area_def( 
0194       fastjet::active_area_explicit_ghosts, fastjet::GhostedAreaSpec(ghost_max_rap, 1, ghost_R));
0195 
0196   fastjet::JetDefinition jet_def_antikt(fastjet::antikt_algorithm, jet_R);
0197 
0198 
0199   fastjet::JetDefinition jet_def_bkgd(fastjet::kt_algorithm, jet_R); // <--
0200   fastjet::Selector selector_rm2 = jetrap * (!fastjet::SelectorNHardest(2)); // <--
0201                                                                                                               //
0202   fastjet::JetMedianBackgroundEstimator bge_rm2 {selector_rm2, jet_def_bkgd, area_def};
0203   bge_rm2.set_particles(calo_pseudojets);
0204 
0205   m_rho = bge_rm2.rho();
0206   m_rho_sigma = bge_rm2.sigma();
0207 
0208   mBGE_mean_area    = bge_rm2 .mean_area();
0209   mBGE_empty_area   = bge_rm2 .empty_area();
0210   mBGE_n_empty_jets = bge_rm2 .n_empty_jets();
0211   mBGE_n_jets_used  = bge_rm2 .n_jets_used();
0212 
0213   fastjet::ClusterSequenceArea clustSeq(calo_pseudojets, jet_def_antikt, area_def);
0214   std::vector<fastjet::PseudoJet> jets 
0215     = sorted_by_pt( selection( clustSeq.inclusive_jets(5.) ));
0216 
0217   int fixmek = 0;
0218   if (Verbosity()>100 && m_cent_mdb<10) {
0219     fixmek = 0;
0220     std::cout << "  --- rhoA jets [set("<<outevent<<")]" << std::endl;
0221     for (auto jet : jets) {
0222       if (jet.perp()-m_rho*jet.area() > 5.) {
0223         auto phi = jet.phi();
0224         while (phi > M_PI) phi -= 2*M_PI;
0225          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;
0226       }
0227     }
0228   }
0229 
0230   // check that the leading jet is within R=0.4 of the embedded jet
0231   if (jets.size() == 0) {
0232     std::cout << " no jets reconstructed; something is very wrong" << std::endl;
0233     m_rhoA_ismatched    = false;
0234     m_rhoAJetPhi        = -100.;
0235     m_rhoAJetEta        = -100.;
0236     m_rhoAJetPt         = -100.;
0237     m_rhoAJetArea       = -100.;
0238     m_rhoAJetPtLessRhoA = -100.;
0239     m_rhoAJet_delpt     = -100.;
0240   } else {
0241     auto& lead_jet = jets[0];
0242     m_rhoAJetPhi = lead_jet.phi();
0243     m_rhoAJetEta = lead_jet.eta();
0244     m_rhoAJetPt  = lead_jet.perp();
0245     m_rhoAJetArea = lead_jet.area();
0246     m_rhoAJetPtLessRhoA = m_rhoAJetPt - m_rhoAJetArea*m_rho;
0247 
0248     float  dphi = TMath::Abs(lead_jet.phi_std() - m_1TeV_phi);
0249     while (dphi>M_PI) dphi = TMath::Abs(dphi-2*M_PI);
0250     float dR = TMath::Sq(dphi)+TMath::Sq(lead_jet.eta()-m_1TeV_eta);
0251     if (dR > 0.16) {
0252       m_rhoA_ismatched = false;
0253       m_rhoAJet_delpt = -100.;
0254     } else {
0255       m_rhoA_ismatched = true;
0256       m_rhoAJet_delpt = m_rhoAJetPtLessRhoA - 30.;
0257     }
0258   }
0259 
0260   // --------------------------------------
0261   // now make the jets with Sub1
0262   // --------------------------------------
0263   std::vector<fastjet::PseudoJet> sub1_pseudojets;
0264   for (auto inp : m_input_Sub1) {
0265     if (Verbosity()>TOW_PRINT_INT) std::cout << " Input Tower Set("<<n_inputs++<<"): ";
0266     float sum_pt {0.};
0267     float n_towers {0};
0268     std::vector<Jet *> particles = inp->get_input(topNode);
0269     for (unsigned int i = 0; i < particles.size(); ++i)
0270     {
0271       auto& part = particles[i];
0272       float this_e = part->get_e();
0273       if (this_e == 0.) continue;
0274       float this_px = particles[i]->get_px();
0275       float this_py = particles[i]->get_py();
0276       float this_pz = particles[i]->get_pz();
0277       if (this_e < 0)
0278       {
0279         // make energy = +1 MeV for purposes of clustering
0280         float e_ratio = 0.001 / this_e;
0281         this_e  = this_e * e_ratio;
0282         this_px = this_px * e_ratio;
0283         this_py = this_py * e_ratio;
0284         this_pz = this_pz * e_ratio;
0285       }
0286       fastjet::PseudoJet pseudojet(this_px, this_py, this_pz, this_e);
0287       /* if (pseudojet.pt() < 0.2) continue; */
0288       sub1_pseudojets    .push_back(pseudojet);
0289       if (Verbosity() > TOW_PRINT_INT) {
0290         sum_pt += pseudojet.perp();
0291         ++n_towers;
0292       }
0293     }
0294     if (Verbosity()>TOW_PRINT_INT) std::cout << " sumpt(" << sum_pt <<") from ntowers("<<n_towers <<")" << std::endl;
0295     for (auto &p : particles) delete p;
0296   }
0297   // add the TeV particle
0298   sub1_pseudojets.push_back(embjet);
0299 
0300   fastjet::ClusterSequenceArea clustSeq_sub1(sub1_pseudojets, jet_def_antikt, area_def);
0301   std::vector<fastjet::PseudoJet> jets_sub1 
0302     = sorted_by_pt( selection( clustSeq_sub1.inclusive_jets(5.) ));
0303 
0304   // check that the leading jet is within R=0.4 of the embedded jet
0305   if (jets_sub1.size() == 0) {
0306     std::cout << " no jets_sub1 reconstructed; something is very wrong" << std::endl;
0307     m_Sub1_ismatched    = false;
0308     m_Sub1JetPhi        = -100.;
0309     m_Sub1JetEta        = -100.;
0310     m_Sub1JetPt         = -100.;
0311     m_Sub1Jet_delpt     = -100.;
0312   } else {
0313     auto& lead_jet = jets_sub1[0];
0314     m_Sub1JetPhi = lead_jet.phi();
0315     m_Sub1JetEta = lead_jet.eta();
0316     m_Sub1JetPt  = lead_jet.perp();
0317 
0318     float  dphi = TMath::Abs(lead_jet.phi_std() - m_1TeV_phi);
0319     while (dphi>M_PI) dphi = TMath::Abs(dphi-2*M_PI);
0320     float dR = TMath::Sq(dphi)+TMath::Sq(lead_jet.eta()-m_1TeV_eta);
0321     if (dR > 0.16) {
0322       m_Sub1_ismatched = false;
0323       m_Sub1Jet_delpt = -100.;
0324     } else {
0325       m_Sub1_ismatched = true;
0326       m_Sub1Jet_delpt = m_Sub1JetPt - 30.;
0327     }
0328   }
0329 
0330   m_T->Fill();
0331   return Fun4AllReturnCodes::EVENT_OK;
0332 }
0333
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!