packages/jetbackground/CopyAndSubtractJets.cc

0001 #include "CopyAndSubtractJets.h"
0002
0003 #include "TowerBackground.h"
0004
0005 #include <calobase/RawTower.h>
0006 #include <calobase/RawTowerContainer.h>
0007 #include <calobase/RawTowerDefs.h>
0008 #include <calobase/RawTowerGeom.h>
0009 #include <calobase/RawTowerGeomContainer.h>
0010 #include <calobase/TowerInfo.h>
0011 #include <calobase/TowerInfoContainer.h>
0012
0013 #include <jetbase/Jet.h>
0014 #include <jetbase/JetContainer.h>
0015 #include <jetbase/JetContainerv1.h>
0016
0017 #include <fun4all/Fun4AllReturnCodes.h>
0018 #include <fun4all/SubsysReco.h>
0019
0020 #include <phool/PHCompositeNode.h>
0021 #include <phool/PHIODataNode.h>
0022 #include <phool/PHNode.h>
0023 #include <phool/PHNodeIterator.h>
0024 #include <phool/PHObject.h>
0025 #include <phool/getClass.h>
0026 #include <phool/phool.h>
0027
0028 // standard includes
0029 #include <cmath>
0030 #include <cstdlib>
0031 #include <iostream>
0032 #include <utility>
0033 #include <vector>
0034
0035 CopyAndSubtractJets::CopyAndSubtractJets(const std::string &name)
0036   : SubsysReco(name)
0037 {
0038 }
0039
0040 int CopyAndSubtractJets::InitRun(PHCompositeNode *topNode)
0041 {
0042   CreateNode(topNode);
0043
0044   return Fun4AllReturnCodes::EVENT_OK;
0045 }
0046
0047 int CopyAndSubtractJets::process_event(PHCompositeNode *topNode)
0048 {
0049   if (Verbosity() > 0)
0050   {
0051     std::cout << "CopyAndSubtractJets::process_event: entering, with _use_flow_modulation = " << _use_flow_modulation << std::endl;
0052   }
0053
0054   // pull out needed calo tower info
0055   RawTowerContainer *towersEM3 = nullptr;
0056   RawTowerContainer *towersIH3 = nullptr;
0057   RawTowerContainer *towersOH3 = nullptr;
0058   TowerInfoContainer *towerinfosEM3 = nullptr;
0059   TowerInfoContainer *towerinfosIH3 = nullptr;
0060   TowerInfoContainer *towerinfosOH3 = nullptr;
0061   if (m_use_towerinfo)
0062   {
0063     EMTowerName = m_towerNodePrefix + "_CEMC_RETOWER";
0064     IHTowerName = m_towerNodePrefix + "_HCALIN";
0065     OHTowerName = m_towerNodePrefix + "_HCALOUT";
0066     towerinfosEM3 = findNode::getClass<TowerInfoContainer>(topNode, EMTowerName);
0067     towerinfosIH3 = findNode::getClass<TowerInfoContainer>(topNode, IHTowerName);
0068     towerinfosOH3 = findNode::getClass<TowerInfoContainer>(topNode, OHTowerName);
0069     if (!towerinfosEM3)
0070     {
0071       std::cout << "CopyAndSubtractJets::process_event: Cannot find node " << EMTowerName << std::endl;
0072       exit(1);
0073     }
0074     if (!towerinfosIH3)
0075     {
0076       std::cout << "CopyAndSubtractJets::process_event: Cannot find node " << IHTowerName << std::endl;
0077       exit(1);
0078     }
0079     if (!towerinfosOH3)
0080     {
0081       std::cout << "CopyAndSubtractJets::process_event: Cannot find node " << OHTowerName << std::endl;
0082       exit(1);
0083     }
0084   }
0085   else
0086   {
0087     towersEM3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC_RETOWER");
0088     towersIH3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN");
0089     towersOH3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT");
0090   }
0091
0092   RawTowerGeomContainer *geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0093   RawTowerGeomContainer *geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0094
0095   // pull out jets and background
0096   JetContainer *unsub_jets;
0097   JetContainer *sub_jets;
0098   TowerBackground *background;
0099   if (m_use_towerinfo)
0100   {
0101     unsub_jets = findNode::getClass<JetContainer>(topNode, "AntiKt_TowerInfo_HIRecoSeedsRaw_r02");
0102     sub_jets = findNode::getClass<JetContainer>(topNode, "AntiKt_TowerInfo_HIRecoSeedsSub_r02");
0103     background = findNode::getClass<TowerBackground>(topNode, "TowerInfoBackground_Sub1");
0104   }
0105   else
0106   {
0107     unsub_jets = findNode::getClass<JetContainer>(topNode, "AntiKt_Tower_HIRecoSeedsRaw_r02");
0108     sub_jets = findNode::getClass<JetContainer>(topNode, "AntiKt_Tower_HIRecoSeedsSub_r02");
0109     background = findNode::getClass<TowerBackground>(topNode, "TowerBackground_Sub1");
0110   }
0111   std::vector<float> background_UE_0 = background->get_UE(0);
0112   std::vector<float> background_UE_1 = background->get_UE(1);
0113   std::vector<float> background_UE_2 = background->get_UE(2);
0114
0115   float background_v2 = background->get_v2();
0116   float background_Psi2 = background->get_Psi2();
0117
0118   if (Verbosity() > 0)
0119   {
0120     std::cout << "CopyAndSubtractJets::process_event: entering with # unsubtracted jets = " << unsub_jets->size() << std::endl;
0121     std::cout << "CopyAndSubtractJets::process_event: entering with # subtracted jets = " << sub_jets->size() << std::endl;
0122   }
0123
0124   // iterate over old jets
0125   int ijet = 0;
0126   for (auto *this_jet : *unsub_jets)
0127   {
0128     float this_pt = this_jet->get_pt();
0129     float this_phi = this_jet->get_phi();
0130     float this_eta = this_jet->get_eta();
0131
0132     float new_total_px = 0;
0133     float new_total_py = 0;
0134     float new_total_pz = 0;
0135     float new_total_e = 0;
0136
0137     // if (this_jet->get_pt() < 5) continue;
0138
0139     if (Verbosity() > 1 && this_jet->get_pt() > 5)
0140     {
0141       std::cout << "CopyAndSubtractJets::process_event: unsubtracted jet with pt / eta / phi = " << this_pt << " / " << this_eta << " / " << this_phi << std::endl;
0142     }
0143
0144     for (const auto &comp : this_jet->get_comp_vec())
0145     {
0146       RawTower *tower = nullptr;
0147       RawTowerGeom *tower_geom = nullptr;
0148       TowerInfo *towerinfo = nullptr;
0149
0150       double comp_e = 0;
0151       double comp_eta = 0;
0152       double comp_phi = 0;
0153
0154       int comp_ieta = 0;
0155
0156       double comp_background = 0;
0157
0158       if (m_use_towerinfo)
0159       {
0160         if (comp.first == 5 || comp.first == 26)
0161         {
0162           towerinfo = towerinfosIH3->get_tower_at_channel(comp.second);
0163           unsigned int towerkey = towerinfosIH3->encode_key(comp.second);
0164           comp_ieta = towerinfosIH3->getTowerEtaBin(towerkey);
0165           int comp_iphi = towerinfosIH3->getTowerPhiBin(towerkey);
0166           const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, comp_ieta, comp_iphi);
0167
0168           tower_geom = geomIH->get_tower_geometry(key);
0169           comp_background = background_UE_1.at(comp_ieta);
0170         }
0171         else if (comp.first == 7 || comp.first == 27)
0172         {
0173           towerinfo = towerinfosOH3->get_tower_at_channel(comp.second);
0174           unsigned int towerkey = towerinfosOH3->encode_key(comp.second);
0175           comp_ieta = towerinfosOH3->getTowerEtaBin(towerkey);
0176           int comp_iphi = towerinfosOH3->getTowerPhiBin(towerkey);
0177           const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, comp_ieta, comp_iphi);
0178           tower_geom = geomOH->get_tower_geometry(key);
0179           comp_background = background_UE_2.at(comp_ieta);
0180         }
0181         else if (comp.first == 13 || comp.first == 28)
0182         {
0183           towerinfo = towerinfosEM3->get_tower_at_channel(comp.second);
0184           unsigned int towerkey = towerinfosEM3->encode_key(comp.second);
0185           comp_ieta = towerinfosEM3->getTowerEtaBin(towerkey);
0186           int comp_iphi = towerinfosEM3->getTowerPhiBin(towerkey);
0187           const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, comp_ieta, comp_iphi);
0188
0189           tower_geom = geomIH->get_tower_geometry(key);
0190           comp_background = background_UE_0.at(comp_ieta);
0191         }
0192         if (towerinfo)
0193         {
0194           comp_e = towerinfo->get_energy();
0195         }
0196       }
0197       else
0198       {
0199         if (comp.first == 5)
0200         {
0201           tower = towersIH3->getTower(comp.second);
0202           tower_geom = geomIH->get_tower_geometry(tower->get_key());
0203
0204           comp_ieta = geomIH->get_etabin(tower_geom->get_eta());
0205           comp_background = background_UE_1.at(comp_ieta);
0206         }
0207         else if (comp.first == 7)
0208         {
0209           tower = towersOH3->getTower(comp.second);
0210           tower_geom = geomOH->get_tower_geometry(tower->get_key());
0211
0212           comp_ieta = geomOH->get_etabin(tower_geom->get_eta());
0213           comp_background = background_UE_2.at(comp_ieta);
0214         }
0215         else if (comp.first == 13)
0216         {
0217           tower = towersEM3->getTower(comp.second);
0218           tower_geom = geomIH->get_tower_geometry(tower->get_key());
0219
0220           comp_ieta = geomIH->get_etabin(tower_geom->get_eta());
0221           comp_background = background_UE_0.at(comp_ieta);
0222         }
0223         if (tower)
0224         {
0225           comp_e = tower->get_energy();
0226         }
0227       }
0228
0229       if (tower_geom)
0230       {
0231         comp_eta = tower_geom->get_eta();
0232         comp_phi = tower_geom->get_phi();
0233       }
0234
0235       if (Verbosity() > 4 && this_jet->get_pt() > 5)
0236       {
0237         std::cout << "CopyAndSubtractJets::process_event: --> constituent in layer " << comp.first << ", has unsub E = " << comp_e << ", is at ieta #" << comp_ieta << ", and has UE = " << comp_background << std::endl;
0238       }
0239
0240       // flow modulate background if turned on
0241       if (_use_flow_modulation)
0242       {
0243         comp_background = comp_background * (1 + 2 * background_v2 * cos(2 * (comp_phi - background_Psi2)));
0244         if (Verbosity() > 4 && this_jet->get_pt() > 5)
0245         {
0246           std::cout << "CopyAndSubtractJets::process_event: --> --> flow mod, at phi = " << comp_phi << ", v2 and Psi2 are = " << background_v2 << " , " << background_Psi2 << ", UE after modulation = " << comp_background << std::endl;
0247         }
0248       }
0249
0250       // update constituent energy based on the background
0251       double comp_sub_e = comp_e - comp_background;
0252
0253       // now define new kinematics
0254
0255       double comp_px = comp_sub_e / cosh(comp_eta) * cos(comp_phi);
0256       double comp_py = comp_sub_e / cosh(comp_eta) * sin(comp_phi);
0257       double comp_pz = comp_sub_e * tanh(comp_eta);
0258
0259       new_total_px += comp_px;
0260       new_total_py += comp_py;
0261       new_total_pz += comp_pz;
0262       new_total_e += comp_sub_e;
0263     }
0264
0265     auto *new_jet = sub_jets->add_jet();  // returns a new Jet_v2
0266
0267     new_jet->set_px(new_total_px);
0268     new_jet->set_py(new_total_py);
0269     new_jet->set_pz(new_total_pz);
0270     new_jet->set_e(new_total_e);
0271     new_jet->set_id(ijet);
0272
0273     if (Verbosity() > 1 && this_pt > 5)
0274     {
0275       std::cout << "CopyAndSubtractJets::process_event: old jet #" << ijet << ", old px / py / pz / e = " << this_jet->get_px() << " / " << this_jet->get_py() << " / " << this_jet->get_pz() << " / " << this_jet->get_e() << std::endl;
0276       std::cout << "CopyAndSubtractJets::process_event: new jet #" << ijet << ", new px / py / pz / e = " << new_jet->get_px() << " / " << new_jet->get_py() << " / " << new_jet->get_pz() << " / " << new_jet->get_e() << std::endl;
0277     }
0278
0279     ijet++;
0280   }
0281
0282   if (Verbosity() > 0)
0283   {
0284     std::cout << "CopyAndSubtractJets::process_event: exiting with # subtracted jets = " << sub_jets->size() << std::endl;
0285   }
0286
0287   return Fun4AllReturnCodes::EVENT_OK;
0288 }
0289
0290 int CopyAndSubtractJets::End(PHCompositeNode * /*topNode*/)
0291 {
0292   return Fun4AllReturnCodes::EVENT_OK;
0293 }
0294
0295 int CopyAndSubtractJets::CreateNode(PHCompositeNode *topNode)
0296 {
0297   PHNodeIterator iter(topNode);
0298
0299   // Looking for the DST node
0300   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0301   if (!dstNode)
0302   {
0303     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0304     return Fun4AllReturnCodes::ABORTRUN;
0305   }
0306
0307   // Looking for the ANTIKT node
0308   PHCompositeNode *antiktNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "ANTIKT"));
0309   if (!antiktNode)
0310   {
0311     std::cout << PHWHERE << "ANTIKT node not found, doing nothing." << std::endl;
0312   }
0313
0314   // Looking for the TOWER node
0315   PHCompositeNode *towerNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "TOWER"));
0316   if (!towerNode)
0317   {
0318     std::cout << PHWHERE << "TOWER node not found, doing nothing." << std::endl;
0319   }
0320
0321   // store the new jet collection
0322   JetContainer *test_jets;
0323   if (m_use_towerinfo)
0324   {
0325     test_jets = findNode::getClass<JetContainer>(topNode, "AntiKt_TowerInfo_HIRecoSeedsSub_r02");
0326   }
0327   else
0328   {
0329     test_jets = findNode::getClass<JetContainer>(topNode, "AntiKt_Tower_HIRecoSeedsSub_r02");
0330   }
0331   if (!test_jets)
0332   {
0333     JetContainer *sub_jets = new JetContainerv1();
0334     PHIODataNode<PHObject> *subjetNode;
0335     if (m_use_towerinfo)
0336     {
0337       if (Verbosity() > 0)
0338       {
0339         std::cout << "CopyAndSubtractJets::CreateNode : creating AntiKt_TowerInfo_HIRecoSeedsSub_r02 node " << std::endl;
0340       }
0341       subjetNode = new PHIODataNode<PHObject>(sub_jets, "AntiKt_TowerInfo_HIRecoSeedsSub_r02", "PHObject");
0342     }
0343     else
0344     {
0345       if (Verbosity() > 0)
0346       {
0347         std::cout << "CopyAndSubtractJets::CreateNode : creating AntiKt_Tower_HIRecoSeedsSub_r02 node " << std::endl;
0348       }
0349       subjetNode = new PHIODataNode<PHObject>(sub_jets, "AntiKt_Tower_HIRecoSeedsSub_r02", "PHObject");
0350     }
0351     towerNode->addNode(subjetNode);
0352   }
0353   else
0354   {
0355     std::cout << "CopyAndSubtractJets::CreateNode : AntiKt_Tower_HIRecoSeedsSub_r02 already exists! " << std::endl;
0356   }
0357
0358   return Fun4AllReturnCodes::EVENT_OK;
0359 }