Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:13:23 

0001 
0002 //module for producing a TTree with jet information for doing jet validation studies
0003 // for questions/bugs please contact Virginia Bailey vbailey13@gsu.edu
0004 #include <fun4all/Fun4AllBase.h>
0005 #include <JetValidationTC.h>
0006 #include <fun4all/Fun4AllReturnCodes.h>
0007 #include <fun4all/PHTFileServer.h>
0008 #include <fun4all/Fun4AllHistoManager.h>
0009 #include <phool/PHCompositeNode.h>
0010 #include <phool/getClass.h>
0011 #include <jetbase/JetMap.h>
0012 #include <jetbase/JetContainer.h>
0013 #include <jetbase/Jetv2.h>
0014 #include <jetbase/Jetv1.h>
0015 #include <centrality/CentralityInfo.h>
0016 #include <globalvertex/GlobalVertex.h>
0017 #include <globalvertex/GlobalVertexMap.h>
0018 #include <calobase/RawTower.h>
0019 #include <calobase/RawTowerContainer.h>
0020 #include <calobase/RawTowerGeom.h>
0021 #include <calobase/RawTowerGeomContainer.h>
0022 #include <calobase/TowerInfoContainer.h>
0023 #include <calobase/TowerInfo.h>
0024 #include <calobase/RawClusterContainer.h>
0025 #include <calobase/RawCluster.h>
0026 #include <calobase/RawClusterUtility.h>
0027 
0028 #include <ffarawobjects/Gl1Packet.h>
0029 
0030 #include <jetbackground/TowerBackground.h>
0031 
0032 #include <TTree.h>
0033 #include <iostream>
0034 #include <sstream>
0035 #include <iomanip>
0036 #include <cmath>
0037 #include <vector>
0038 #include <TLorentzVector.h>
0039 
0040 #include "fastjet/ClusterSequence.hh"
0041 #include "fastjet/contrib/SoftDrop.hh"                                                                                          
0042 
0043 using namespace fastjet;
0044 
0045 // ROOT, for histogramming.                                                                                                                                                                                 
0046 
0047 #include "TH1.h"
0048 #include "TH2.h"
0049 #include "TFile.h"
0050 #include <TTree.h>
0051 
0052 //____________________________________________________________________________..
0053 JetValidationTC::JetValidationTC(const std::string& recojetname, const std::string& truthjetname, const std::string& outputfilename):
0054   SubsysReco("JetValidation_" + recojetname + "_" + truthjetname)
0055   , m_recoJetName(recojetname)
0056   , m_truthJetName(truthjetname)
0057   , m_outputFileName(outputfilename)
0058   , m_etaRange(-1, 1)
0059   , m_ptRange(5, 100)
0060   , m_doTruthJets(0)
0061   , m_doSeeds(0)
0062   , m_doUnsubJet(0)
0063   , m_removeJetClusters(0)
0064   , m_T(nullptr)
0065   , m_event(-1)
0066   , m_nTruthJet(-1)
0067   , m_nJet(-1)
0068   , m_id()
0069   , m_nComponent()
0070   , m_eta()
0071   , m_phi()
0072   , m_e()
0073   , m_pt()
0074   , m_zg()
0075   , m_rg()
0076   , m_fe()
0077   , m_fh()
0078   , m_sub_et()
0079   , m_truthID()
0080   , m_truthNComponent()
0081   , m_truthEta()
0082   , m_truthPhi()
0083   , m_truthE()
0084   , m_truthPt()
0085   , m_eta_rawseed()
0086   , m_phi_rawseed()
0087   , m_pt_rawseed()
0088   , m_e_rawseed()
0089   , m_rawseed_cut()
0090   , m_eta_subseed()
0091   , m_phi_subseed()
0092   , m_pt_subseed()
0093   , m_e_subseed()
0094   , m_subseed_cut()
0095   
0096 {
0097   std::cout << "JetValidationTC::JetValidationTC(const std::string &name) Calling ctor" << std::endl;
0098 }
0099 
0100 //____________________________________________________________________________..
0101 JetValidationTC::~JetValidationTC()
0102 {
0103   std::cout << "JetValidationTC::~JetValidationTC() Calling dtor" << std::endl;
0104 }
0105 
0106 //____________________________________________________________________________..
0107 //____________________________________________________________________________..
0108 int JetValidationTC::Init(PHCompositeNode *topNode)
0109 {
0110   std::cout << "JetValidationTC::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0111   //PHTFileServer::get().open(m_outputFileName, "RECREATE");
0112   m_histoFileName =  m_outputFileName + "_histos.root"; 
0113   m_outputFileName = m_outputFileName + ".root";
0114   hm = new Fun4AllHistoManager(Name());
0115   outfile = new TFile(m_outputFileName.c_str(),"RECREATE");
0116   std::cout << "JetValidationTC::Init - Output to " << m_outputFileName << std::endl;
0117   m_hasJetAboveThreshold = 0;
0118 
0119   
0120   hm->setOutfileName("JetValidation_histos.root");
0121 
0122   h_pi0M = new TH1F("h_pi0M","h_pi0M",4000,0,20);
0123   h_npions = new TH1F("h_npions","h_npions",100,0,100);
0124   h_eta_phi_clusters = new TH2F("h_eta_phi_clusters","h_eta_phi_clusters",22,-1.1,1.1,63,-3.14,3.14);
0125   h_jetPionMult = new TH1F("h_jetPionMult","h_jetPionMult",100,0,100);
0126   h_photonEnergy = new TH1F("h_photonEnergy","h_photonEnergy",400,0,20);
0127 
0128 
0129   
0130 
0131   // configure Tree
0132   m_T = new TTree("T", "MyJetAnalysis Tree");
0133   m_T->Branch("m_event", &m_event, "event/I");
0134   m_T->Branch("nJet", &m_nJet, "nJet/I");
0135   m_T->Branch("cent", &m_centrality);
0136   m_T->Branch("zvtx", &m_zvtx);
0137   m_T->Branch("b", &m_impactparam);
0138   m_T->Branch("id", &m_id);
0139   m_T->Branch("nComponent", &m_nComponent);
0140 
0141   m_T->Branch("eta", &m_eta);
0142   m_T->Branch("phi", &m_phi);
0143   m_T->Branch("e", &m_e);
0144   m_T->Branch("pt", &m_pt);
0145   m_T->Branch("zg", &m_zg);
0146   m_T->Branch("rg", &m_rg);
0147   m_T->Branch("triggers",&m_triggers);
0148 
0149   if (m_doClusters) {
0150     m_T->Branch("fe",&m_fe);
0151     m_T->Branch("fh",&m_fh);
0152     m_T->Branch("f_Et_emcal", &m_fEt_emcal);
0153     m_T->Branch("f_Et_hcal", &m_fEt_hcal);
0154     m_T->Branch("constE",&m_constE);
0155     m_T->Branch("hcalE",&m_hcalE);
0156     m_T->Branch("ecalE",&m_ecalE);
0157     m_T->Branch("constEt",&m_constEt);
0158     m_T->Branch("hcalClusterEt",&m_cluster_hcalEt);
0159     m_T->Branch("ecalClusterEt",&m_cluster_ecalEt);
0160     m_T->Branch("Et_emcal",&m_Et_emcal);
0161     m_T->Branch("Et_hcal",&m_Et_hcal);
0162     m_T->Branch("ClusterTowerE",&m_ClusterTowE);
0163     m_T->Branch("EMCalClusterMult",&m_emcal_cluster_mult);
0164     m_T->Branch("HCalClusterMult",&m_hcal_cluster_mult);
0165     m_T->Branch("ClusterMult",&m_cluster_mult);
0166     m_T->Branch("HasJetAboveThreshold", &m_hasJetAboveThreshold);
0167     m_T->Branch("pion_Z",&m_pionZ);
0168     m_T->Branch("jetPionMult", &m_jetPionMult);
0169     m_T->Branch("jetPionPt", &m_jetPionPt);
0170     m_T->Branch("jetPionEta",&m_jetPionEta);
0171     m_T->Branch("jetPionPhi",&m_jetPionPhi);
0172     m_T->Branch("jetPionMass",&m_jetPionMass);
0173     m_T->Branch("jetPionEnergy",&m_jetPionEnergy);
0174     //m_T->Branch("pionMassJetsRemoved", &m_pionMassJetsRemoved);
0175     m_T->Branch("pi0Mass",&pi0cand_mass);
0176     m_T->Branch("pi0pt",&pi0cand_pt);
0177     m_T->Branch("pi0eta",&pi0cand_eta);
0178     m_T->Branch("pi0phi",&pi0cand_phi);
0179     m_T->Branch("pi0energy",&pi0cand_energy);
0180   }
0181   if(m_doUnsubJet)
0182     {
0183       m_T->Branch("pt_unsub", &m_unsub_pt);
0184       m_T->Branch("subtracted_et", &m_sub_et);
0185     }
0186   if(m_doTruthJets){
0187     m_T->Branch("nTruthJet", &m_nTruthJet);
0188     m_T->Branch("truthID", &m_truthID);
0189     m_T->Branch("truthNComponent", &m_truthNComponent);
0190     m_T->Branch("truthEta", &m_truthEta);
0191     m_T->Branch("truthPhi", &m_truthPhi);
0192     m_T->Branch("truthE", &m_truthE);
0193     m_T->Branch("truthPt", &m_truthPt);
0194   }
0195 
0196   if(m_doSeeds){
0197     m_T->Branch("rawseedEta", &m_eta_rawseed);
0198     m_T->Branch("rawseedPhi", &m_phi_rawseed);
0199     m_T->Branch("rawseedPt", &m_pt_rawseed);
0200     m_T->Branch("rawseedE", &m_e_rawseed);
0201     m_T->Branch("rawseedCut", &m_rawseed_cut);
0202     m_T->Branch("subseedEta", &m_eta_subseed);
0203     m_T->Branch("subseedPhi", &m_phi_subseed);
0204     m_T->Branch("subseedPt", &m_pt_subseed);
0205     m_T->Branch("subseedE", &m_e_subseed);
0206     m_T->Branch("subseedCut", &m_subseed_cut);
0207   }
0208   
0209   return Fun4AllReturnCodes::EVENT_OK;
0210 }
0211 
0212 //____________________________________________________________________________..
0213 int JetValidationTC::InitRun(PHCompositeNode *topNode)
0214 {
0215   std::cout << "JetValidationTC::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0216   
0217   return Fun4AllReturnCodes::EVENT_OK;
0218 }
0219 
0220 //____________________________________________________________________________..
0221 int JetValidationTC::process_event(PHCompositeNode *topNode)
0222 {
0223   //  std::cout << "JetValidationTC::process_event(PHCompositeNode *topNode) Processing Event" << std::endl;
0224   ++m_event;
0225 
0226   // interface to reco jets
0227   JetContainer* jets = findNode::getClass<JetContainer>(topNode, m_recoJetName);
0228   if (!jets)
0229     {
0230       std::cout
0231     << "MyJetAnalysis::process_event - Error can not find DST Reco JetContainer node "
0232     << m_recoJetName << std::endl;
0233       exit(-1);
0234     }
0235 
0236   //interface to truth jets
0237   //JetMap* jetsMC = findNode::getClass<JetMap>(topNode, m_truthJetName);
0238   JetContainer* jetsMC = findNode::getClass<JetContainer>(topNode, m_truthJetName);
0239   if (!jetsMC && m_doTruthJets)
0240     {
0241       std::cout
0242     << "MyJetAnalysis::process_event - Error can not find DST Truth JetMap node "
0243     << m_truthJetName << std::endl;
0244       exit(-1);
0245     }
0246   
0247   // interface to jet seeds
0248   JetContainer* seedjetsraw = findNode::getClass<JetContainer>(topNode, "AntiKt_TowerInfo_HIRecoSeedsRaw_r02");
0249   if (!seedjetsraw && m_doSeeds)
0250     {
0251       std::cout
0252     << "MyJetAnalysis::process_event - Error can not find DST raw seed jets "
0253 << std::endl;
0254       exit(-1);
0255     }
0256 
0257   JetContainer* seedjetssub = findNode::getClass<JetContainer>(topNode, "AntiKt_TowerInfo_HIRecoSeedsSub_r02");
0258   if (!seedjetssub && m_doSeeds)
0259     {
0260       std::cout
0261 << "MyJetAnalysis::process_event - Error can not find DST subtracted seed jets "
0262 << std::endl;
0263       exit(-1);
0264     }
0265 
0266   //centrality
0267   CentralityInfo* cent_node = findNode::getClass<CentralityInfo>(topNode, "CentralityInfo");
0268   if (!cent_node)
0269     {
0270       std::cout
0271         << "MyJetAnalysis::process_event - Error can not find centrality node "
0272         << std::endl;
0273       exit(-1);
0274     }
0275   
0276   //zvertex
0277   GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0278   if (!vertexmap)
0279     {
0280       std::cout
0281         << "MyJetAnalysis::process_event - Error can not find global vertex  node "
0282         << std::endl;
0283       exit(-1);
0284     }
0285   if (vertexmap->empty())
0286     {
0287       std::cout
0288         << "MyJetAnalysis::process_event - global vertex node is empty "
0289         << std::endl;
0290       return Fun4AllReturnCodes::ABORTEVENT;
0291     }
0292 
0293     GlobalVertex *vtx = vertexmap->begin()->second;
0294   m_zvtx = vtx->get_z();
0295   CLHEP::Hep3Vector vertex;
0296   vertex.set(vtx->get_x(),vtx->get_y(),vtx->get_z());
0297 
0298   
0299  
0300   //calorimeter towers
0301   TowerInfoContainer *towersEM3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC_RETOWER_SUB1");
0302   TowerInfoContainer *towersIH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN_SUB1");
0303   TowerInfoContainer *towersOH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT_SUB1");
0304   RawTowerGeomContainer *tower_geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0305   RawTowerGeomContainer *tower_geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0306   RawClusterContainer *clustersHCAL = findNode::getClass<RawClusterContainer>(topNode, "TOPOCLUSTER_HCAL");
0307   RawClusterContainer *clustersEMCAL = findNode::getClass<RawClusterContainer>(topNode, "TOPOCLUSTER_EMCAL");
0308 
0309   if(!towersEM3 || !towersIH3 || !towersOH3){
0310     std::cout
0311       <<"MyJetAnalysis::process_event - Error can not find raw tower node "
0312       << std::endl;
0313     exit(-1);
0314   }
0315 
0316   if(!tower_geom || !tower_geomOH){
0317     std::cout
0318       <<"MyJetAnalysis::process_event - Error can not find raw tower geometry "
0319       << std::endl;
0320     exit(-1);
0321   }
0322 
0323   //underlying event
0324   TowerBackground *background = findNode::getClass<TowerBackground>(topNode, "TowerInfoBackground_Sub2");
0325   if(!background){
0326     std::cout<<"Can't get background. Exiting"<<std::endl;
0327     return Fun4AllReturnCodes::EVENT_OK;
0328   }
0329 
0330   //get the event centrality/impact parameter from HIJING
0331   m_centrality =  cent_node->get_centile(CentralityInfo::PROP::bimp);
0332   m_impactparam =  cent_node->get_quantity(CentralityInfo::PROP::bimp);
0333 
0334   //get reco jets
0335   m_nJet = 0;
0336   m_nPionJets = 0;
0337   float background_v2 = 0;
0338   float background_Psi2 = 0;
0339   if(m_doUnsubJet)
0340     {
0341       background_v2 = background->get_v2();
0342       background_Psi2 = background->get_Psi2();
0343     }
0344 
0345   reconstruct_pi0s(topNode);
0346 
0347   for (auto jet : *jets)
0348     {
0349 
0350       if(jet->get_pt() < 1) continue; // to remove noise jets
0351       if(jet->get_pt() < m_ptRange.first) continue; // remove jets outside of pt range
0352 
0353       m_id.push_back(jet->get_id());
0354       m_nComponent.push_back(jet->size_comp());
0355       m_eta.push_back(jet->get_eta());
0356       m_phi.push_back(jet->get_phi());
0357       m_e.push_back(jet->get_e());
0358       m_pt.push_back(jet->get_pt());
0359       m_zg.push_back(jet->get_property(jets->property_index(Jet::PROPERTY::prop_zg)));
0360       m_rg.push_back(jet->get_property(jets->property_index(Jet::PROPERTY::prop_Rg)));
0361       m_Et.push_back(jet->get_et());
0362 
0363       //get the pions within dR of the jet axis
0364       std::cout << "JetValidationTC::process_event: looking for pions" << std::endl;
0365       findJetPions(jet);
0366 
0367       //get the jet constituents and calculate signal ratios
0368 
0369       float _fe = 0.;
0370       float _fh = 0.;
0371       float _fEt_h = 0.;
0372       float _fEt_e = 0.;
0373 
0374       if (m_doClusters) {
0375 
0376         
0377 
0378         if (!clustersEMCAL || !clustersHCAL) {
0379           std::cout << "JetValidationTC::process_event - Error can not find RawClusterContainer"
0380           << std::endl;
0381           exit(-1);
0382         }
0383 
0384         std::vector<float> this_jet_hcal_et;
0385         std::vector<float> this_jet_emcal_et;
0386         std::vector<float> this_jet_const_et;
0387 
0388         int this_hcal_cluster_mult = 0;
0389         int this_emcal_cluster_mult = 0;
0390         int this_cluster_mult = 0;
0391 
0392         for (auto comp: jet->get_comp_vec()) {
0393 
0394           this_cluster_mult++;
0395 
0396           auto clusterType = comp.first;
0397 
0398           unsigned int clusterID = comp.second;
0399 
0400           const RawCluster *cluster;
0401 
0402           float cluster_energy;
0403           float cluster_Et;
0404 
0405           
0406 
0407           if (clusterType == Jet::SRC::HCAL_TOPO_CLUSTER) {
0408             this_hcal_cluster_mult++;
0409             
0410             cluster = clustersHCAL->getCluster(clusterID);
0411             cluster_energy = cluster->get_energy();
0412             cluster_Et = RawClusterUtility::GetET(*cluster,vertex);
0413 
0414             _fh += cluster_energy;
0415             _fEt_h += cluster_Et;
0416             
0417             m_hcalE.push_back(cluster_energy);
0418             this_jet_hcal_et.push_back(cluster_Et);
0419             m_constE.push_back(cluster_energy);
0420             this_jet_const_et.push_back(cluster_Et);
0421 
0422             RawCluster::TowerConstRange towers = cluster -> get_towers();
0423             RawCluster::TowerConstIterator toweriter;
0424 
0425             for(toweriter = towers.first; toweriter != towers.second; toweriter++) {
0426               float towE = toweriter -> second;
0427               m_ClusterTowE.push_back(towE);
0428             }
0429 
0430           } else if (clusterType == Jet::SRC::ECAL_TOPO_CLUSTER) {
0431             this_emcal_cluster_mult++;
0432 
0433             cluster = clustersEMCAL->getCluster(clusterID);
0434             cluster_energy = cluster->get_energy();
0435             cluster_Et = RawClusterUtility::GetET(*cluster,vertex);
0436 
0437             _fe += cluster_energy;
0438             _fEt_e += cluster_Et;
0439 
0440             m_ecalE.push_back(cluster_energy);
0441             this_jet_emcal_et.push_back(cluster_Et);
0442             m_constE.push_back(cluster_energy);
0443             this_jet_const_et.push_back(cluster_Et);
0444 
0445             RawCluster::TowerConstRange towers = cluster -> get_towers();
0446             RawCluster::TowerConstIterator toweriter;
0447 
0448             for(toweriter = towers.first; toweriter != towers.second; toweriter++) {
0449               float towE = toweriter -> second;
0450               m_ClusterTowE.push_back(towE);
0451             }
0452 
0453 
0454 
0455           }
0456 
0457           
0458 
0459           
0460           
0461 
0462         }
0463 
0464         m_emcal_cluster_mult.push_back(this_emcal_cluster_mult);
0465         m_hcal_cluster_mult.push_back(this_hcal_cluster_mult);
0466         m_cluster_mult.push_back(this_cluster_mult);
0467 
0468         float jet_energy = jet->get_e();
0469         float jet_et = jet->get_et();
0470 
0471         m_fe.push_back(_fe / jet_energy);
0472         m_fh.push_back(_fh / jet_energy);
0473         m_fEt_emcal.push_back(_fEt_e / jet_et);
0474         m_fEt_hcal.push_back(_fEt_h / jet_et);
0475 
0476         m_Et_emcal.push_back(_fEt_e);
0477         m_Et_hcal.push_back(_fEt_h);
0478 
0479         m_constEt.push_back(this_jet_const_et);
0480         m_cluster_hcalEt.push_back(this_jet_hcal_et);
0481         m_cluster_ecalEt.push_back(this_jet_emcal_et);
0482 
0483         
0484 
0485 
0486 
0487 
0488       }
0489       
0490       if(m_doUnsubJet)
0491     {
0492       Jet* unsubjet = new Jetv1();
0493       
0494       float totalPx = 0;
0495       float totalPy = 0;
0496       float totalPz = 0;
0497       float totalE = 0;
0498       int nconst = 0;
0499         
0500       for (auto comp: jet->get_comp_vec())
0501         {
0502           TowerInfo *tower;
0503           nconst++;
0504           unsigned int channel = comp.second;
0505                 
0506           if (comp.first == 15 ||  comp.first == 30)
0507         {
0508           tower = towersIH3->get_tower_at_channel(channel);
0509           if(!tower || !tower_geom){
0510             continue;
0511           }
0512           unsigned int calokey = towersIH3->encode_key(channel);
0513           int ieta = towersIH3->getTowerEtaBin(calokey);
0514           int iphi = towersIH3->getTowerPhiBin(calokey);
0515           const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0516           float UE = background->get_UE(1).at(ieta);
0517           float tower_phi = tower_geom->get_tower_geometry(key)->get_phi();
0518           float tower_eta = tower_geom->get_tower_geometry(key)->get_eta();
0519 
0520           UE = UE * (1 + 2 * background_v2 * cos(2 * (tower_phi - background_Psi2)));
0521           totalE += tower->get_energy() + UE;
0522           double pt = (tower->get_energy() + UE) / cosh(tower_eta);
0523           totalPx += pt * cos(tower_phi);
0524           totalPy += pt * sin(tower_phi);
0525           totalPz += pt * sinh(tower_eta);
0526         }
0527           else if (comp.first == 16 || comp.first == 31)
0528         {
0529           tower = towersOH3->get_tower_at_channel(channel);
0530           if(!tower || !tower_geomOH)
0531             {
0532               continue;
0533             }
0534             
0535           unsigned int calokey = towersOH3->encode_key(channel);
0536           int ieta = towersOH3->getTowerEtaBin(calokey);
0537           int iphi = towersOH3->getTowerPhiBin(calokey);
0538           const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, ieta, iphi);
0539           float UE = background->get_UE(2).at(ieta);
0540           float tower_phi = tower_geomOH->get_tower_geometry(key)->get_phi();
0541           float tower_eta = tower_geomOH->get_tower_geometry(key)->get_eta();
0542             
0543           UE = UE * (1 + 2 * background_v2 * cos(2 * (tower_phi - background_Psi2)));
0544           totalE +=tower->get_energy() + UE;
0545           double pt = (tower->get_energy() + UE) / cosh(tower_eta);
0546           totalPx += pt * cos(tower_phi);
0547           totalPy += pt * sin(tower_phi);
0548           totalPz += pt * sinh(tower_eta);
0549         }
0550           else if (comp.first == 14 || comp.first == 29)
0551         {
0552           tower = towersEM3->get_tower_at_channel(channel);
0553           if(!tower || !tower_geom)
0554             {
0555               continue;
0556             }
0557             
0558           unsigned int calokey = towersEM3->encode_key(channel);
0559           int ieta = towersEM3->getTowerEtaBin(calokey);
0560           int iphi = towersEM3->getTowerPhiBin(calokey);
0561           const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0562           float UE = background->get_UE(0).at(ieta);
0563           float tower_phi = tower_geom->get_tower_geometry(key)->get_phi();
0564           float tower_eta = tower_geom->get_tower_geometry(key)->get_eta();
0565             
0566           UE = UE * (1 + 2 * background_v2 * cos(2 * (tower_phi - background_Psi2)));
0567           totalE +=tower->get_energy() + UE;
0568           double pt = (tower->get_energy() + UE) / cosh(tower_eta);
0569           totalPx += pt * cos(tower_phi);
0570           totalPy += pt * sin(tower_phi);
0571           totalPz += pt * sinh(tower_eta);
0572             
0573         }
0574         }
0575       //get unsubtracted jet
0576       unsubjet->set_px(totalPx);
0577       unsubjet->set_py(totalPy);
0578       unsubjet->set_pz(totalPz);
0579       unsubjet->set_e(totalE);
0580       m_unsub_pt.push_back(unsubjet->get_pt());
0581       m_sub_et.push_back(unsubjet->get_et() - jet->get_et());
0582     }
0583    
0584   
0585       m_nJet++;
0586     }
0587 
0588   //get truth jets
0589   if(m_doTruthJets)
0590     {
0591       m_nTruthJet = 0;
0592       //for (JetMap::Iter iter = jetsMC->begin(); iter != jetsMC->end(); ++iter)
0593       for (auto truthjet : *jetsMC)
0594     {
0595       //Jet* truthjet = iter->second;
0596         
0597       bool eta_cut = (truthjet->get_eta() >= m_etaRange.first) and (truthjet->get_eta() <= m_etaRange.second);
0598       bool pt_cut = (truthjet->get_pt() >= m_ptRange.first) and (truthjet->get_pt() <= m_ptRange.second);
0599       if ((not eta_cut) or (not pt_cut)) continue;
0600       m_truthID.push_back(truthjet->get_id());
0601       m_truthNComponent.push_back(truthjet->size_comp());
0602       m_truthEta.push_back(truthjet->get_eta());
0603       m_truthPhi.push_back(truthjet->get_phi());
0604       m_truthE.push_back(truthjet->get_e());
0605       m_truthPt.push_back(truthjet->get_pt());
0606       m_nTruthJet++;
0607     }
0608     }
0609   
0610   //get seed jets
0611   if(m_doSeeds)
0612     {
0613       for (auto jet : *seedjetsraw)
0614     {
0615       int passesCut = jet->get_property(seedjetsraw->property_index(Jet::PROPERTY::prop_SeedItr));
0616       m_eta_rawseed.push_back(jet->get_eta());
0617       m_phi_rawseed.push_back(jet->get_phi());
0618       m_e_rawseed.push_back(jet->get_e());
0619       m_pt_rawseed.push_back(jet->get_pt());
0620       m_rawseed_cut.push_back(passesCut);
0621     }
0622       
0623       for (auto jet : *seedjetssub)
0624     {
0625       int passesCut = jet->get_property(seedjetssub->property_index(Jet::PROPERTY::prop_SeedItr));
0626       m_eta_subseed.push_back(jet->get_eta());
0627       m_phi_subseed.push_back(jet->get_phi());
0628       m_e_subseed.push_back(jet->get_e());
0629       m_pt_subseed.push_back(jet->get_pt());
0630       m_subseed_cut.push_back(passesCut);
0631     }
0632     }
0633 
0634   Gl1Packet *gl1Packet = findNode::getClass<Gl1Packet>(topNode, "GL1Packet");
0635 
0636   if(gl1Packet)
0637   {
0638     auto scaled_vector = gl1Packet->getScaledVector();
0639     for(int i = 0; i < 32; i++)
0640     {
0641       if((scaled_vector & (int)std::pow(2,i)) != 0)
0642       {
0643         m_triggers.push_back(i);
0644       }
0645     }
0646   }
0647   else
0648   {
0649     std::cout << "gl1Packet not found" << std::endl;
0650   }
0651   //fill the tree
0652   m_T->Fill();
0653 
0654   return Fun4AllReturnCodes::EVENT_OK;
0655 }
0656 
0657 //____________________________________________________________________________..
0658 int JetValidationTC::ResetEvent(PHCompositeNode *topNode)
0659 {
0660 
0661   m_hasJetAboveThreshold = 0;
0662   //std::cout << "JetValidationTC::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
0663   m_id.clear();
0664   m_nComponent.clear();
0665   m_eta.clear();
0666   m_phi.clear();
0667   m_e.clear();
0668   m_pt.clear();
0669   m_zg.clear();
0670   m_rg.clear();
0671   
0672   m_unsub_pt.clear();
0673   m_sub_et.clear();
0674 
0675   m_fe.clear();
0676   m_fh.clear();
0677   m_constE.clear();
0678   m_hcalE.clear();
0679   m_ecalE.clear();
0680   m_fEt_emcal.clear();
0681   m_fEt_hcal.clear();
0682   m_constEt.clear();
0683   m_cluster_hcalEt.clear();
0684   m_cluster_ecalEt.clear();
0685   m_Et_emcal.clear();
0686   m_Et_hcal.clear();
0687   m_ClusterTowE.clear();
0688   m_emcal_cluster_mult.clear();
0689   m_hcal_cluster_mult.clear();
0690   m_cluster_mult.clear();
0691 
0692   cluster_energy.clear();
0693   cluster_eta.clear();
0694   cluster_phi.clear();
0695   cluster_prob.clear();
0696   cluster_chi2.clear();
0697 
0698   pi0cand_pt.clear();
0699   pi0cand_eta.clear();
0700   pi0cand_phi.clear();
0701   pi0cand_mass.clear();
0702   pi0cand_energy.clear();
0703   m_jetPionMult.clear();
0704   m_pionZ.clear();
0705   m_jetPionPt.clear();
0706   m_jetPionEta.clear();
0707   m_jetPionPhi.clear();
0708   m_jetPionMass.clear();
0709   m_jetPionEnergy.clear();
0710   //m_pionMassJetsRemoved.clear();
0711 
0712   m_triggers.clear();
0713 
0714   m_truthID.clear();
0715   m_truthNComponent.clear();
0716   m_truthEta.clear();
0717   m_truthPhi.clear();
0718   m_truthE.clear();
0719   m_truthPt.clear();
0720   m_truthdR.clear();
0721 
0722   m_eta_subseed.clear();
0723   m_phi_subseed.clear();
0724   m_e_subseed.clear();
0725   m_pt_subseed.clear();
0726   m_subseed_cut.clear();
0727 
0728   m_eta_rawseed.clear();
0729   m_phi_rawseed.clear();
0730   m_e_rawseed.clear();
0731   m_pt_rawseed.clear();
0732   m_rawseed_cut.clear();
0733   return Fun4AllReturnCodes::EVENT_OK;
0734 }
0735 
0736 //____________________________________________________________________________..
0737 int JetValidationTC::EndRun(const int runnumber)
0738 {
0739   std::cout << "JetValidationTC::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0740   return Fun4AllReturnCodes::EVENT_OK;
0741 }
0742 
0743 //____________________________________________________________________________..
0744 int JetValidationTC::End(PHCompositeNode *topNode)
0745 {
0746   
0747   std::cout << "JetValidationTC::End: Dumping Histograms" << std::endl;
0748 
0749   //h_pi0M->SetDirectory(nullptr);
0750   //h_npions->SetDirectory(nullptr);
0751 
0752   hm->registerHisto(h_pi0M);
0753   hm->registerHisto(h_npions);
0754   hm->registerHisto(h_eta_phi_clusters);
0755   hm->registerHisto(h_jetPionMult);
0756   hm->registerHisto(h_photonEnergy);
0757 
0758   
0759 
0760   hm->dumpHistos(m_histoFileName.c_str(),"RECREATE");
0761   
0762   std::cout << "JetValidationTC::End - Output to " << m_outputFileName << std::endl;
0763   outfile->cd();
0764   m_T->Write();
0765   //m_T->cd();
0766   outfile->Write();
0767   outfile->Close();
0768   //delete outfile;
0769   //PHTFileServer::get().cd(m_outputFileName);
0770 
0771  
0772   
0773 
0774   
0775   std::cout << "JetValidationTC::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0776   return Fun4AllReturnCodes::EVENT_OK;
0777 }
0778 
0779 //____________________________________________________________________________..
0780 int JetValidationTC::Reset(PHCompositeNode *topNode)
0781 {
0782   std::cout << "JetValidationTC::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0783   return Fun4AllReturnCodes::EVENT_OK;
0784 }
0785 
0786 //____________________________________________________________________________..
0787 void JetValidationTC::Print(const std::string &what) const
0788 {
0789   std::cout << "JetValidationTC::Print(const std::string &what) const Printing info for " << what << std::endl;
0790 }
0791 
0792 
0793 void JetValidationTC::reconstruct_pi0s(PHCompositeNode *topNode) {
0794   
0795   m_npi0 = 0;
0796 
0797   //check whether there are any jets above jet threshold in event
0798 
0799   JetContainer* jets = findNode::getClass<JetContainer>(topNode, m_recoJetName);
0800   if (!jets)
0801   {
0802     std::cout
0803     << "MyJetAnalysis::process_event - Error can not find DST Reco JetContainer node "
0804     << m_recoJetName << std::endl;
0805     exit(-1);
0806   }
0807 
0808   std::vector<float> these_jetsEta;
0809   std::vector<float> these_jetsPhi;
0810 
0811   for (auto jet : *jets) {
0812     if(jet->get_pt() < 1) continue; // to remove noise jets
0813     if(jet->get_pt() < m_ptRange.first) continue; // remove jets outside of pt range
0814 
0815     if (jet->get_pt() > m_pi0_jetThreshold) {
0816       m_hasJetAboveThreshold = 1;
0817       continue;
0818     }
0819 
0820   }
0821   
0822   if (m_hasJetAboveThreshold != 1) {
0823     std::cout << "No jets above threshold. Not reconstructing pi0s" << std::endl;
0824     return;  
0825   }
0826 
0827   nclus = 0;
0828   
0829   //reconstruct pions now that we've found a jet > threshold
0830   RawClusterContainer *clustersEMCAL = findNode::getClass<RawClusterContainer>(topNode, m_clusterType);
0831   //RawClusterContainer *clustersEMCAL = findNode::getClass<RawClusterContainer>(topNode,"CLUSTERINFO_POS_COR_CEMC");
0832   //zvertex
0833   GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0834   if (!vertexmap)
0835   {
0836     std::cout
0837     << "MyJetAnalysis::reconstruct_pi0s - Error can not find global vertex  node "
0838     << std::endl;
0839     exit(-1);
0840   }
0841   if (vertexmap->empty())
0842   {
0843     std::cout
0844     << "MyJetAnalysis::reconstruct_pi0s - global vertex node is empty "
0845     << std::endl;
0846     exit(-1);
0847   }
0848 
0849   GlobalVertex *vtx = vertexmap->begin()->second;
0850   m_zvtx = vtx->get_z();
0851   CLHEP::Hep3Vector vertex;
0852   vertex.set(vtx->get_x(),vtx->get_y(),vtx->get_z());
0853 
0854 
0855   auto clusters = clustersEMCAL->getClusters();
0856   RawClusterContainer::ConstIterator clusIter;
0857 
0858   for (clusIter = clusters.first; clusIter != clusters.second; clusIter++) {
0859     const RawCluster *cluster = clusIter->second;
0860 
0861     if (cluster->get_energy() > m_minClusterE) {
0862       CLHEP::Hep3Vector cluster_vector = RawClusterUtility::GetECoreVec(*cluster, vertex);
0863       /*
0864       cluster_energy.push_back(cluster->get_energy());
0865       cluster_eta.push_back(cluster_vector.pseudoRapidity());
0866       std::cout << "JetValidationTC::reconstruct_pi0: cluster eta = " << cluster_vector.pseudoRapidity() << std::endl;
0867       cluster_phi.push_back(cluster->get_phi());
0868       cluster_prob.push_back(cluster->get_prob());
0869       cluster_chi2.push_back(cluster->get_chi2());*/
0870 
0871       float this_eta = RawClusterUtility::GetPseudorapidity(*cluster,vertex);
0872       float this_phi = RawClusterUtility::GetAzimuthAngle(*cluster,vertex);
0873 
0874       bool clusterIsWithinJetCone = false;
0875 
0876       for(auto jet: *jets) {
0877         float this_jet_eta = jet->get_eta();
0878         float this_jet_phi = jet->get_phi();
0879 
0880         float dR = sqrt(pow((this_eta - this_jet_eta),2) + pow((this_phi - this_jet_phi),2));
0881 
0882         if (dR < m_maxdR) {
0883           clusterIsWithinJetCone = true;
0884           continue;
0885         }
0886 
0887       }
0888 
0889       if (m_removeJetClusters && clusterIsWithinJetCone) {
0890         continue;
0891       }
0892 
0893       nclus++;
0894 
0895       cluster_energy.push_back(cluster->get_energy());
0896       h_photonEnergy->Fill(cluster->get_energy());
0897       
0898       cluster_eta.push_back(this_eta);
0899       cluster_phi.push_back(this_phi);
0900 
0901       h_eta_phi_clusters->Fill(this_eta,this_phi);
0902 
0903       std::cout << "JetValidationTC::reconstruct_pi0: cluster eta = " << RawClusterUtility::GetPseudorapidity(*cluster,vertex) << std::endl;
0904       std::cout << "JetValidationTC::reconstruct_pi0: cluster phi = " << RawClusterUtility::GetAzimuthAngle(*cluster,vertex) << std::endl;
0905       
0906     }
0907     
0908   }
0909 
0910 
0911   for (int i = 0; i < nclus; i++) {
0912     TLorentzVector v1;
0913     v1.SetPtEtaPhiM(cluster_energy[i] / cosh(cluster_eta[i]), cluster_eta[i], cluster_phi[i], 0.0);
0914     
0915     for (int j = i + 1; j < nclus; j++) {
0916       /*float alpha = fabs(cluster_energy[i] - cluster_energy[j]) / (cluster_energy[i] + cluster_energy[j]);
0917 
0918       
0919       if (alpha > 0.5) {
0920         continue;
0921       }*/
0922 
0923       TLorentzVector v2;
0924       v2.SetPtEtaPhiM(cluster_energy[j] / cosh(cluster_eta[j]), cluster_eta[j], cluster_phi[j], 0.0);
0925       
0926       if (v1.DeltaR(v2) < m_mindR) {
0927         continue;
0928       }
0929 
0930       TLorentzVector res;
0931       res = v1 + v2;
0932 
0933       pi0cand_pt.push_back(res.Pt());
0934       pi0cand_eta.push_back(res.Eta());
0935       pi0cand_phi.push_back(res.Phi());
0936       pi0cand_mass.push_back(res.M());
0937       pi0cand_energy.push_back(res.E());
0938       h_pi0M->Fill(res.M());
0939       //m_pionMassJetsRemoved.push_back(res.M()); // get all of the masses and remove those in jets later
0940       std::cout << "JetValidationTC::reconstruct_pi0: pi0 candidate mass = " << res.M() << std::endl;
0941       std::cout << "JetValidationTC::reconstruct_pi0: pi0 candidate energy = " << res.E() << std::endl;
0942       
0943       double mag = res.Mag2();
0944       
0945       std::cout << "JetValidationTC::reconstruct_pi0: pi0 candidate mag2 =  " << mag << std::endl;
0946       m_npi0++;
0947 
0948     }
0949   }
0950 
0951   h_npions->Fill(m_npi0);
0952 
0953 }
0954 
0955 
0956 void JetValidationTC::findJetPions(Jet* jet) {
0957   std::cout << "JetValidationTC::findJetPions: looking for pions in jet with pT = " << jet->get_pt() << std::endl;
0958   int pions_in_this_jet = 0;
0959 
0960   float jetEta = jet->get_eta();
0961   float jetPhi = jet->get_phi();
0962   float jetpT = jet->get_pt();
0963 
0964   for (int pion = 0; pion < m_npi0; pion++){
0965     float dR = sqrt(pow((jetEta - pi0cand_eta[pion]),2) + pow((jetPhi - pi0cand_phi[pion]),2));
0966     if(dR < m_maxdR){
0967       float z = pi0cand_pt[pion] / jetpT;
0968 
0969       m_pionZ.push_back(z);
0970 
0971       m_jetPionPt.push_back(pi0cand_pt[pion]);
0972       m_jetPionEta.push_back(pi0cand_eta[pion]);
0973       m_jetPionPhi.push_back(pi0cand_phi[pion]);
0974       m_jetPionMass.push_back(pi0cand_mass[pion]);
0975       m_jetPionEnergy.push_back(pi0cand_energy[pion]);
0976 
0977       /*
0978       auto pionIdx = std::find(m_pionMassJetsRemoved.begin(),m_pionMassJetsRemoved.end(),pi0cand_mass[pion]);
0979 
0980       if (pionIdx != m_pionMassJetsRemoved.end()) {
0981         std::cout << "Found pion in jet" << std::endl;
0982         m_pionMassJetsRemoved.erase(pionIdx);
0983       } else {
0984         std::cout << "Could not find the pion!" << std::endl;
0985       }*/
0986 
0987 
0988       pions_in_this_jet++;
0989     }
0990   }
0991   if( pions_in_this_jet != 0) {
0992     m_nPionJets ++;
0993   }
0994 
0995   std::cout << "JetValidationTC::findJetPions: Number of pions associated with jet = " << pions_in_this_jet << std::endl;
0996   
0997 
0998   m_jetPionMult.push_back(pions_in_this_jet);
0999   h_jetPionMult->Fill(pions_in_this_jet);
1000 }
[ 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!