packages/jetbase/JetCalib.cc

0001 #include "JetCalib.h"
0002
0003 #include "JetContainerv1.h"
0004
0005 #include <cdbobjects/CDBTF.h>  // for CDBTF1
0006
0007 #include <ffamodules/CDBInterface.h>
0008 #include <ffaobjects/EventHeader.h>
0009
0010 #include <globalvertex/GlobalVertex.h>
0011 #include <globalvertex/GlobalVertexMap.h>
0012
0013 #include <fun4all/Fun4AllReturnCodes.h>
0014 #include <fun4all/SubsysReco.h>  // for SubsysReco
0015
0016 #include <phool/PHCompositeNode.h>
0017 #include <phool/PHIODataNode.h>    // for PHIODataNode
0018 #include <phool/PHNode.h>          // for PHNode
0019 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0020 #include <phool/PHObject.h>        // for PHObject
0021 #include <phool/getClass.h>
0022 #include <phool/phool.h>
0023 #include <phool/recoConsts.h>
0024
0025 #include <TF1.h>
0026 #include <TSystem.h>
0027
0028 #include <boost/format.hpp>
0029
0030 #include <cmath>
0031 #include <cstdlib>    // for exit
0032 #include <exception>  // for exception
0033 #include <iostream>   // for operator<<, basic_ostream
0034 #include <stdexcept>  // for runtime_error
0035
0036 JetCalib::JetCalib(const std::string &name)
0037   : SubsysReco(name)
0038 {
0039   if (Verbosity() > 0)
0040   {
0041     std::cout << "JetCalib::JetCalib(const std::string &name) Calling ctor." << std::endl;
0042   }
0043 }
0044
0045 JetCalib::~JetCalib()
0046 {
0047   if (Verbosity() > 0)
0048   {
0049     std::cout << "JetCalib::~JetCalib() : Calling dtor." << std::endl;
0050   }
0051 }
0052
0053 int JetCalib::InitRun(PHCompositeNode *topNode)
0054 {
0055   // Create calib jet node.
0056   CreateNodeTree(topNode);
0057
0058   // Load calibration file and function.
0059   if (fetchCalibDir("Default").empty())
0060   {
0061     std::cout << "JetCalib::InitRun(PHCompositeNode *topNode) : No default calibration available! Will apply calib factor 1." << std::endl;
0062   }
0063   else
0064   {
0065     m_JetCalibFile = new CDBTF(fetchCalibDir("Default"));
0066     if (m_JetCalibFile)
0067     {
0068       m_JetCalibFile->LoadCalibrations();
0069       std::string JetCalibFunc_name = "JES_Calib_Func_R0" + std::to_string((int)(jet_radius * 10));
0070       if (!ApplyZvrtxDependentCalib && !ApplyEtaDependentCalib)
0071       {
0072         int nZvrtxBins = 1;
0073         int nEtaBins = 1;
0074         m_JetCalibFunc.resize(nZvrtxBins);
0075         for (auto &row : m_JetCalibFunc)
0076         {
0077           row.resize(nEtaBins, nullptr);
0078         }
0079         TF1 *func_temp = m_JetCalibFile->getTF(JetCalibFunc_name + "_Default");
0080         if (!func_temp)
0081         {
0082           std::cout << "JetCalib::InitRun(PHCompositeNode *topNode) : Could not find calibration function: " << JetCalibFunc_name << "_Default" << std::endl;
0083           exit(1);
0084         }
0085         m_JetCalibFunc[0][0] = func_temp;
0086       }
0087       else if (ApplyZvrtxDependentCalib && !ApplyEtaDependentCalib)
0088       {
0089         int nZvrtxBins = 5;
0090         int nEtaBins = 1;
0091         m_JetCalibFunc.resize(nZvrtxBins);
0092         for (auto &row : m_JetCalibFunc)
0093         {
0094           row.resize(nEtaBins, nullptr);
0095         }
0096         for (int iz = 0; iz < nZvrtxBins; ++iz)
0097         {
0098           TF1 *func_temp = m_JetCalibFile->getTF(JetCalibFunc_name + "_Z" + std::to_string(iz));
0099           if (!func_temp)
0100           {
0101             std::cout << "JetCalib::InitRun(PHCompositeNode *topNode) : Could not find calibration function: " << JetCalibFunc_name << "_Z" << iz << std::endl;
0102             exit(1);
0103           }
0104           m_JetCalibFunc[iz][0] = func_temp;
0105         }
0106       }
0107       else if (ApplyEtaDependentCalib)
0108       {
0109         if (!ApplyZvrtxDependentCalib)
0110         {
0111           std::cout << "JetCalib::InitRun(PHCompositeNode *topNode) : Must apply Zvrtx dependent calibration to apply eta dependent calibration. Applying Zvrtx + eta dependent calibration." << std::endl;
0112         }
0113         int nZvrtxBins = 5;
0114         int nEtaBins = 4;
0115         m_JetCalibFunc.resize(nZvrtxBins);
0116         for (auto &row : m_JetCalibFunc)
0117         {
0118           row.resize(nEtaBins, nullptr);
0119         }
0120         for (int iz = 0; iz < nZvrtxBins; ++iz)
0121         {
0122           if (iz < 3) {
0123             for (int ieta = 0; ieta < nEtaBins; ++ieta)
0124             {
0125               TF1 *func_temp = m_JetCalibFile->getTF(JetCalibFunc_name + "_Z" + std::to_string(iz) + "_Eta" + std::to_string(ieta));
0126               if (!func_temp)
0127               {
0128                 std::cout << "JetCalib::InitRun(PHCompositeNode *topNode) : Could not find calibration function: " << JetCalibFunc_name << "_Z" << iz << "_Eta" << ieta << std::endl;
0129                 exit(1);
0130               }
0131               m_JetCalibFunc[iz][ieta] = func_temp;
0132             }
0133           }
0134           else
0135           {
0136             TF1 *func_temp = m_JetCalibFile->getTF(JetCalibFunc_name + "_Z" + std::to_string(iz));
0137             if (!func_temp)
0138             {
0139               std::cout << "JetCalib::InitRun(PHCompositeNode *topNode) : Could not find calibration function: " << JetCalibFunc_name << "_Z" << iz << std::endl;
0140               exit(1);
0141             }
0142             m_JetCalibFunc[iz][0] = func_temp;
0143           }
0144         }
0145       }
0146     }
0147     else
0148     {
0149       std::cout << "JetCalib::InitRun(PHCompositeNode *topNode) : Could not open calibration file!" << std::endl;
0150     }
0151   }
0152
0153   if (Verbosity() > 0)
0154   {
0155     std::cout << "JetCalib::InitRun(PHCompositeNode *topNode) - InitRun with inputNode: " << m_inputNode << " outputNode: " << m_outputNode << std::endl;
0156   }
0157   return Fun4AllReturnCodes::EVENT_OK;
0158 }
0159
0160 int JetCalib::process_event(PHCompositeNode *topNode)
0161 {
0162   // Get raw and calib jet nodes.
0163   JetContainer *_raw_jets = findNode::getClass<JetContainer>(topNode, m_inputNode);
0164   if (!_raw_jets)
0165   {
0166     std::cout << "JetCalib::process_event(PHCompositeNode *topNode) : Raw jet node not found! Cannot apply calibration." << std::endl;
0167     return Fun4AllReturnCodes::ABORTRUN;
0168   }
0169   JetContainer *_calib_jets = findNode::getClass<JetContainer>(topNode, m_outputNode);
0170   if (!_calib_jets)
0171   {
0172     std::cout << "JetCalib::process_event(PHCompositeNode *topNode) : Calib jet node not found! Cannot apply calibration." << std::endl;
0173     return Fun4AllReturnCodes::ABORTRUN;
0174   }
0175   GlobalVertexMap *vertexmap = nullptr;
0176   if (ApplyZvrtxDependentCalib)
0177   {
0178     vertexmap = findNode::getClass<GlobalVertexMap>(topNode, m_zvrtxNode);
0179     if (!vertexmap)
0180     {
0181       std::cout << "JetCalib::process_event(PHCompositeNode *topNode) : GlobalVertexMap node not found! Cannot apply Z-vertex dependent calibration." << std::endl;
0182       return Fun4AllReturnCodes::ABORTRUN;
0183     }
0184   }
0185   // Apply calibration to jets and add to calib jet node.
0186   int ijet = 0;
0187   for (auto *jet : *_raw_jets)
0188   {
0189     float pt = jet->get_pt();
0190     float eta = jet->get_eta();
0191     float phi = jet->get_phi();
0192     float zvertex = -999.0;
0193     if (ApplyZvrtxDependentCalib)
0194     {
0195       if (!vertexmap->empty())
0196       {
0197         GlobalVertex *vtx = vertexmap->begin()->second;
0198         if (vtx)
0199         {
0200           zvertex = vtx->get_z();
0201         }
0202       }
0203       else
0204       {
0205         if (Verbosity() > 0)
0206         {
0207           std::cout << "JetCalib::process_event(PHCompositeNode *topNode) : GlobalVertexMap is empty. Assign zvertex = 0." << std::endl;
0208         }
0209       }
0210     }
0211     float calib_pt = doCalibration(m_JetCalibFunc, pt, zvertex, eta);
0212     auto *calib_jet = _calib_jets->add_jet();
0213     calib_jet->set_px(calib_pt * std::cos(phi));
0214     calib_jet->set_py(calib_pt * std::sin(phi));
0215     calib_jet->set_pz(calib_pt * std::sinh(eta));
0216     calib_jet->set_id(ijet);
0217     calib_jet->set_isCalib(1);
0218     ijet++;
0219   }
0220
0221   if (Verbosity() > 0)
0222   {
0223     std::cout << "JetCalib::process_event(PHCompositeNode *topNode) : nRawJets: " << _raw_jets->size() << " nCalibJets: " << _calib_jets->size() << std::endl;
0224     if (_calib_jets->size() != _raw_jets->size())
0225     {
0226       std::cout << "JetCalib::process_event(PHCompositeNode *topNode) : different number of raw jets vs. calib jets! Something is amiss! " << std::endl;
0227     }
0228   }
0229
0230   return Fun4AllReturnCodes::EVENT_OK;
0231 }
0232
0233 int JetCalib::CreateNodeTree(PHCompositeNode *topNode)
0234 {
0235   // Check nodes.
0236   PHNodeIterator iter(topNode);
0237   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0238   if (!dstNode)
0239   {
0240     std::cout << PHWHERE << "JetCalib::CreateNodeTree : DST Node missing, aborting!" << std::endl;
0241     return Fun4AllReturnCodes::ABORTRUN;
0242   }
0243   PHCompositeNode *antiktNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "ANTIKT"));
0244   if (!antiktNode)
0245   {
0246     std::cout << PHWHERE << "JetCalib::CreateNodeTree : ANTIKT node missing, aborting!" << std::endl;
0247     return Fun4AllReturnCodes::ABORTRUN;
0248   }
0249   PHCompositeNode *towerNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "TOWER"));
0250   if (!towerNode)
0251   {
0252     std::cout << PHWHERE << "TOWER node not found, aborting!" << std::endl;
0253     return Fun4AllReturnCodes::ABORTRUN;
0254   }
0255
0256   // Create calib jet node if it doesn't exist.
0257   JetContainer *test_jets = findNode::getClass<JetContainer>(topNode, m_outputNode);
0258   if (!test_jets)
0259   {
0260     JetContainer *calib_jets = new JetContainerv1();
0261     PHIODataNode<PHObject> *calibjetNode;
0262     if (Verbosity() > 0)
0263     {
0264       std::cout << "JetCalib::CreateNode : creating " << m_outputNode << std::endl;
0265     }
0266     calibjetNode = new PHIODataNode<PHObject>(calib_jets, m_outputNode, "PHObject");
0267     towerNode->addNode(calibjetNode);
0268   }
0269   else
0270   {
0271     std::cout << "JetCalib::CreateNode : " << m_outputNode << " already exists! " << std::endl;
0272   }
0273
0274   return Fun4AllReturnCodes::EVENT_OK;
0275 }
0276
0277 std::string JetCalib::fetchCalibDir(const char *calibType)
0278 {
0279   std::string calibName = std::string("JES_Calib_") + calibType;
0280   return CDBInterface::instance()->getUrl(calibName);
0281 }
0282
0283 int getZvrtxBin(float zvrtx)
0284 {
0285   if (zvrtx >= -60.0 && zvrtx < -30.0)
0286   {
0287     return 1;  // -60 to -30
0288   }
0289   else if (zvrtx >= -30.0 && zvrtx < 30.0)
0290   {
0291     return 0;  // -30 to 30
0292   }
0293   else if (zvrtx >= 30.0 && zvrtx < 60)
0294   {
0295     return 2;  // 30 to 60
0296   }
0297   else if ((zvrtx < -60.0 && zvrtx >= -900) || (zvrtx >= 60.0 && zvrtx < 900))
0298   {
0299     return 3;  // -inf to -60 or 60 to inf
0300   }
0301   else if (zvrtx < -900)
0302   {
0303     return 4;  // -999 no zvertex
0304   }
0305
0306   return 0;  // Default case, should not happen
0307 }
0308
0309 int getEtaBin(int zvrtxbin, float eta, float jet_radius)
0310 {
0311   float eta_low = 0;
0312   float eta_high = 0;
0313   if (zvrtxbin == 0)  // -30 < zvrtx < 30
0314   {
0315     eta_low = -1.2 + jet_radius;
0316     eta_high = 1.2 - jet_radius;
0317   }
0318   else if (zvrtxbin == 1)  // -60 < zvrtx < -30
0319   {
0320     eta_low = -0.95 + jet_radius;
0321     eta_high = 1.25 - jet_radius;
0322   }
0323   else if (zvrtxbin == 2)  // 30 < zvrtx < inf
0324   {
0325     eta_low = -1.25 + jet_radius;
0326     eta_high = 0.95 - jet_radius;
0327   }
0328   else if (zvrtxbin == 3 || zvrtxbin == 4)
0329   {
0330     return 0;
0331   }
0332
0333   float threshold1 = eta_low + ((eta_high - eta_low) / 4.0);
0334   float threshold2 = eta_low + ((eta_high - eta_low) / 2.0);
0335   float threshold3 = eta_low + (3.0 * (eta_high - eta_low) / 4.0);
0336
0337   if (eta < threshold1)
0338   {
0339     return 0;  // -inf to threshold1
0340   }
0341   if (eta < threshold2)
0342   {
0343     return 1;  // threshold1 to threshold2
0344   }
0345   if (eta < threshold3)
0346   {
0347     return 2;  // threshold2 to threshold3
0348   }
0349
0350   return 3;  // threshold3 to inf
0351 }
0352
0353 float JetCalib::doCalibration(const std::vector<std::vector<TF1 *>> &JetCalibFunc, float jetPt, float zvrtx, float eta) const
0354 {
0355   float calib = 1;
0356   int zvrtxbin = 0;
0357   int etabin = 0;
0358   if (ApplyZvrtxDependentCalib || ApplyEtaDependentCalib)
0359   {
0360     zvrtxbin = getZvrtxBin(zvrtx);
0361     if (ApplyEtaDependentCalib)
0362     {
0363       etabin = getEtaBin(zvrtxbin, eta, jet_radius);
0364     }
0365   }
0366   calib = JetCalibFunc[zvrtxbin][etabin]->Eval(jetPt);
0367   return calib;
0368 }