Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:12:46 

0001 #include "SoftLeptonTaggingTruth.h"
0002 
0003 #include <fun4all/SubsysReco.h>
0004 #include <fun4all/Fun4AllServer.h>
0005 #include <fun4all/PHTFileServer.h>
0006 #include <phool/PHCompositeNode.h>
0007 #include <fun4all/Fun4AllReturnCodes.h>
0008 #include <phool/getClass.h>
0009 
0010 #include <phool/PHCompositeNode.h>
0011 
0012 #include <g4main/PHG4TruthInfoContainer.h>
0013 #include <g4main/PHG4VtxPoint.h>
0014 #include <g4main/PHG4Particle.h>
0015 #include <g4main/PHG4HitDefs.h>
0016 
0017 #include <g4eval/JetEvalStack.h>
0018 #include <g4eval/JetTruthEval.h>
0019 
0020 #include <TString.h>
0021 #include <TFile.h>
0022 #include <TMath.h>
0023 #include <TH1F.h>
0024 #include <TH2F.h>
0025 #include <TVector3.h>
0026 #include <TLorentzVector.h>
0027 #include <TAxis.h>
0028 #include <TLine.h>
0029 #include <TDatabasePDG.h>
0030 
0031 #include <exception>
0032 #include <stdexcept>
0033 #include <iostream>
0034 #include <vector>
0035 #include <set>
0036 #include <algorithm>
0037 #include <cassert>
0038 #include <cmath>
0039 
0040 using namespace std;
0041 
0042 SoftLeptonTaggingTruth::SoftLeptonTaggingTruth(const std::string & truth_jet,
0043     enu_flags flags) :
0044     SubsysReco("SoftLeptonTaggingTruth_" + truth_jet), //
0045     _jetevalstacks(), //
0046     _truth_jet(truth_jet), _reco_jets(), _truth_container(NULL), _flags(flags), //
0047     eta_range(-1, 1), //
0048     _jet_match_dR(.7), _jet_match_dca(.2), _jet_match_E_Ratio(.0)
0049 {
0050 
0051 }
0052 
0053 SoftLeptonTaggingTruth::~SoftLeptonTaggingTruth()
0054 {
0055 }
0056 
0057 int
0058 SoftLeptonTaggingTruth::InitRun(PHCompositeNode *topNode)
0059 {
0060   _truth_container = findNode::getClass<PHG4TruthInfoContainer>(topNode,
0061       "G4TruthInfo");
0062   if (!_truth_container)
0063     {
0064       cout << "SoftLeptonTaggingTruth::InitRun - Fatal Error - "
0065           << "unable to find DST node " << "G4TruthInfo" << endl;
0066       assert(_truth_container);
0067     }
0068 
0069   if (flag(kProcessTruthSpectrum))
0070     {
0071       if (not _jettrutheval)
0072         _jettrutheval = shared_ptr < JetTruthEval
0073             > (new JetTruthEval(topNode, _truth_jet));
0074 
0075       assert(_jettrutheval);
0076       _jettrutheval->set_strict(true);
0077       _jettrutheval->set_verbosity(verbosity + 1);
0078     }
0079 
0080   return Fun4AllReturnCodes::EVENT_OK;
0081 }
0082 
0083 int
0084 SoftLeptonTaggingTruth::End(PHCompositeNode *topNode)
0085 {
0086 
0087   return Fun4AllReturnCodes::EVENT_OK;
0088 }
0089 
0090 int
0091 SoftLeptonTaggingTruth::Init(PHCompositeNode *topNode)
0092 {
0093 
0094   Fun4AllHistoManager *hm = getHistoManager();
0095   assert(hm);
0096 
0097   if (flag(kProcessTruthSpectrum))
0098     {
0099       if (verbosity >= 1)
0100         cout << "SoftLeptonTaggingTruth::Init - Process TruthSpectrum "
0101             << _truth_jet << endl;
0102       Init_Spectrum(topNode, _truth_jet);
0103     }
0104 
0105   return Fun4AllReturnCodes::EVENT_OK;
0106 }
0107 
0108 int
0109 SoftLeptonTaggingTruth::process_event(PHCompositeNode *topNode)
0110 {
0111 
0112   if (verbosity > 2)
0113     cout << "SoftLeptonTaggingTruth::process_event() entered" << endl;
0114 
0115   for (jetevalstacks_map::iterator it_jetevalstack = _jetevalstacks.begin();
0116       it_jetevalstack != _jetevalstacks.end(); ++it_jetevalstack)
0117     {
0118       assert(it_jetevalstack->second);
0119       it_jetevalstack->second->next_event(topNode);
0120     }
0121   if (_jettrutheval)
0122     _jettrutheval->next_event(topNode);
0123 
0124   if (flag(kProcessTruthSpectrum))
0125     {
0126       if (verbosity >= 1)
0127         cout << "SoftLeptonTaggingTruth::process_event - Process TruthSpectrum "
0128             << _truth_jet << endl;
0129       process_Spectrum(topNode, _truth_jet, false);
0130     }
0131 
0132   return Fun4AllReturnCodes::EVENT_OK;
0133 }
0134 
0135 //! set eta range
0136 void
0137 SoftLeptonTaggingTruth::set_eta_range(double low, double high)
0138 {
0139   if (low > high)
0140     swap(low, high);
0141   assert(low < high);
0142   // eliminate zero range
0143 
0144   eta_range.first = low;
0145   eta_range.second = high;
0146 }
0147 
0148 //! string description of eta range
0149 //! @return TString as ROOT likes
0150 TString
0151 SoftLeptonTaggingTruth::get_eta_range_str(const char * eta_name) const
0152 {
0153   assert(eta_name);
0154   return TString(
0155       Form("%.1f < %s < %.1f", eta_range.first, eta_name, eta_range.second));
0156 }
0157 
0158 //! acceptance cut on jet object
0159 bool
0160 SoftLeptonTaggingTruth::jet_acceptance_cut(const Jet * jet) const
0161 {
0162   assert(jet);
0163   bool eta_cut = (jet->get_eta() >= eta_range.first)
0164       and (jet->get_eta() <= eta_range.second);
0165   return eta_cut;
0166 }
0167 
0168 string
0169 SoftLeptonTaggingTruth::get_histo_prefix(const std::string & src_jet_name,
0170     const std::string & reco_jet_name)
0171 {
0172   std::string histo_prefix = "h_QAG4SimJet_";
0173 
0174   if (src_jet_name.length() > 0)
0175     {
0176       histo_prefix += src_jet_name;
0177       histo_prefix += "_";
0178     }
0179   if (reco_jet_name.length() > 0)
0180     {
0181       histo_prefix += reco_jet_name;
0182       histo_prefix += "_";
0183     }
0184 
0185   return histo_prefix;
0186 }
0187 
0188 int
0189 SoftLeptonTaggingTruth::Init_Spectrum(PHCompositeNode *topNode,
0190     const std::string & jet_name)
0191 {
0192 
0193   Fun4AllHistoManager *hm = getHistoManager();
0194   assert(hm);
0195 
0196   // normalization plot with counts
0197   const int norm_size = 3;
0198 
0199   TH1D * h_norm = new TH1D(
0200       TString(get_histo_prefix(jet_name)) + "Normalization",
0201       " Normalization;Item;Count", norm_size, .5, norm_size + .5);
0202   int i = 1;
0203 
0204   h_norm->GetXaxis()->SetBinLabel(i++, "Event");
0205   h_norm->GetXaxis()->SetBinLabel(i++, "Inclusive Jets");
0206   h_norm->GetXaxis()->SetBinLabel(i++, "Leading Jets");
0207   assert(norm_size >= i - 1);
0208   h_norm->GetXaxis()->LabelsOption("v");
0209   hm->registerHisto(h_norm);
0210 
0211   hm->registerHisto(
0212       new TH1F(
0213           //
0214           TString(get_histo_prefix(jet_name)) + "Leading_Et", //
0215           TString(jet_name) + " leading jet Et, " + get_eta_range_str()
0216               + ";E_{T} (GeV)", 100, 0, 100));
0217   hm->registerHisto(
0218       new TH1F(
0219           //
0220           TString(get_histo_prefix(jet_name)) + "Leading_eta", //
0221           TString(jet_name) + " leading jet #eta, " + get_eta_range_str()
0222               + ";#eta", 50, -1, 1));
0223   hm->registerHisto(
0224       new TH1F(
0225           //
0226           TString(get_histo_prefix(jet_name)) + "Leading_phi", //
0227           TString(jet_name) + " leading jet #phi, " + get_eta_range_str()
0228               + ";#phi", 50, -M_PI, M_PI));
0229 
0230   TH1F * lcomp = new TH1F(
0231       //
0232       TString(get_histo_prefix(jet_name)) + "Leading_CompSize", //
0233       TString(jet_name) + " leading jet # of component, " + get_eta_range_str()
0234           + ";Number of component;", 100, 1, 3000);
0235   useLogBins(lcomp->GetXaxis());
0236   hm->registerHisto(lcomp);
0237 
0238   hm->registerHisto(
0239       new TH1F(
0240           //
0241           TString(get_histo_prefix(jet_name)) + "Leading_Mass", //
0242           TString(jet_name) + " leading jet mass, " + get_eta_range_str()
0243               + ";Jet Mass (GeV);", 100, 0, 20));
0244   hm->registerHisto(
0245       new TH1F(
0246           //
0247           TString(get_histo_prefix(jet_name)) + "Leading_CEMC_Ratio", //
0248           TString(jet_name) + " leading jet EMCal ratio, " + get_eta_range_str()
0249               + ";Energy ratio CEMC/Total;", 100, 0, 1.01));
0250   hm->registerHisto(
0251       new TH1F(
0252           //
0253           TString(get_histo_prefix(jet_name)) + "Leading_CEMC_HCalIN_Ratio", //
0254           TString(jet_name) + " leading jet EMCal+HCal_{IN} ratio, "
0255               + get_eta_range_str() + ";Energy ratio (CEMC + HCALIN)/Total;",
0256           100, 0, 1.01));
0257 
0258   // reco jet has no definition for leakages, since leakage is not reconstructed as part of jet energy.
0259   // It is only available for truth jets
0260   hm->registerHisto(
0261       new TH1F(
0262           //
0263           TString(get_histo_prefix(jet_name)) + "Leading_Leakage_Ratio", //
0264           TString(jet_name) + " leading jet leakage ratio, "
0265               + get_eta_range_str() + ";Energy ratio, Back leakage/Total;", 100,
0266           0, 1.01));
0267 
0268   TH1F * h = new TH1F(
0269       //
0270       TString(get_histo_prefix(jet_name)) + "Inclusive_E", //
0271       TString(jet_name) + " inclusive jet E, " + get_eta_range_str()
0272           + ";Total jet energy (GeV)", 100, 1e-3, 100);
0273   useLogBins(h->GetXaxis());
0274   hm->registerHisto(h);
0275   hm->registerHisto(
0276       new TH1F(
0277           //
0278           TString(get_histo_prefix(jet_name)) + "Inclusive_eta", //
0279           TString(jet_name) + " inclusive jet #eta, " + get_eta_range_str()
0280               + ";#eta;Jet energy density", 50, -1, 1));
0281   hm->registerHisto(
0282       new TH1F(
0283           //
0284           TString(get_histo_prefix(jet_name)) + "Inclusive_phi", //
0285           TString(jet_name) + " inclusive jet #phi, " + get_eta_range_str()
0286               + ";#phi;Jet energy density", 50, -M_PI, M_PI));
0287 
0288   TH2F * h2 = NULL;
0289 
0290   h2 = new TH2F(
0291       //
0292       TString(get_histo_prefix(jet_name)) + "Leading_Trk_PtRel", //
0293       TString(jet_name) + " leading jet, " + get_eta_range_str()
0294           + ";j_{T} (GeV/c);Particle Selection", 100, 0, 10, 4, 0.5, 4.5);
0295   i = 1;
0296   h2->GetYaxis()->SetBinLabel(i++, "Electron");
0297   h2->GetYaxis()->SetBinLabel(i++, "Muon");
0298   h2->GetYaxis()->SetBinLabel(i++, "Other Charged");
0299   h2->GetYaxis()->SetBinLabel(i++, "Other Neutral");
0300   hm->registerHisto(h2);
0301 
0302   h2 = new TH2F(
0303       //
0304       TString(get_histo_prefix(jet_name)) + "Leading_Trk_DCA2D", //
0305       TString(jet_name) + " leading jet, " + get_eta_range_str()
0306           + ";DCA_{2D} (cm);Particle Selection", 200, -.2, .2, 4, 0.5, 4.5);
0307   i = 1;
0308   h2->GetYaxis()->SetBinLabel(i++, "Electron");
0309   h2->GetYaxis()->SetBinLabel(i++, "Muon");
0310   h2->GetYaxis()->SetBinLabel(i++, "Other Charged");
0311   h2->GetYaxis()->SetBinLabel(i++, "Other Neutral");
0312   hm->registerHisto(h2);
0313 
0314   h2 = new TH2F(
0315       //
0316       TString(get_histo_prefix(jet_name)) + "Leading_Trk_PoverETotal", //
0317       TString(jet_name) + " leading jet, " + get_eta_range_str()
0318           + ";|P_{part.}|c/E_{Jet};Particle Selection", 120, 0, 1.2, 4, 0.5,
0319       4.5);
0320   i = 1;
0321   h2->GetYaxis()->SetBinLabel(i++, "Electron");
0322   h2->GetYaxis()->SetBinLabel(i++, "Muon");
0323   h2->GetYaxis()->SetBinLabel(i++, "Other Charged");
0324   h2->GetYaxis()->SetBinLabel(i++, "Other Neutral");
0325   hm->registerHisto(h2);
0326 
0327   h2 = new TH2F(
0328       //
0329       TString(get_histo_prefix(jet_name)) + "Leading_Trk_dR", //
0330       TString(jet_name) + " leading jet, " + get_eta_range_str()
0331           + ";#DeltaR;Particle Selection", 100, 0, 1, 4, 0.5, 4.5);
0332   i = 1;
0333   h2->GetYaxis()->SetBinLabel(i++, "Electron");
0334   h2->GetYaxis()->SetBinLabel(i++, "Muon");
0335   h2->GetYaxis()->SetBinLabel(i++, "Other Charged");
0336   h2->GetYaxis()->SetBinLabel(i++, "Other Neutral");
0337   hm->registerHisto(h2);
0338 
0339   return Fun4AllReturnCodes::EVENT_OK;
0340 }
0341 
0342 int
0343 SoftLeptonTaggingTruth::process_Spectrum(PHCompositeNode *topNode,
0344     const std::string & jet_name, const bool is_reco_jet)
0345 {
0346   JetMap* jets = findNode::getClass<JetMap>(topNode, jet_name.c_str());
0347   if (!jets)
0348     {
0349       cout
0350           << "SoftLeptonTaggingTruth::process_Spectrum - Error can not find DST JetMap node "
0351           << jet_name << endl;
0352       exit(-1);
0353     }
0354 
0355   Fun4AllHistoManager *hm = getHistoManager();
0356   assert(hm);
0357 
0358   TH1D* h_norm = dynamic_cast<TH1D*>(hm->getHisto(
0359       get_histo_prefix(jet_name) + "Normalization"));
0360   assert(h_norm);
0361   h_norm->Fill("Event", 1);
0362   h_norm->Fill("Inclusive Jets", jets->size());
0363 
0364   TH1F * ie = dynamic_cast<TH1F*>(hm->getHisto(
0365       (get_histo_prefix(jet_name)) + "Inclusive_E" //
0366           ));
0367   assert(ie);
0368   TH1F * ieta = dynamic_cast<TH1F*>(hm->getHisto(
0369       (get_histo_prefix(jet_name)) + "Inclusive_eta" //
0370           ));
0371   assert(ieta);
0372   TH1F * iphi = dynamic_cast<TH1F*>(hm->getHisto(
0373       (get_histo_prefix(jet_name)) + "Inclusive_phi" //
0374           ));
0375   assert(iphi);
0376 
0377   Jet* leading_jet = NULL;
0378   double max_et = 0;
0379   for (JetMap::Iter iter = jets->begin(); iter != jets->end(); ++iter)
0380     {
0381       Jet* jet = iter->second;
0382       assert(jet);
0383 
0384       if (not jet_acceptance_cut(jet))
0385         continue;
0386 
0387       if (jet->get_et() > max_et)
0388         {
0389           leading_jet = jet;
0390           max_et = jet->get_et();
0391         }
0392 
0393       ie->Fill(jet->get_e());
0394       ieta->Fill(jet->get_eta());
0395       iphi->Fill(jet->get_phi());
0396     }
0397 
0398   if (leading_jet)
0399     {
0400       if (verbosity)
0401         {
0402           cout
0403               << "SoftLeptonTaggingTruth::process_Spectrum - processing leading jet with # comp = "
0404               << leading_jet->size_comp() << endl;
0405           leading_jet->identify();
0406         }
0407 
0408       h_norm->Fill("Leading Jets", 1);
0409 
0410       TH1F * let = dynamic_cast<TH1F*>(hm->getHisto(
0411           (get_histo_prefix(jet_name)) + "Leading_Et" //
0412               ));
0413       assert(let);
0414       TH1F * leta = dynamic_cast<TH1F*>(hm->getHisto(
0415           (get_histo_prefix(jet_name)) + "Leading_eta" //
0416               ));
0417       assert(leta);
0418       TH1F * lphi = dynamic_cast<TH1F*>(hm->getHisto(
0419           (get_histo_prefix(jet_name)) + "Leading_phi" //
0420               ));
0421       assert(lphi);
0422 
0423       TH1F * lcomp = dynamic_cast<TH1F*>(hm->getHisto(
0424           (get_histo_prefix(jet_name)) + "Leading_CompSize" //
0425               ));
0426       assert(lcomp);
0427       TH1F * lmass = dynamic_cast<TH1F*>(hm->getHisto(
0428           (get_histo_prefix(jet_name)) + "Leading_Mass" //
0429               ));
0430       assert(lmass);
0431       TH1F * lcemcr = dynamic_cast<TH1F*>(hm->getHisto(
0432           (get_histo_prefix(jet_name)) + "Leading_CEMC_Ratio" //
0433               ));
0434       assert(lcemcr);
0435       TH1F * lemchcalr = dynamic_cast<TH1F*>(hm->getHisto(
0436           (get_histo_prefix(jet_name)) + "Leading_CEMC_HCalIN_Ratio" //
0437               ));
0438       assert(lemchcalr);
0439       TH1F * lleak = dynamic_cast<TH1F*>(hm->getHisto(
0440           (get_histo_prefix(jet_name)) + "Leading_Leakage_Ratio" //
0441               ));
0442       assert(lleak);
0443 
0444       let->Fill(leading_jet->get_et());
0445       leta->Fill(leading_jet->get_eta());
0446       lphi->Fill(leading_jet->get_phi());
0447       lcomp->Fill(leading_jet->size_comp());
0448       lmass->Fill(leading_jet->get_mass());
0449 
0450       if (is_reco_jet)
0451         { // this is a reco jet
0452           jetevalstacks_map::iterator it_stack = _jetevalstacks.find(jet_name);
0453           assert(it_stack != _jetevalstacks.end());
0454           shared_ptr<JetEvalStack> eval_stack = it_stack->second;
0455           assert(eval_stack);
0456           JetRecoEval* recoeval = eval_stack->get_reco_eval();
0457           assert(recoeval);
0458 
0459           lcemcr->Fill( //
0460               (recoeval->get_energy_contribution(leading_jet, Jet::CEMC_TOWER) //
0461               +//
0462                   recoeval->get_energy_contribution(leading_jet,
0463                       Jet::CEMC_CLUSTER) //
0464               ) / leading_jet->get_e());
0465           lemchcalr->Fill( //
0466               (recoeval->get_energy_contribution(leading_jet, Jet::CEMC_TOWER) //
0467               +//
0468                   recoeval->get_energy_contribution(leading_jet,
0469                       Jet::CEMC_CLUSTER) //
0470                   +//
0471                   recoeval->get_energy_contribution(leading_jet,
0472                       Jet::HCALIN_TOWER) //
0473                   +//
0474                   recoeval->get_energy_contribution(leading_jet,
0475                       Jet::HCALIN_CLUSTER) //
0476               ) / leading_jet->get_e());
0477           // reco jet has no definition for leakages, since leakage is not reconstructed as part of jet energy.
0478           // It is only available for truth jets
0479         }
0480       else
0481         { // this is a truth jet
0482           assert(_jettrutheval);
0483 
0484           double cemc_e = 0;
0485           double hcalin_e = 0;
0486           double bh_e = 0;
0487 
0488           set<PHG4Shower*> showers = _jettrutheval->all_truth_showers(
0489               leading_jet);
0490 
0491           for (set<PHG4Shower*>::const_iterator it = showers.begin();
0492               it != showers.end(); ++it)
0493             {
0494               cemc_e += (*it)->get_edep(
0495                   PHG4HitDefs::get_volume_id("G4HIT_CEMC"));
0496               cemc_e += (*it)->get_edep(
0497                   PHG4HitDefs::get_volume_id("G4HIT_CEMC_ELECTRONICS"));
0498               cemc_e += (*it)->get_edep(
0499                   PHG4HitDefs::get_volume_id("G4HIT_ABSORBER_CEMC"));
0500 
0501               hcalin_e += (*it)->get_edep(
0502                   PHG4HitDefs::get_volume_id("G4HIT_HCALIN"));
0503               hcalin_e += (*it)->get_edep(
0504                   PHG4HitDefs::get_volume_id("G4HIT_ABSORBER_HCALIN"));
0505 
0506               bh_e += (*it)->get_edep(PHG4HitDefs::get_volume_id("G4HIT_BH_1"));
0507               bh_e += (*it)->get_edep(
0508                   PHG4HitDefs::get_volume_id("G4HIT_BH_FORWARD_PLUS"));
0509               bh_e += (*it)->get_edep(
0510                   PHG4HitDefs::get_volume_id("G4HIT_BH_FORWARD_NEG"));
0511             }
0512 
0513           lcemcr->Fill(cemc_e / leading_jet->get_e());
0514           lemchcalr->Fill((cemc_e + hcalin_e) / leading_jet->get_e());
0515           lleak->Fill(bh_e / leading_jet->get_e());
0516 
0517         }
0518 
0519       TH2F * lptrel = dynamic_cast<TH2F*>(hm->getHisto(
0520           (get_histo_prefix(jet_name)) + "Leading_Trk_PtRel" //
0521               ));
0522       assert(lptrel);
0523       TH2F * ldca2d = dynamic_cast<TH2F*>(hm->getHisto(
0524           (get_histo_prefix(jet_name)) + "Leading_Trk_DCA2D" //
0525               ));
0526       assert(ldca2d);
0527       TH2F * lpoe = dynamic_cast<TH2F*>(hm->getHisto(
0528           (get_histo_prefix(jet_name)) + "Leading_Trk_PoverETotal" //
0529               ));
0530       assert(lpoe);
0531       TH2F * ldr = dynamic_cast<TH2F*>(hm->getHisto(
0532           (get_histo_prefix(jet_name)) + "Leading_Trk_dR" //
0533               ));
0534       assert(ldr);
0535 
0536       // jet vector
0537       TVector3 j_v(leading_jet->get_px(), leading_jet->get_py(),
0538           leading_jet->get_pz());
0539 
0540       PHG4TruthInfoContainer::ConstRange primary_range =
0541           _truth_container->GetPrimaryParticleRange();
0542       for (PHG4TruthInfoContainer::ConstIterator particle_iter =
0543           primary_range.first; particle_iter != primary_range.second;
0544           ++particle_iter)
0545         {
0546           PHG4Particle *particle = particle_iter->second;
0547           assert(particle);
0548           const PHG4VtxPoint* primary_vtx = //
0549               _truth_container->GetPrimaryVtx(particle->get_vtx_id());
0550           assert(primary_vtx);
0551 
0552           TVector3 p_v(particle->get_px(), particle->get_py(),
0553               particle->get_pz());
0554 
0555           TVector3 v_v(primary_vtx->get_x(), primary_vtx->get_y(),
0556               primary_vtx->get_z());
0557 
0558           const double dR = sqrt(
0559               (p_v.Eta() - j_v.Eta()) * (p_v.Eta() - j_v.Eta())
0560                   + (p_v.Phi() - j_v.Phi()) * (p_v.Phi() - j_v.Phi()));
0561           const double E_Ratio = p_v.Mag() / leading_jet->get_e();
0562 
0563           const double charge3 = TDatabasePDG::Instance()->GetParticle(
0564               particle->get_pid())->Charge();
0565           TVector3 dca_v(cos(p_v.Phi() + TMath::Pi() / 2),
0566               sin(p_v.Phi() + TMath::Pi() / 2), 0);
0567           const double dca2d = v_v.Dot(dca_v) * copysign(1, charge3);
0568 
0569           if (dR > _jet_match_dR)
0570             continue;
0571           if (abs(dca2d) > _jet_match_dca)
0572             continue;
0573           if (E_Ratio < _jet_match_E_Ratio)
0574             continue;
0575 
0576 
0577           int histo_pid = 0;
0578           if (abs(particle->get_pid())
0579               == TDatabasePDG::Instance()->GetParticle("e-")->PdgCode())
0580             histo_pid = ldr->GetYaxis()->FindBin("Electron");
0581           else if (abs(particle->get_pid())
0582               == TDatabasePDG::Instance()->GetParticle("mu-")->PdgCode())
0583             histo_pid = ldr->GetYaxis()->FindBin("Muon");
0584           else if (charge3 != 0)
0585             histo_pid = ldr->GetYaxis()->FindBin("Other Charged");
0586           else
0587             histo_pid = ldr->GetYaxis()->FindBin("Other Neutral");
0588           assert(histo_pid);
0589 
0590           if (verbosity)
0591             {
0592               cout
0593                   << "SoftLeptonTaggingTruth::process_Spectrum - particle hist ID "
0594                   << histo_pid << ", charge x 3 = " << charge3 << ", dR = "
0595                   << dR << ", E_Ratio = " << E_Ratio << "/"
0596                   << _jet_match_E_Ratio << " for ";
0597               particle->identify();
0598             }
0599 
0600           const double pT_rel = p_v.Dot(j_v) / j_v.Mag();
0601 
0602           ldr->Fill(dR, histo_pid);
0603           lpoe->Fill(E_Ratio, histo_pid);
0604           lptrel->Fill(pT_rel, histo_pid);
0605           ldca2d->Fill(dca2d, histo_pid);
0606         } //       for (PHG4TruthInfoContainer::ConstIterator particle_iter =
0607 
0608     } //   if (leading_jet)
0609 
0610   return Fun4AllReturnCodes::EVENT_OK;
0611 }
0612 
0613 //! Get a pointer to the default hist manager for QA modules
0614 Fun4AllHistoManager *
0615 SoftLeptonTaggingTruth::getHistoManager()
0616 {
0617 
0618   Fun4AllServer *se = Fun4AllServer::instance();
0619   Fun4AllHistoManager *hm = se->getHistoManager("histSoftLeptonTaggingTruth");
0620 
0621   if (not hm)
0622     {
0623 //        cout
0624 //            << "QAHistManagerDef::get_HistoManager - Making Fun4AllHistoManager EMCalAna_HISTOS"
0625 //            << endl;
0626       hm = new Fun4AllHistoManager("histSoftLeptonTaggingTruth");
0627       se->registerHistoManager(hm);
0628     }
0629 
0630   assert(hm);
0631 
0632   return hm;
0633 }
0634 
0635 //! utility function to
0636 void
0637 SoftLeptonTaggingTruth::useLogBins(TAxis * axis)
0638 {
0639   assert(axis);
0640   assert(axis->GetXmin()>0);
0641   assert( axis->GetXmax()>0);
0642 
0643   const int bins = axis->GetNbins();
0644 
0645   Axis_t from = log10(axis->GetXmin());
0646   Axis_t to = log10(axis->GetXmax());
0647   Axis_t width = (to - from) / bins;
0648   vector<Axis_t> new_bins(bins + 1);
0649 
0650   for (int i = 0; i <= bins; i++)
0651     {
0652       new_bins[i] = TMath::Power(10, from + i * width);
0653     }
0654   axis->Set(bins, new_bins.data());
0655 
0656 }
[ 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!