packages/jetbackground/TimingCut.cc

0001 #include "TimingCut.h"
0002
0003 #include <mbd/MbdOut.h>
0004
0005 #include <jetbase/Jet.h>
0006 #include <jetbase/JetContainer.h>
0007
0008 #include <phparameter/PHParameters.h>
0009
0010 #include <fun4all/Fun4AllReturnCodes.h>
0011
0012 #include <phool/PHCompositeNode.h>
0013 #include <phool/getClass.h>
0014
0015 #include <cmath>
0016 #include <iostream>  // for basic_ostream, operator<<
0017 #include <map>       // for _Rb_tree_iterator, opera...
0018 #include <utility>   // for pair
0019 #include <vector>    // for vector
0020 //____________________________________________________________________________..
0021 TimingCut::TimingCut(const std::string &jetNodeName, const std::string &name, const bool doAbort)
0022   : SubsysReco(name)
0023   , _doAbort(doAbort)
0024   , _jetNodeName(jetNodeName)
0025   , _cutParams(name)
0026 {
0027   SetDefaultParams();
0028 }
0029
0030 //____________________________________________________________________________..
0031 int TimingCut::Init(PHCompositeNode *topNode)
0032 {
0033   if(CreateNodeTree(topNode))
0034     {
0035       return Fun4AllReturnCodes::ABORTRUN;
0036     }
0037
0038   return Fun4AllReturnCodes::EVENT_OK;
0039 }
0040
0041 int TimingCut::CreateNodeTree(PHCompositeNode *topNode)
0042 {
0043   PHNodeIterator iter(topNode);
0044
0045   PHCompositeNode *parNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "PAR"));
0046   if (!parNode)
0047   {
0048     std::cout << "No RUN node found; cannot create PHParameters for storing cut results!";
0049     return 1;
0050   }
0051
0052   _cutParams.SaveToNodeTree(parNode, "TimingCutParams");
0053   return 0;
0054 }
0055
0056 //____________________________________________________________________________..
0057 int TimingCut::process_event(PHCompositeNode *topNode)
0058 {
0059   JetContainer *jets = findNode::getClass<JetContainer>(topNode, _jetNodeName);
0060   if (!jets)
0061   {
0062     if (Verbosity() > 0 && !_missingInfoWarningPrinted)
0063     {
0064       std::cout << "Missing jets; abort event. Further warnings will be suppressed." << std::endl;
0065     }
0066     _missingInfoWarningPrinted = true;
0067     return Fun4AllReturnCodes::ABORTEVENT;
0068   }
0069
0070   float maxJetpT = 0;
0071   float subJetpT = 0;
0072   float maxJett = std::numeric_limits<float>::quiet_NaN();
0073   float subJett = std::numeric_limits<float>::quiet_NaN();
0074   float maxJetPhi = std::numeric_limits<float>::quiet_NaN();
0075   float subJetPhi = std::numeric_limits<float>::quiet_NaN();
0076
0077   if (jets)
0078   {
0079     int tocheck = jets->size();
0080     if (Verbosity() > 2)
0081     {
0082       std::cout << "Found " << tocheck << " jets to check..." << std::endl;
0083     }
0084     for (int i = 0; i < tocheck; ++i)
0085     {
0086       float jetpT = 0;
0087       float jett = std::numeric_limits<float>::quiet_NaN();
0088       float jetPhi = std::numeric_limits<float>::quiet_NaN();
0089       Jet *jet = jets->get_jet(i);
0090       if (jet)
0091       {
0092         jetpT = jet->get_pt();
0093     jett = jet->get_property(Jet::PROPERTY::prop_t);
0094     jetPhi = jet->get_phi();
0095       }
0096       else
0097       {
0098         continue;
0099       }
0100       if (jetpT > maxJetpT)
0101       {
0102         if (maxJetpT)
0103         {
0104           subJetpT = maxJetpT;
0105           subJett = maxJett;
0106       subJetPhi = maxJetPhi;
0107         }
0108         maxJetpT = jetpT;
0109         maxJett = jett;
0110     maxJetPhi = jetPhi;
0111       }
0112       else if (jetpT > subJetpT)
0113       {
0114         subJetpT = jetpT;
0115         subJett = jett;
0116     subJetPhi = jetPhi;
0117       }
0118     }
0119   }
0120   else
0121   {
0122     if (Verbosity() > 0)
0123     {
0124       std::cout << "No jet node!" << std::endl;
0125     }
0126     return Fun4AllReturnCodes::ABORTEVENT;
0127   }
0128
0129   bool passDeltat = Pass_Delta_t(maxJett, subJett, maxJetPhi, subJetPhi);
0130   bool passLeadt = Pass_Lead_t(maxJett);
0131
0132   MbdOut * mbdout = static_cast<MbdOut*>(findNode::getClass<MbdOut>(topNode,"MbdOut"));
0133   float m_mbd_t0 = std::numeric_limits<float>::quiet_NaN();
0134   float m_mbd_ts = std::numeric_limits<float>::quiet_NaN();
0135   float m_mbd_tn = std::numeric_limits<float>::quiet_NaN();
0136   if(mbdout)
0137     {
0138       m_mbd_t0 = mbdout->get_t0();
0139       m_mbd_ts = mbdout->get_time(0); // south side
0140       m_mbd_tn = mbdout->get_time(1); // north side
0141     }
0142
0143   float mbd_time = std::numeric_limits<float>::quiet_NaN();
0144   if(!std::isnan(m_mbd_t0))
0145     {
0146       mbd_time = m_mbd_t0;
0147     }
0148   else if(!std::isnan(m_mbd_tn))
0149     {
0150       mbd_time = m_mbd_tn;
0151     }
0152   else if(!std::isnan(m_mbd_ts))
0153     {
0154       mbd_time = m_mbd_ts;
0155     }
0156
0157   bool passMbdt = false;
0158   if(!std::isnan(mbd_time))
0159     {
0160       passMbdt = Pass_Mbd_dt(maxJett, mbd_time);
0161     }
0162
0163   bool failAnyCut = !passDeltat || !passLeadt || !passMbdt;
0164
0165   if (failAnyCut && _doAbort)
0166   {
0167     return Fun4AllReturnCodes::ABORTEVENT;
0168   }
0169   PHNodeIterator iter(topNode);
0170   PHCompositeNode *parNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "PAR"));
0171   _cutParams.set_int_param("passLeadtCut",passLeadt);
0172   _cutParams.set_int_param("passDeltatCut",passDeltat);
0173   _cutParams.set_int_param("passMbdDtCut",passMbdt);
0174   _cutParams.set_int_param("failAnyTimeCut", failAnyCut);
0175   _cutParams.set_double_param("maxJett",maxJett);
0176   _cutParams.set_double_param("subJett",subJett);
0177   _cutParams.set_double_param("mbd_time",mbd_time);
0178   _cutParams.set_double_param("dPhi",calc_dphi(maxJetPhi, subJetPhi));
0179   _cutParams.UpdateNodeTree(parNode, "TimingCutParams");
0180
0181   return Fun4AllReturnCodes::EVENT_OK;
0182 }
0183 //____________________________________________________________________________..
0184 int TimingCut::ResetEvent(PHCompositeNode * /*topNode*/)
0185 {
0186   if (Verbosity() > 0)
0187   {
0188     std::cout << "TimingCut::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
0189   }
0190   return Fun4AllReturnCodes::EVENT_OK;
0191 }
0192
0193 //____________________________________________________________________________..
0194 int TimingCut::End(PHCompositeNode * /*topNode*/)
0195 {
0196   if (Verbosity() > 0)
0197   {
0198     std::cout << "TimingCut::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0199   }
0200   return Fun4AllReturnCodes::EVENT_OK;
0201 }
0202
0203 //____________________________________________________________________________..
0204 int TimingCut::Reset(PHCompositeNode * /*topNode*/)
0205 {
0206   if (Verbosity() > 0)
0207   {
0208     std::cout << "TimingCut::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0209   }
0210   return Fun4AllReturnCodes::EVENT_OK;
0211 }
0212
0213 //____________________________________________________________________________..
0214 void TimingCut::Print(const std::string &what) const
0215 {
0216   std::cout << "TimingCut::Print(const std::string &what) const Printing info for " << what << std::endl;
0217 }