analysis/TreeMaker/TreeMakerGetJets.C

0001 #include <TreeMaker.h>
0002
0003 #include <phool/getClass.h>
0004 #include <phool/PHCompositeNode.h>
0005
0006 #include <TLorentzVector.h>
0007 #include <g4jets/JetMap.h>
0008 #include <g4jets/Jet.h>
0009 #include <jetbackground/TowerBackground.h>
0010
0011
0012
0013 int TreeMaker::GetTruthJets(PHCompositeNode *topNode)
0014 {
0015
0016   // ------------------------------------------------------------------------------
0017   // --- get truth jet information
0018   // ------------------------------------------------------------------------------
0019
0020   JetMap* truth2_jets = findNode::getClass<JetMap>(topNode,"AntiKt_Truth_r02");
0021   //cout << "R = 0.2 truth jets has size " << truth2_jets->size() << endl;
0022   for (JetMap::Iter iter = truth2_jets->begin(); iter != truth2_jets->end(); ++iter)
0023     {
0024       Jet* this_jet = iter->second;
0025       float this_e = this_jet->get_e();
0026       float this_pt = this_jet->get_pt();
0027       float this_phi = this_jet->get_phi();
0028       float this_eta = this_jet->get_eta();
0029       if (this_jet->get_pt() < 5 || fabs(this_eta) > 5) continue;
0030       b_truthjet2_e[ b_truthjet2_n ] = this_e;
0031       b_truthjet2_pt[ b_truthjet2_n ] = this_pt;
0032       b_truthjet2_eta[ b_truthjet2_n ] = this_eta;
0033       b_truthjet2_phi[ b_truthjet2_n ] = this_phi;
0034       //cout << " truth R=0.2 jet # " << b_truthjet2_n << " pt / eta / phi / m = " << this_pt << " / " << this_eta << " / " << this_phi << endl;
0035       b_truthjet2_n++;
0036     } // loop over R = 0.2 truth jets
0037
0038   JetMap* truth3_jets = findNode::getClass<JetMap>(topNode,"AntiKt_Truth_r03");
0039   //cout << "R = 0.3 truth jets has size " << truth3_jets->size() << endl;
0040   for (JetMap::Iter iter = truth3_jets->begin(); iter != truth3_jets->end(); ++iter)
0041     {
0042       Jet* this_jet = iter->second;
0043       float this_e = this_jet->get_e();
0044       float this_pt = this_jet->get_pt();
0045       float this_phi = this_jet->get_phi();
0046       float this_eta = this_jet->get_eta();
0047       if (this_jet->get_pt() < 5 || fabs(this_eta) > 5) continue;
0048       b_truthjet3_e[ b_truthjet3_n ] = this_e;
0049       b_truthjet3_pt[ b_truthjet3_n ] = this_pt;
0050       b_truthjet3_eta[ b_truthjet3_n ] = this_eta;
0051       b_truthjet3_phi[ b_truthjet3_n ] = this_phi;
0052       //cout << " truth R=0.3 jet # " << b_truthjet3_n << " pt / eta / phi / m = " << this_pt << " / " << this_eta << " / " << this_phi << endl;
0053       b_truthjet3_n++;
0054     } // loop over R = 0.3 truth jets
0055
0056   JetMap* truth4_jets = findNode::getClass<JetMap>(topNode,"AntiKt_Truth_r04");
0057   //cout << "R = 0.4 truth jets has size " << truth4_jets->size() << endl;
0058   for (JetMap::Iter iter = truth4_jets->begin(); iter != truth4_jets->end(); ++iter)
0059     {
0060       Jet* this_jet = iter->second;
0061       float this_e = this_jet->get_e();
0062       float this_pt = this_jet->get_pt();
0063       float this_phi = this_jet->get_phi();
0064       float this_eta = this_jet->get_eta();
0065       if (this_jet->get_pt() < 5 || fabs(this_eta) > 5) continue;
0066       b_truthjet4_e[ b_truthjet4_n ] = this_e;
0067       b_truthjet4_pt[ b_truthjet4_n ] = this_pt;
0068       b_truthjet4_eta[ b_truthjet4_n ] = this_eta;
0069       b_truthjet4_phi[ b_truthjet4_n ] = this_phi;
0070       //cout << " truth R=0.4 jet # " << b_truthjet4_n << " pt / eta / phi / m = " << this_pt << " / " << this_eta << " / " << this_phi << endl;
0071       b_truthjet4_n++;
0072     } // loop over R = 0.4 truth jets
0073
0074   JetMap* truth5_jets = findNode::getClass<JetMap>(topNode,"AntiKt_Truth_r05");
0075   //cout << "R = 0.5 truth jets has size " << truth5_jets->size() << endl;
0076   for (JetMap::Iter iter = truth5_jets->begin(); iter != truth5_jets->end(); ++iter)
0077     {
0078       Jet* this_jet = iter->second;
0079       float this_e = this_jet->get_e();
0080       float this_pt = this_jet->get_pt();
0081       float this_phi = this_jet->get_phi();
0082       float this_eta = this_jet->get_eta();
0083       if (this_jet->get_pt() < 5 || fabs(this_eta) > 5) continue;
0084       b_truthjet5_e[ b_truthjet5_n ] = this_e;
0085       b_truthjet5_pt[ b_truthjet5_n ] = this_pt;
0086       b_truthjet5_eta[ b_truthjet5_n ] = this_eta;
0087       b_truthjet5_phi[ b_truthjet5_n ] = this_phi;
0088       //cout << " truth R=0.5 jet # " << b_truthjet5_n << " pt / eta / phi / m = " << this_pt << " / " << this_eta << " / " << this_phi << endl;
0089       b_truthjet5_n++;
0090     } // loop over R = 0.5 truth jets
0091
0092   return 0;
0093
0094 }
0095
0096
0097
0098 int TreeMaker::GetRecoJets(PHCompositeNode *topNode)
0099 {
0100
0101   // ------------------------------------------------------------------------------
0102   // --- get reconstructed jet information
0103   // ------------------------------------------------------------------------------
0104
0105   JetMap* nsreco2_jets = findNode::getClass<JetMap>(topNode,"AntiKt_Tower_r02");
0106   //cout << "reco jets R=0.2 has size " << nsreco2_jets->size() << endl;
0107   for (JetMap::Iter iter = nsreco2_jets->begin(); iter != nsreco2_jets->end(); ++iter)
0108     {
0109       Jet* this_jet = iter->second;
0110       float this_e = this_jet->get_e();
0111       float this_pt = this_jet->get_pt();
0112       float this_phi = this_jet->get_phi();
0113       float this_eta = this_jet->get_eta();
0114       if (this_jet->get_pt() < 10 || fabs(this_eta) > 5) continue; // stricter pt cut for unsubtracted
0115       b_jet2_e[ b_jet2_n ] = this_e;
0116       b_jet2_pt[ b_jet2_n ] = this_pt;
0117       b_jet2_eta[ b_jet2_n ] = this_eta;
0118       b_jet2_phi[ b_jet2_n ] = this_phi;
0119       //cout << " pp reco R=0.2 jet #" << b_jet2_n << ", pt / eta / phi = " << this_pt << " / " << this_eta << " / " << this_phi << endl;
0120       b_jet2_n++;
0121     } // loop over R=0.2 jets
0122
0123   JetMap* nsreco3_jets = findNode::getClass<JetMap>(topNode,"AntiKt_Tower_r03");
0124   //cout << "reco jets R=0.3 has size " << nsreco3_jets->size() << endl;
0125   for (JetMap::Iter iter = nsreco3_jets->begin(); iter != nsreco3_jets->end(); ++iter)
0126     {
0127       Jet* this_jet = iter->second;
0128       float this_e = this_jet->get_e();
0129       float this_pt = this_jet->get_pt();
0130       float this_phi = this_jet->get_phi();
0131       float this_eta = this_jet->get_eta();
0132       if (this_jet->get_pt() < 10 || fabs(this_eta) > 5) continue; // stricter pt cut for unsubtracted
0133       b_jet3_e[ b_jet3_n ] = this_e;
0134       b_jet3_pt[ b_jet3_n ] = this_pt;
0135       b_jet3_eta[ b_jet3_n ] = this_eta;
0136       b_jet3_phi[ b_jet3_n ] = this_phi;
0137       //cout << " pp reco R=0.3 jet #" << b_jet3_n << ", pt / eta / phi = " << this_pt << " / " << this_eta << " / " << this_phi << endl;
0138       b_jet3_n++;
0139     } // loop over R=0.3 jets
0140
0141   JetMap* nsreco4_jets = findNode::getClass<JetMap>(topNode,"AntiKt_Tower_r04");
0142   //cout << "reco jets R=0.4 has size " << nsreco4_jets->size() << endl;
0143   for (JetMap::Iter iter = nsreco4_jets->begin(); iter != nsreco4_jets->end(); ++iter)
0144     {
0145       Jet* this_jet = iter->second;
0146       float this_e = this_jet->get_e();
0147       float this_pt = this_jet->get_pt();
0148       float this_phi = this_jet->get_phi();
0149       float this_eta = this_jet->get_eta();
0150       if (this_jet->get_pt() < 10 || fabs(this_eta) > 5) continue; // stricter pt cut for unsubtracted
0151       b_jet4_e[ b_jet4_n ] = this_e;
0152       b_jet4_pt[ b_jet4_n ] = this_pt;
0153       b_jet4_eta[ b_jet4_n ] = this_eta;
0154       b_jet4_phi[ b_jet4_n ] = this_phi;
0155       //cout << " pp reco R=0.4 jet #" << b_jet4_n << ", pt / eta / phi = " << this_pt << " / " << this_eta << " / " << this_phi << endl;
0156       b_jet4_n++;
0157     } // loop over R=0.4 jets
0158
0159   JetMap* nsreco5_jets = findNode::getClass<JetMap>(topNode,"AntiKt_Tower_r05");
0160   //cout << "reco jets R=0.5 has size " << nsreco5_jets->size() << endl;
0161   for (JetMap::Iter iter = nsreco5_jets->begin(); iter != nsreco5_jets->end(); ++iter)
0162     {
0163       Jet* this_jet = iter->second;
0164       float this_e = this_jet->get_e();
0165       float this_pt = this_jet->get_pt();
0166       float this_phi = this_jet->get_phi();
0167       float this_eta = this_jet->get_eta();
0168       if (this_jet->get_pt() < 10 || fabs(this_eta) > 5) continue; // stricter pt cut for unsubtracted
0169       b_jet5_e[ b_jet5_n ] = this_e;
0170       b_jet5_pt[ b_jet5_n ] = this_pt;
0171       b_jet5_eta[ b_jet5_n ] = this_eta;
0172       b_jet5_phi[ b_jet5_n ] = this_phi;
0173       //cout << " pp reco R=0.2 jet #" << b_jet5_n << ", pt / eta / phi = " << this_pt << " / " << this_eta << " / " << this_phi << endl;
0174       b_jet5_n++;
0175     } // loop over R=0.5 jets
0176
0177
0178
0179   // ------------------------------------------------------------------------------
0180   // --- get reconstructed jet information (modified)
0181   // ------------------------------------------------------------------------------
0182
0183   JetMap* modnsreco2_jets = findNode::getClass<JetMap>(topNode,"AntiKt_Tower_Mod_r02");
0184   //cout << "reco jets R=0.2 has size " << modnsreco2_jets->size() << endl;
0185   for (JetMap::Iter iter = modnsreco2_jets->begin(); iter != modnsreco2_jets->end(); ++iter)
0186     {
0187       Jet* this_jet = iter->second;
0188       float this_e = this_jet->get_e();
0189       float this_pt = this_jet->get_pt();
0190       float this_phi = this_jet->get_phi();
0191       float this_eta = this_jet->get_eta();
0192       if (this_jet->get_pt() < 10 || fabs(this_eta) > 5) continue; // stricter pt cut for unsubtracted
0193       b_modjet2_e[ b_modjet2_n ] = this_e;
0194       b_modjet2_pt[ b_modjet2_n ] = this_pt;
0195       b_modjet2_eta[ b_modjet2_n ] = this_eta;
0196       b_modjet2_phi[ b_modjet2_n ] = this_phi;
0197       //cout << " pp reco R=0.2 jet #" << b_jet2_n << ", pt / eta / phi = " << this_pt << " / " << this_eta << " / " << this_phi << endl;
0198       b_modjet2_n++;
0199     } // loop over R=0.2 jets
0200
0201
0202
0203   return 0;
0204
0205 }