TopoclusterJets/src/EMJetVal.cc

0001 // Modified version of JetValidation.cc for the purpose of subjet analysis -- Jennifer James jennifer.l.james@vanderbilt.eduuti
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 and myself
0004 #include <EMJetVal.h>
0005 #include <fun4all/Fun4AllBase.h>
0006 #include <fun4all/Fun4AllReturnCodes.h>
0007 #include <fun4all/PHTFileServer.h>
0008 #include <phool/PHCompositeNode.h>
0009 #include <phool/getClass.h>
0010 #include <jetbase/JetMap.h>
0011 #include <jetbase/JetContainer.h>
0012 #include <jetbase/Jetv2.h>
0013 #include <jetbase/Jetv1.h>
0014 #include <jetbase/Jet.h>
0015 #include <jetbase/JetAlgo.h>
0016 #include <jetbase/FastJetAlgo.h>
0017 #include <jetbase/JetInput.h>
0018 #include <jetbase/TowerJetInput.h>
0019 #include <jetbase/JetMapv1.h>
0020 #include <jetbase/JetContainer.h>
0021 #include <jetbase/JetContainerv1.h>
0022 #include <centrality/CentralityInfo.h>
0023 #include <globalvertex/GlobalVertex.h>
0024 #include <globalvertex/GlobalVertexMapv1.h>
0025 #include <calobase/RawTower.h>
0026 #include <calobase/RawTowerContainer.h>
0027 #include <calobase/RawTowerGeom.h>
0028 #include <calobase/RawTowerGeomContainer.h>
0029 #include <calobase/TowerInfoContainer.h>
0030 #include <calobase/TowerInfo.h>
0031
0032 #include <mbd/MbdOut.h>
0033 #include "fastjet/AreaDefinition.hh"
0034 #include "fastjet/ClusterSequenceArea.hh"
0035 #include "fastjet/Selector.hh"
0036 #include "fastjet/tools/BackgroundEstimatorBase.hh"
0037 #include "fastjet/tools/JetMedianBackgroundEstimator.hh"
0038 #include <fastjet/JetDefinition.hh>
0039 #include "fastjet/ClusterSequence.hh"
0040 #include "fastjet/contrib/SoftDrop.hh" // In external code, this should be fastjet/contrib/SoftDrop.hh
0041 #include <map>
0042 #include <utility>
0043 #include <cstdlib>  // for exit
0044 #include <memory>  // for allocator_traits<>::value_type
0045 #include <Pythia8/Pythia.h> // Include the Pythia header
0046 #include <jetbackground/TowerBackground.h>
0047
0048 #include <TTree.h>
0049 #include <iostream>
0050 #include <fastjet/PseudoJet.hh>
0051 #include <sstream>
0052 #include <iomanip>
0053 #include <cmath>
0054 #include <vector>
0055
0056 using namespace fastjet;
0057
0058 // ROOT, for histogramming.
0059 #include "TH1.h"
0060 #include "TH2.h"
0061
0062 // ROOT, for interactive graphics.
0063 #include "TVirtualPad.h"
0064 #include "TApplication.h"
0065
0066 // ROOT, for saving file.
0067 #include "TFile.h"
0068
0069
0070
0071 using std::cout;
0072 using std::endl;
0073
0074 //____________________________________________________________________________..
0075 EMJetVal::EMJetVal(const std::string& recojetname, const std::string& truthjetname, const std::string& outputfilename)
0076   :SubsysReco("EMJetVal_" + recojetname + "_" + truthjetname)
0077   , m_recoJetName(recojetname)
0078   , m_vtxZ_cut(10.0)
0079   , m_truthJetName(truthjetname)
0080   , m_outputFileName(outputfilename)
0081   , m_etaRange(-1, 1)
0082   , m_ptRange(5, 100)
0083   , m_doTruthJets(0)
0084   , m_doSeeds(0)
0085   , m_doUnsubJet(0)
0086   , m_T(nullptr)
0087   , m_event(-1)
0088   , m_nTruthJet(-1)
0089   , m_nJet(-1)
0090   , m_id()
0091   , m_nComponent()
0092   , m_eta()
0093   , m_phi()
0094   , m_e()
0095   , m_pt()
0096   , m_sub_et()
0097   , m_truthID()
0098   , m_truthNComponent()
0099   , m_truthEta()
0100   , m_truthPhi()
0101   , m_truthE()
0102   , m_truthPt()
0103   , m_eta_rawseed()
0104   , m_phi_rawseed()
0105   , m_pt_rawseed()
0106   , m_e_rawseed()
0107   , m_rawseed_cut()
0108   , m_eta_subseed()
0109   , m_phi_subseed()
0110   , m_pt_subseed()
0111   , m_e_subseed()
0112   , m_subseed_cut()
0113 {
0114   std::cout << "EMJetVal::EMJetVal(const std::string &name) Calling ctor" << std::endl;
0115 }
0116
0117
0118 //____________________________________________________________________________..
0119 EMJetVal::~EMJetVal()
0120 {
0121   std::cout << "EMJetVal::~EMJetVal() Calling dtor" << std::endl;
0122 }
0123
0124 //____________________________________________________________________________..
0125 int EMJetVal::Init(PHCompositeNode *topNode)
0126 {
0127   std::cout << "EMJetVal::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0128   //PHTFileServer::get().open(m_outputFileName, "RECREATE");
0129   std::cout << "EMJetVal::Init - Output to " << m_outputFileName << std::endl;
0130   //Analysis hists
0131   outFile = new TFile(m_outputFileName.c_str(), "RECREATE");
0132   _h_R04_z_sj_10_20= new TH1F("R04_z_sj_10_20","z_sj in subjets 1 & 2", 10, 0, 0.5);
0133   _h_R04_theta_sj_10_20= new TH1F("R04_theta_sj_10_20","theta_sj in subjets 1 & 2", 10, 0, 0.5);
0134   // softdrop hists
0135    _h_R04_z_g_10_20= new TH1F("R04_z_g_10_20","z_g in subjets 1 & 2", 10, 0, 0.5);
0136    _h_R04_theta_g_10_20= new TH1F("R04_theta_g_10_20","theta_g in subjets 1 & 2", 10, 0, 0.5);
0137   // multiplicity hists
0138   _hmult_R04= new TH1F("mult_R04","total number of constituents inside R=0.4 jets", 30, 0, 30);
0139   _hmult_R04_pT_10_20GeV= new TH1F("mult_R04_pT_10_20GeV","total number of constituents inside R=0.4 jets with 10 < p_{T} < 20", 30, 0, 30);
0140   _hjetpT_R04 = new TH1F("jetpT_R04","jet transverse momentum for R=0.4 jets", 100, 0, 100);
0141   _hjeteta_R04 = new TH1F("jeteta_R04","jet pseudorapidity for R=0.4 jets", 50, -0.6, 0.6);
0142   // corellation hists
0143   correlation_theta_10_20 = new TH2D("correlation_theta_10_20", "Correlation Plot 10 < p_{T} < 20 [GeV/c]; R_{g}; #theta_{sj}", 20, 0, 0.5, 20, 0, 0.5);
0144   correlation_z_10_20 = new TH2D("correlation_z_10_20", "Correlation Plot; z_{g}; z_{sj}", 20, 0, 0.5, 20, 0, 0.5);
0145   // configure Tree
0146
0147   //PHTFileServer::get().open("hist_jets.root", "RECREATE");
0148
0149   m_T = new TTree("T", "MyJetAnalysis Tree");
0150   m_T->Branch("m_event", &m_event, "event/I");
0151   m_T->Branch("nJet", &m_nJet, "nJet/I");
0152   m_T->Branch("cent", &m_centrality);
0153   m_T->Branch("b", &m_impactparam);
0154   m_T->Branch("id", &m_id);
0155   m_T->Branch("nComponent", &m_nComponent);
0156
0157   m_T->Branch("eta", &m_eta);
0158   m_T->Branch("phi", &m_phi);
0159   m_T->Branch("e", &m_e);
0160   m_T->Branch("pt", &m_pt);
0161   if(m_doUnsubJet)
0162     {
0163       m_T->Branch("pt_unsub", &m_unsub_pt);
0164       m_T->Branch("subtracted_et", &m_sub_et);
0165     }
0166   if(m_doTruthJets){
0167     m_T->Branch("nTruthJet", &m_nTruthJet);
0168     m_T->Branch("truthID", &m_truthID);
0169     m_T->Branch("truthNComponent", &m_truthNComponent);
0170     m_T->Branch("truthEta", &m_truthEta);
0171     m_T->Branch("truthPhi", &m_truthPhi);
0172     m_T->Branch("truthE", &m_truthE);
0173     m_T->Branch("truthPt", &m_truthPt);
0174   }
0175
0176   if(m_doSeeds){
0177     m_T->Branch("rawseedEta", &m_eta_rawseed);
0178     m_T->Branch("rawseedPhi", &m_phi_rawseed);
0179     m_T->Branch("rawseedPt", &m_pt_rawseed);
0180     m_T->Branch("rawseedE", &m_e_rawseed);
0181     m_T->Branch("rawseedCut", &m_rawseed_cut);
0182     m_T->Branch("subseedEta", &m_eta_subseed);
0183     m_T->Branch("subseedPhi", &m_phi_subseed);
0184     m_T->Branch("subseedPt", &m_pt_subseed);
0185     m_T->Branch("subseedE", &m_e_subseed);
0186     m_T->Branch("subseedCut", &m_subseed_cut);
0187   }
0188
0189   //  std::cout << "finished declaring histos" << std::endl;
0190
0191   return Fun4AllReturnCodes::EVENT_OK;
0192 }
0193 //____________________________________________________________________________..
0194
0195 int EMJetVal::retrieveEvent(const fastjet::PseudoJet& jet) {
0196   // Add the PseudoJet to the event vector
0197   eventVector.push_back(jet);
0198
0199   return 0; // Return 0 for success
0200 }
0201
0202 //____________________________________________________________________________..
0203 int EMJetVal::InitRun(PHCompositeNode *topNode)
0204 {
0205   std::cout << "EMJetVal::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0206   return Fun4AllReturnCodes::EVENT_OK;
0207 }
0208
0209 //____________________________________________________________________________..
0210 int EMJetVal::process_event(PHCompositeNode *topNode)
0211 {
0212   // std::cout << "EMJetVal::process_event(PHCompositeNode *topNode) Processing Event " << m_event << std::endl;
0213    ++m_event;
0214
0215   // interface to reco jets
0216   JetContainer* jets = findNode::getClass<JetContainer>(topNode, m_recoJetName);
0217   if (!jets)
0218     {
0219       std::cout
0220     << "MyJetAnalysis::process_event - Error can not find DST Reco JetContainer node "
0221     << m_recoJetName << std::endl;
0222       exit(-1);
0223     }
0224
0225   //interface to truth jets
0226   JetMap* jetsMC = findNode::getClass<JetMap>(topNode, m_truthJetName);
0227   if (!jetsMC && m_doTruthJets)
0228     {
0229       std::cout
0230     << "MyJetAnalysis::process_event - Error can not find DST Truth JetMap node "
0231     << m_truthJetName << std::endl;
0232       exit(-1);
0233     }
0234
0235
0236   // interface to jet seeds
0237   JetContainer* seedjetsraw = findNode::getClass<JetContainer>(topNode, "AntiKt_TowerInfo_HIRecoSeedsRaw_r02");
0238   if (!seedjetsraw && m_doSeeds)
0239     {
0240       std::cout
0241     << "MyJetAnalysis::process_event - Error can not find DST raw seed jets "
0242     << std::endl;
0243       exit(-1);
0244     }
0245
0246   JetContainer* seedjetssub = findNode::getClass<JetContainer>(topNode, "AntiKt_TowerInfo_HIRecoSeedsSub_r02");
0247   if (!seedjetssub && m_doSeeds)
0248     {
0249       std::cout
0250     << "MyJetAnalysis::process_event - Error can not find DST subtracted seed jets "
0251     << std::endl;
0252       exit(-1);
0253     }
0254
0255   //centrality
0256   CentralityInfo* cent_node = findNode::getClass<CentralityInfo>(topNode, "CentralityInfo");
0257   if (!cent_node)
0258     {
0259       std::cout
0260     << "MyJetAnalysis::process_event - Error can not find centrality node "
0261     << std::endl;
0262       exit(-1);
0263     }
0264
0265   //calorimeter towers
0266   TowerInfoContainer *towersEM3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC_RETOWER_SUB1");
0267   TowerInfoContainer *towersIH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN_SUB1");
0268   TowerInfoContainer *towersOH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT_SUB1");
0269   RawTowerGeomContainer *tower_geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0270   RawTowerGeomContainer *tower_geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0271   if(!towersEM3 || !towersIH3 || !towersOH3){
0272     std::cout
0273       <<"MyJetAnalysis::process_event - Error can not find raw tower node "
0274       << std::endl;
0275     exit(-1);
0276   }
0277
0278   if(!tower_geom || !tower_geomOH){
0279     std::cout
0280       <<"MyJetAnalysis::process_event - Error can not find raw tower geometry "
0281       << std::endl;
0282     exit(-1);
0283   }
0284
0285   //underlying event
0286   TowerBackground *background = findNode::getClass<TowerBackground>(topNode, "TowerInfoBackground_Sub2");
0287   if(!background){
0288     //  std::cout<<"Can't get background. Exiting"<<std::endl;
0289     return Fun4AllReturnCodes::EVENT_OK;
0290   }
0291
0292   //get the event centrality/impact parameter from HIJING
0293   m_centrality =  cent_node->get_centile(CentralityInfo::PROP::bimp);
0294   m_impactparam =  cent_node->get_quantity(CentralityInfo::PROP::bimp);
0295
0296   //get reco jets
0297   m_nJet = 0;
0298   float background_v2 = 0;
0299   float background_Psi2 = 0;
0300   if(m_doUnsubJet)
0301     {
0302       background_v2 = background->get_v2();
0303       background_Psi2 = background->get_Psi2();
0304     }
0305   for (auto jet : *jets)// jets //JetMap::Iter iter = jets->begin(); iter != jets->end(); ++iter)
0306     {
0307
0308       // std::cout << "working on original jet " << jet->get_id() << " out of " << jets->size() << std::endl;
0309
0310       if(jet->get_pt() < 1) continue; // to remove noise jets
0311
0312
0313       m_id.push_back(jet->get_id());
0314       m_nComponent.push_back(jet->size_comp());
0315       m_eta.push_back(jet->get_eta());
0316       m_phi.push_back(jet->get_phi());
0317       m_e.push_back(jet->get_e());
0318       m_pt.push_back(jet->get_pt());
0319
0320       std::vector<PseudoJet> particles;
0321       particles.clear();
0322
0323       float totalPx = 0;
0324       float totalPy = 0;
0325       float totalPz = 0;
0326       float totalE = 0;
0327       int nconst = 0;
0328
0329       for (auto comp : jet->get_comp_vec()) //Jet::ConstIter comp = jet->begin_comp(); comp != jet->end_comp(); ++comp)
0330         {
0331           TowerInfo *tower;
0332           nconst++;
0333           unsigned int channel = comp.second;
0334
0335           if (comp.first == 15 || comp.first == 30)
0336         {
0337           tower = towersIH3->get_tower_at_channel(channel);
0338           if(!tower || !tower_geom){
0339             continue;
0340           }
0341           unsigned int calokey = towersIH3->encode_key(channel);
0342           int ieta = towersIH3->getTowerEtaBin(calokey);
0343           int iphi = towersIH3->getTowerPhiBin(calokey);
0344           const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0345           float UE = background->get_UE(1).at(ieta);
0346           float tower_phi = tower_geom->get_tower_geometry(key)->get_phi();
0347           float tower_eta = tower_geom->get_tower_geometry(key)->get_eta();
0348
0349           UE = UE * (1 + 2 * background_v2 * cos(2 * (tower_phi - background_Psi2)));
0350           totalE += tower->get_energy() + UE;
0351           double pt = (tower->get_energy() + UE) / cosh(tower_eta);
0352           totalPx += pt * cos(tower_phi);
0353           totalPy += pt * sin(tower_phi);
0354           totalPz += pt * sinh(tower_eta);
0355         }
0356           else if (comp.first == 16 || comp.first == 31)
0357         {
0358           tower = towersOH3->get_tower_at_channel(channel);
0359           if(!tower || !tower_geomOH)
0360             {
0361               continue;
0362             }
0363
0364           unsigned int calokey = towersOH3->encode_key(channel);
0365           int ieta = towersOH3->getTowerEtaBin(calokey);
0366           int iphi = towersOH3->getTowerPhiBin(calokey);
0367           const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, ieta, iphi);
0368           float UE = background->get_UE(2).at(ieta);
0369           float tower_phi = tower_geomOH->get_tower_geometry(key)->get_phi();
0370           float tower_eta = tower_geomOH->get_tower_geometry(key)->get_eta();
0371
0372           UE = UE * (1 + 2 * background_v2 * cos(2 * (tower_phi - background_Psi2)));
0373           totalE +=tower->get_energy() + UE;
0374           double pt = (tower->get_energy() + UE) / cosh(tower_eta);
0375           totalPx += pt * cos(tower_phi);
0376           totalPy += pt * sin(tower_phi);
0377           totalPz += pt * sinh(tower_eta);
0378         }
0379           else if (comp.first == 14 || comp.first == 29)
0380         {
0381           tower = towersEM3->get_tower_at_channel(channel);
0382           if(!tower || !tower_geom)
0383             {
0384               continue;
0385             }
0386
0387           unsigned int calokey = towersEM3->encode_key(channel);
0388           int ieta = towersEM3->getTowerEtaBin(calokey);
0389           int iphi = towersEM3->getTowerPhiBin(calokey);
0390           const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0391           float UE = background->get_UE(0).at(ieta);
0392           float tower_phi = tower_geom->get_tower_geometry(key)->get_phi();
0393           float tower_eta = tower_geom->get_tower_geometry(key)->get_eta();
0394
0395
0396           UE = UE * (1 + 2 * background_v2 * cos(2 * (tower_phi - background_Psi2)));
0397           double ETmp = tower->get_energy();
0398           double pt = tower->get_energy() / cosh(tower_eta);
0399           if(m_doUnsubJet){
0400             ETmp += UE;
0401             pt = ETmp / cosh(tower_eta);
0402           }
0403           double pxTmp = pt * cos(tower_phi);
0404           double pyTmp = pt * sin(tower_phi);
0405           double pzTmp = pt * sinh(tower_eta);
0406
0407           totalE += ETmp;
0408           totalPx += pxTmp;
0409           totalPy += pyTmp;
0410           totalPz += pzTmp;
0411
0412
0413           particles.push_back( PseudoJet( pxTmp, pyTmp, pzTmp, ETmp) );
0414
0415         }//end if over sources
0416
0417         }//end loop over constituents
0418       //Insert jenn's analysis
0419           double radius[5] = {0.05, 0.1, 0.2, 0.4, 0.6}; // jet radius
0420           double pseudorapidity = -999.; // pseudorapidity
0421           double theta_sj = -1.; // delta radius (value describes an unachievable value)
0422           double z_sj = -1.; // delta radius (value describes an unachievable value)
0423
0424           //  std::cout << "passed here -397 -my jet defs" << std::endl
0425           // Set up FastJet jet finders and modified mass-drop tagger.
0426           JetDefinition jetDefAKT_R01( antikt_algorithm, radius[1]);
0427           JetDefinition jetDefAKT_R04( antikt_algorithm, radius[3]);
0428           JetDefinition jetDefCA(cambridge_algorithm, radius[3]);
0429
0430           //  vector<PseudoJet> particles;
0431
0432           // Run Fastjet anti-kT algorithm and sort jets in pT order.
0433           ClusterSequence clustSeq_R04( particles, jetDefAKT_R04 );
0434           std::vector<PseudoJet> sortedJets_R04 = sorted_by_pt( clustSeq_R04.inclusive_jets() );
0435
0436           //! create a loop to run over the jets -
0437
0438           for (int j = 0; j < (int)sortedJets_R04.size(); j++) {
0439
0440             //      std::cout << "working on reclustered jet " << j << " of " << sortedJets_R04.size() << std::endl;
0441
0442             PseudoJet jet_reco = sortedJets_R04.at(j);
0443             if(fabs(jet_reco.eta()) > 0.6)
0444               continue;
0445
0446             // std::cout << "jet within eta of 0.6" << std::endl;
0447
0448             ClusterSequence clustSeq_R01_con(jet_reco.constituents() , jetDefAKT_R01 );
0449             //  std::cout << "made R=0.1 cluster sequence" << std::endl;
0450             std:: vector<PseudoJet> sortedJets_R01_con = sorted_by_pt( clustSeq_R01_con.inclusive_jets() );
0451             //  std::cout << "ran R=0.1 and sorted by pT" << std::endl;
0452             // std::cout << "number of subjets: " << sortedJets_R01_con.size() << std::endl;
0453             // std::cout << "passed here -476 -pseudojet" << std::endl
0454             if (sortedJets_R01_con.size() < 2){
0455               //  std::cout << "not enough subjets" << std::endl;
0456               continue;
0457             }
0458
0459             //std:: cout << "at least 2 subjets" << std::endl;
0460
0461             PseudoJet sj1 = sortedJets_R01_con.at(0);
0462             PseudoJet sj2 = sortedJets_R01_con.at(1);
0463
0464             // std::cout << "sj1 pt=" << sj1.pt() << std::endl;
0465             // std::cout << "sj2 pt=" << sj2.pt() << std::endl;
0466             if (sj1.pt() < 3 || sj2.pt() < 3 )
0467               continue;
0468             // std::cout << "sj1 pt=" << sj1.pt() << std::endl;
0469             // std::cout << "sj2 pt=" << sj2.pt() << std::endl;
0470             // std::cout << "both are above 3" << std::endl;
0471
0472             theta_sj = sj1.delta_R(sj2);
0473             z_sj = sj2.pt()/(sj2.pt()+sj1.pt());
0474
0475             // std::cout << "theta_sj = " << theta_sj << "   z_sj = " << z_sj << std::endl;
0476
0477             // 10 to 20 pT
0478             if (jet_reco.pt() > 10 && jet_reco.pt() < 20 ){
0479               // std::cout<<"sorted jets at "<<j<<" the pT = " << jet.pt()<<endl;
0480               // std::cout << "jet pt >10 & <20: " << jet_reco.pt() << std::endl;
0481               _hjetpT_R04->Fill(jet_reco.perp());
0482               pseudorapidity = jet_reco.eta();
0483               _hjeteta_R04->Fill(pseudorapidity);
0484               _hmult_R04->Fill(jet_reco.constituents().size());
0485
0486               ClusterSequence clustSeqCA(jet_reco.constituents(), jetDefCA);
0487               std::vector<PseudoJet> cambridgeJets = sorted_by_pt(clustSeqCA.inclusive_jets());
0488
0489               // std::cout << "have CA sort jets: " << cambridgeJets.size() << std::endl;
0490
0491               // SoftDrop parameters
0492               double z_cut = 0.10;
0493               double beta = 0.0;
0494               contrib::SoftDrop sd(beta, z_cut);
0495               //! get subjets
0496               if (!isnan(theta_sj) && !isnan(z_sj) && !isinf(theta_sj) && !isinf(z_sj)){
0497               _h_R04_z_sj_10_20->Fill(z_sj);
0498               _h_R04_theta_sj_10_20->Fill(theta_sj);
0499             }
0500
0501               // std::cout << "filled some histos" << std::endl;
0502               // Apply SoftDrop to the jet
0503               PseudoJet sd_jet = sd(jet_reco);
0504               // std::cout << "ran sd" << std::endl;
0505               if (sd_jet == 0)
0506
0507             continue;
0508               // std::cout << "sd jet exists" << std::endl;
0509                double delta_R_subjets = sd_jet.structure_of<contrib::SoftDrop>().delta_R();
0510                double z_subjets = sd_jet.structure_of<contrib::SoftDrop>().symmetry();
0511
0512               _h_R04_z_g_10_20->Fill(z_subjets);
0513               _h_R04_theta_g_10_20->Fill(delta_R_subjets);
0514
0515               correlation_theta_10_20->Fill(delta_R_subjets, theta_sj);
0516               correlation_z_10_20->Fill(z_subjets, z_sj);
0517
0518               // SoftDrop failed, handle the case as needed
0519               // e.g., skip this jet or perform alternative analysis
0520             } else {
0521               _hmult_R04_pT_10_20GeV->Fill(jet_reco.constituents().size());
0522             }
0523
0524             //! jet loop
0525             // Rjet = 0.4 End
0526
0527             //filled nEvent in histogram
0528             _hmult_R04->Fill(m_nJet);
0529
0530             //          std::cout << "iZ = " << nEvent << std::endl;
0531             // std::cout << "finished jet " << j << std::endl;
0532
0533           }//! event loop ends for pT
0534
0535
0536           //std::cout << "finished loop over reclustered jets" << std::endl;
0537
0538           // End of event loop. Statistics. Histograms. Done.
0539
0540
0541           m_nJet++;
0542           // std::cout << "m_nJet: " << m_nJet << std::endl;
0543     }
0544
0545   // std::cout << "finised loop over original jets" << std::endl;
0546
0547 //get truth jets
0548 if(m_doTruthJets)
0549   {
0550     m_nTruthJet = 0;
0551     for (JetMap::Iter iter = jetsMC->begin(); iter != jetsMC->end(); ++iter)
0552       {
0553     Jet* truthjet = iter->second;
0554
0555     bool eta_cut = (truthjet->get_eta() >= m_etaRange.first) and (truthjet->get_eta() <= m_etaRange.second);
0556     bool pt_cut = (truthjet->get_pt() >= m_ptRange.first) and (truthjet->get_pt() <= m_ptRange.second);
0557     if ((not eta_cut) or (not pt_cut)) continue;
0558     m_truthID.push_back(truthjet->get_id());
0559     m_truthNComponent.push_back(truthjet->size_comp());
0560     m_truthEta.push_back(truthjet->get_eta());
0561     m_truthPhi.push_back(truthjet->get_phi());
0562     m_truthE.push_back(truthjet->get_e());
0563     m_truthPt.push_back(truthjet->get_pt());
0564     m_nTruthJet++;
0565       }
0566   }
0567   //get seed jets
0568   if(m_doSeeds)
0569     {
0570       for (auto jet : *seedjetsraw)
0571     {
0572       int passesCut = jet->get_property(Jet::PROPERTY::prop_SeedItr);
0573       m_eta_rawseed.push_back(jet->get_eta());
0574       m_phi_rawseed.push_back(jet->get_phi());
0575       m_e_rawseed.push_back(jet->get_e());
0576       m_pt_rawseed.push_back(jet->get_pt());
0577       m_rawseed_cut.push_back(passesCut);
0578     }
0579
0580       for (auto jet : *seedjetssub) //JetMap::Iter iter = seedjetssub->begin(); iter != seedjetssub->end(); ++iter)
0581     {
0582       int passesCut = jet->get_property(Jet::PROPERTY::prop_SeedItr);
0583       m_eta_subseed.push_back(jet->get_eta());
0584       m_phi_subseed.push_back(jet->get_phi());
0585       m_e_subseed.push_back(jet->get_e());
0586       m_pt_subseed.push_back(jet->get_pt());
0587       m_subseed_cut.push_back(passesCut);
0588     }
0589     }
0590
0591   //fill the tree
0592
0593   m_T->Fill();
0594   //  std::cout << "filled TTree" << std::endl;
0595
0596   return Fun4AllReturnCodes::EVENT_OK;
0597 }
0598
0599
0600     //____________________________________________________________________________..
0601   int EMJetVal::ResetEvent(PHCompositeNode *topNode)
0602   {
0603     //std::cout << "EMJetVal::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
0604     m_id.clear();
0605     m_nComponent.clear();
0606     m_eta.clear();
0607     m_phi.clear();
0608     m_e.clear();
0609     m_pt.clear();
0610     m_unsub_pt.clear();
0611     m_sub_et.clear();
0612
0613     m_truthID.clear();
0614     m_truthNComponent.clear();
0615     m_truthEta.clear();
0616     m_truthPhi.clear();
0617     m_truthE.clear();
0618     m_truthPt.clear();
0619     m_truthdR.clear();
0620
0621     m_eta_subseed.clear();
0622     m_phi_subseed.clear();
0623     m_e_subseed.clear();
0624     m_pt_subseed.clear();
0625     m_subseed_cut.clear();
0626
0627     m_eta_rawseed.clear();
0628     m_phi_rawseed.clear();
0629     m_e_rawseed.clear();
0630     m_pt_rawseed.clear();
0631     m_rawseed_cut.clear();
0632   return Fun4AllReturnCodes::EVENT_OK;
0633
0634   }
0635   //____________________________________________________________________________..
0636   int EMJetVal::EndRun(const int runnumber)
0637   {
0638     std::cout << "EMJetVal::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0639     return Fun4AllReturnCodes::EVENT_OK;
0640   }
0641
0642   //____________________________________________________________________________..
0643   int EMJetVal::End(PHCompositeNode *topNode)
0644   {
0645
0646     outFile->cd();
0647     outFile->Write();
0648     outFile->Close();
0649
0650     std::cout << "EMJetVal::End - Output to " << m_outputFileName << std::endl;
0651
0652
0653     // PHTFileServer::get().cd(m_outputFileName);
0654
0655     //m_T->Write();
0656     std::cout << "EMJetVal::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0657     return Fun4AllReturnCodes::EVENT_OK;
0658   }
0659
0660   //____________________________________________________________________________..
0661   int EMJetVal::Reset(PHCompositeNode *topNode)
0662   {
0663     std::cout << "EMJetVal::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0664     return Fun4AllReturnCodes::EVENT_OK;
0665   }
0666
0667   //____________________________________________________________________________..
0668   void EMJetVal::Print(const std::string &what) const
0669   {
0670      std::cout << "EMJetVal::Print(const std::string &what) const Printing info for " << what << std::endl;
0671   }
0672
0673