packages/trigger/MinimumBiasClassifier.cc

0001 #include "MinimumBiasClassifier.h"
0002
0003 #include "MinimumBiasInfov1.h"
0004
0005 #include <zdcinfo/Zdcinfo.h>
0006
0007 #include <cdbobjects/CDBTTree.h>
0008 #include <ffamodules/CDBInterface.h>
0009
0010 #include <globalvertex/GlobalVertex.h>
0011 #include <globalvertex/GlobalVertexMap.h>
0012 #include <mbd/MbdPmtContainer.h>
0013 #include <mbd/MbdPmtHit.h>
0014
0015 #include <fun4all/Fun4AllReturnCodes.h>
0016
0017 #include <phool/PHCompositeNode.h>
0018 #include <phool/PHIODataNode.h>
0019 #include <phool/PHNode.h>
0020 #include <phool/PHNodeIterator.h>
0021 #include <phool/PHObject.h>
0022 #include <phool/getClass.h>
0023 #include <phool/phool.h>
0024
0025 #include <cmath>
0026 #include <filesystem>
0027 #include <iostream>
0028 #include <map>  // for _Rb_tree_iterator
0029 #include <string>
0030 #include <utility>  // for pair
0031
0032 MinimumBiasClassifier::MinimumBiasClassifier(const std::string &name)
0033   : SubsysReco(name)
0034 {
0035 }
0036 int MinimumBiasClassifier::InitRun(PHCompositeNode *topNode)
0037 {
0038   if (Verbosity() > 1)
0039   {
0040     std::cout << __FILE__ << " :: " << __FUNCTION__ << std::endl;
0041   }
0042   CDBInterface *m_cdb = CDBInterface::instance();
0043
0044   std::string centscale_url = m_cdb->getUrl("CentralityScale");
0045   if (m_overwrite_scale)
0046   {
0047     centscale_url = m_overwrite_url_scale;
0048     std::cout << " Overwriting Scale to " << m_overwrite_url_scale << std::endl;
0049   }
0050
0051   if (Download_centralityScale(centscale_url))
0052   {
0053     return Fun4AllReturnCodes::ABORTRUN;
0054   }
0055
0056   std::string vertexscale_url = m_cdb->getUrl("CentralityVertexScale");
0057   if (m_overwrite_vtx)
0058   {
0059     vertexscale_url = m_overwrite_url_vtx;
0060     std::cout << " Overwriting Vtx to " << m_overwrite_url_vtx << std::endl;
0061   }
0062
0063   if (Download_centralityVertexScales(vertexscale_url))
0064   {
0065     return Fun4AllReturnCodes::ABORTRUN;
0066   }
0067
0068   CreateNodes(topNode);
0069
0070   m_zdc_energy_sum.fill(0);
0071   m_mbd_charge_sum.fill(0);
0072   m_mbd_hit.fill(0);
0073
0074   return Fun4AllReturnCodes::EVENT_OK;
0075 }
0076
0077 int MinimumBiasClassifier::ResetEvent(PHCompositeNode * /*unused*/)
0078 {
0079   m_zdc_energy_sum.fill(0);
0080   m_mbd_charge_sum.fill(0);
0081   m_mbd_hit.fill(0);
0082
0083   return Fun4AllReturnCodes::EVENT_OK;
0084 }
0085
0086 int MinimumBiasClassifier::FillMinimumBiasInfo()
0087 {
0088   if (Verbosity())
0089   {
0090     std::cout << "Getting Vertex" << std::endl;
0091   }
0092
0093   // if (!m_global_vertex_map)
0094   // {
0095   //   m_mb_info->setIsAuAuMinimumBias(false);
0096   //   return Fun4AllReturnCodes::EVENT_OK;
0097   // }
0098
0099   // if (m_global_vertex_map->empty())
0100   // {
0101   //   m_mb_info->setIsAuAuMinimumBias(false);
0102   //   return Fun4AllReturnCodes::EVENT_OK;
0103   // }
0104
0105   if (m_global_vertex_map->empty())
0106   {
0107     m_mb_info->setIsAuAuMinimumBias(false);
0108     return Fun4AllReturnCodes::EVENT_OK;
0109   }
0110
0111   GlobalVertex *vtx = m_global_vertex_map->begin()->second;
0112   if (!vtx)
0113   {
0114     m_mb_info->setIsAuAuMinimumBias(false);
0115     return Fun4AllReturnCodes::EVENT_OK;
0116   }
0117
0118   if (!vtx->isValid())
0119   {
0120     m_mb_info->setIsAuAuMinimumBias(false);
0121     return Fun4AllReturnCodes::EVENT_OK;
0122   }
0123
0124   bool minbiascheck = true;
0125   ;
0126
0127   m_vertex = vtx->get_z();
0128
0129   m_vertex_scale = getVertexScale();
0130
0131   if (Verbosity())
0132   {
0133     std::cout << "Getting ZDC" << std::endl;
0134   }
0135   if (!m_issim)
0136   {
0137     if (!m_zdcinfo)
0138     {
0139       m_mb_info->setIsAuAuMinimumBias(false);
0140       return Fun4AllReturnCodes::EVENT_OK;
0141     }
0142   }
0143   //  Z vertex is within range
0144   if (std::fabs(m_vertex) > m_z_vtx_cut && minbiascheck)
0145   {
0146     minbiascheck = false;
0147   }
0148   if (Verbosity())
0149   {
0150     std::cout << "Calculating" << std::endl;
0151   }
0152
0153   // calculate charge sum and n hit
0154   for (int i = 0; i < 128; i++)
0155   {
0156     m_mbd_pmt = m_mbd_container->get_pmt(i);
0157     short side = i / 64;
0158     bool pass = passesHitCut(m_mbd_pmt);
0159     if (!pass)
0160     {
0161       continue;
0162     }
0163     m_mbd_hit[side]++;
0164     m_mbd_charge_sum[side] += m_mbd_pmt->get_q() * m_vertex_scale * m_centrality_scale;
0165   }
0166   if (Verbosity())
0167   {
0168     std::cout << m_mbd_charge_sum[0] << " " << m_mbd_charge_sum[1] << std::endl;
0169   }
0170
0171   // MBD Background cut
0172   if (m_mbd_charge_sum[1] < m_mbd_north_cut && m_mbd_charge_sum[0] > m_mbd_south_cut && minbiascheck)
0173   {
0174     minbiascheck = false;
0175     //    m_mb_info->setIsAuAuMinimumBias(false);
0176     // return Fun4AllReturnCodes::EVENT_OK;
0177   }
0178
0179   // Mbd two hit requirement and ZDC energy sum coincidence requirement
0180   for (int iside = 0; iside < 2; iside++)
0181   {
0182     if (m_mbd_hit[iside] < 2 && minbiascheck)
0183     {
0184       minbiascheck = false;
0185       // m_mb_info->setIsAuAuMinimumBias(false);
0186       // return Fun4AllReturnCodes::EVENT_OK;
0187     }
0188     if (!m_issim)
0189     {
0190       if (m_zdcinfo->get_zdc_energy(iside) <= m_zdc_cut && minbiascheck)
0191       {
0192         minbiascheck = false;
0193         // m_mb_info->setIsAuAuMinimumBias(false);
0194         // return Fun4AllReturnCodes::EVENT_OK;
0195       }
0196     }
0197   }
0198   if ((m_mbd_charge_sum[0] + m_mbd_charge_sum[1]) > 2100 && minbiascheck)
0199   {
0200     minbiascheck = false;
0201     // m_mb_info->setIsAuAuMinimumBias(false);
0202     // return Fun4AllReturnCodes::EVENT_OK;
0203   }
0204
0205   m_mb_info->setIsAuAuMinimumBias(minbiascheck);
0206
0207   return Fun4AllReturnCodes::EVENT_OK;
0208 }
0209 int MinimumBiasClassifier::process_event(PHCompositeNode *topNode)
0210 {
0211   if (Verbosity())
0212   {
0213     std::cout << "------------MinimumBiasClassifier-------------" << std::endl;
0214   }
0215
0216   // Get Nodes from the Tree
0217   if (GetNodes(topNode))
0218   {
0219     return Fun4AllReturnCodes::EVENT_OK;
0220   }
0221
0222   if (FillMinimumBiasInfo())
0223   {
0224     return Fun4AllReturnCodes::EVENT_OK;
0225   }
0226
0227   return Fun4AllReturnCodes::EVENT_OK;
0228 }
0229
0230 int MinimumBiasClassifier::GetNodes(PHCompositeNode *topNode)
0231 {
0232   if (Verbosity() > 1)
0233   {
0234     std::cout << __FILE__ << " :: " << __FUNCTION__ << " :: " << __LINE__ << std::endl;
0235   }
0236
0237   m_mb_info = findNode::getClass<MinimumBiasInfo>(topNode, "MinimumBiasInfo");
0238
0239   if (!m_mb_info)
0240   {
0241     std::cout << "no minimum bias node " << std::endl;
0242     return Fun4AllReturnCodes::ABORTRUN;
0243   }
0244
0245   m_mbd_container = findNode::getClass<MbdPmtContainer>(topNode, "MbdPmtContainer");
0246   if (Verbosity())
0247   {
0248     std::cout << "Getting MBD Tubes" << std::endl;
0249   }
0250
0251   if (!m_mbd_container)
0252   {
0253     std::cout << "no MBD out node " << std::endl;
0254     return Fun4AllReturnCodes::ABORTRUN;
0255   }
0256
0257   if (!m_issim)
0258   {
0259     m_zdcinfo = findNode::getClass<Zdcinfo>(topNode, "Zdcinfo");
0260     if (Verbosity())
0261     {
0262       std::cout << "Getting ZDC Info" << std::endl;
0263     }
0264
0265     if (!m_zdcinfo)
0266     {
0267       std::cout << "no zdc towers node " << std::endl;
0268       return Fun4AllReturnCodes::ABORTRUN;
0269     }
0270   }
0271   if (Verbosity())
0272   {
0273     std::cout << "Getting Vertex Map" << std::endl;
0274   }
0275
0276   m_global_vertex_map = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0277
0278   if (!m_global_vertex_map)
0279   {
0280     std::cout << "no vertex map node " << std::endl;
0281     return Fun4AllReturnCodes::ABORTRUN;
0282   }
0283   return Fun4AllReturnCodes::EVENT_OK;
0284 }
0285
0286 void MinimumBiasClassifier::CreateNodes(PHCompositeNode *topNode)
0287 {
0288   PHNodeIterator iter(topNode);
0289   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0290   if (!dstNode)
0291   {
0292     std::cout << PHWHERE << "DST Node missing doing nothing" << std::endl;
0293   }
0294
0295   PHNodeIterator dstIter(dstNode);
0296
0297   PHCompositeNode *detNode = dynamic_cast<PHCompositeNode *>(dstIter.findFirst("PHCompositeNode", "GLOBAL"));
0298   if (!detNode)
0299   {
0300     std::cout << PHWHERE << "Detector Node missing, making one" << std::endl;
0301     detNode = new PHCompositeNode("GLOBAL");
0302     dstNode->addNode(detNode);
0303   }
0304
0305   MinimumBiasInfo *mb = new MinimumBiasInfov1();
0306
0307   PHIODataNode<PHObject> *mbNode = new PHIODataNode<PHObject>(mb, "MinimumBiasInfo", "PHObject");
0308   detNode->addNode(mbNode);
0309
0310   return;
0311 }
0312
0313 bool MinimumBiasClassifier::passesHitCut(MbdPmtHit *hit) const
0314 {
0315   if (fabs(hit->get_time()) >= m_mbd_time_cut)
0316   {
0317     return false;
0318   }
0319   if (hit->get_q() <= m_mbd_charge_cut)
0320   {
0321     return false;
0322   }
0323
0324   return true;
0325 }
0326 int MinimumBiasClassifier::Download_centralityScale(const std::string &dbfile)
0327 {
0328   m_centrality_scale = 1.00;
0329
0330   std::filesystem::path dbase_file = dbfile;
0331   if (dbase_file.extension() == ".root")
0332   {
0333     CDBTTree *cdbttree = new CDBTTree(dbase_file);
0334     cdbttree->LoadCalibrations();
0335     m_centrality_scale = cdbttree->GetDoubleValue(0, "centralityscale");
0336     if (Verbosity())
0337     {
0338       std::cout << "centscale = " << m_centrality_scale << std::endl;
0339     }
0340     delete cdbttree;
0341   }
0342   else
0343   {
0344     std::cout << PHWHERE << ", ERROR, unknown file type, " << dbfile << std::endl;
0345     return Fun4AllReturnCodes::ABORTRUN;
0346   }
0347
0348   return Fun4AllReturnCodes::EVENT_OK;
0349 }
0350
0351 int MinimumBiasClassifier::Download_centralityVertexScales(const std::string &dbfile)
0352 {
0353   std::filesystem::path dbase_file = dbfile;
0354
0355   if (dbase_file.extension() == ".root")
0356   {
0357     CDBTTree *cdbttree = new CDBTTree(dbase_file);
0358     cdbttree->LoadCalibrations();
0359     if (Verbosity())
0360     {
0361       cdbttree->Print();
0362     }
0363
0364     int nvertexbins = cdbttree->GetIntValue(0, "nvertexbins");
0365
0366     for (int iv = 0; iv < nvertexbins; iv++)
0367     {
0368       float scale = cdbttree->GetDoubleValue(iv, "scale");
0369       float lowvertex = cdbttree->GetDoubleValue(iv, "low_vertex");
0370       float highvertex = cdbttree->GetDoubleValue(iv, "high_vertex");
0371       m_vertex_scales.emplace_back(std::make_pair(lowvertex, highvertex), scale);
0372     }
0373
0374     delete cdbttree;
0375   }
0376   else
0377   {
0378     std::cout << PHWHERE << ", ERROR, unknown file type, " << dbfile << std::endl;
0379     return Fun4AllReturnCodes::ABORTRUN;
0380   }
0381
0382   return Fun4AllReturnCodes::EVENT_OK;
0383 }
0384 float MinimumBiasClassifier::getVertexScale()
0385 {
0386   for (auto v_range_scale : m_vertex_scales)
0387   {
0388     auto v_range = v_range_scale.first;
0389     if (m_vertex > v_range.first && m_vertex <= v_range.second)
0390     {
0391       return v_range_scale.second;
0392     }
0393   }
0394   return 0;
0395 }