g4simulation/g4bbc/MbdVertexFastSimReco.cc

0001
0002 #include "MbdVertexFastSimReco.h"
0003
0004 #include <globalvertex/MbdVertexMapv1.h>
0005 #include <globalvertex/MbdVertexv2.h>
0006
0007 #include <g4main/PHG4TruthInfoContainer.h>
0008 #include <g4main/PHG4VtxPoint.h>
0009
0010 #include <fun4all/Fun4AllReturnCodes.h>
0011 #include <fun4all/SubsysReco.h>  // for SubsysReco
0012
0013 #include <phool/PHCompositeNode.h>
0014 #include <phool/PHIODataNode.h>
0015 #include <phool/PHNode.h>  // for PHNode
0016 #include <phool/PHNodeIterator.h>
0017 #include <phool/PHObject.h>  // for PHObject
0018 #include <phool/PHRandomSeed.h>
0019 #include <phool/getClass.h>
0020 #include <phool/phool.h>  // for PHWHERE
0021
0022 #include <gsl/gsl_randist.h>
0023 #include <gsl/gsl_rng.h>  // for gsl_rng_uniform_pos, gsl_...
0024
0025 #include <cmath>
0026 #include <cstdlib>  // for exit
0027 #include <iostream>
0028 #include <limits>
0029
0030 MbdVertexFastSimReco::MbdVertexFastSimReco(const std::string &name)
0031   : SubsysReco(name)
0032   , RandomGenerator(gsl_rng_alloc(gsl_rng_mt19937))
0033 {
0034 }
0035
0036 MbdVertexFastSimReco::~MbdVertexFastSimReco()
0037 {
0038   gsl_rng_free(RandomGenerator);
0039 }
0040
0041 int MbdVertexFastSimReco::InitRun(PHCompositeNode *topNode)
0042 {
0043   if (std::isnan(m_T_Smear) || std::isnan(m_Z_Smear))
0044   {
0045     std::cout << PHWHERE << "::ERROR - smearing must be defined via setm_T_Smearing(float) and set_z_smeaering(float)" << std::endl;
0046     exit(-1);
0047   }
0048
0049   unsigned int seed = PHRandomSeed();  // fixed random seed handled in PHRandomSeed()
0050   gsl_rng_set(RandomGenerator, seed);
0051
0052   if (Verbosity() > 0)
0053   {
0054     std::cout << "===================== MbdVertexFastSimReco::InitRun() =====================" << std::endl;
0055     std::cout << " t smearing: " << m_T_Smear << " cm " << std::endl;
0056     std::cout << "  z smearing: " << m_Z_Smear << " cm " << std::endl;
0057     std::cout << " random seed: " << seed << std::endl;
0058     std::cout << "===========================================================================" << std::endl;
0059   }
0060
0061   return CreateNodes(topNode);
0062 }
0063
0064 int MbdVertexFastSimReco::process_event(PHCompositeNode *topNode)
0065 {
0066   if (Verbosity() > 1)
0067   {
0068     std::cout << "MbdVertexFastSimReco::process_event -- entered" << std::endl;
0069   }
0070
0071   //---------------------------------
0072   // Get Objects off of the Node Tree
0073   //---------------------------------
0074   PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0075   if (!truthinfo)
0076   {
0077     std::cout << PHWHERE << "::ERROR - cannot find G4TruthInfo" << std::endl;
0078     exit(-1);
0079   }
0080
0081   MbdVertexMap *vertexes = findNode::getClass<MbdVertexMap>(topNode, "MbdVertexMap");
0082   if (!vertexes)
0083   {
0084     std::cout << PHWHERE << "::ERROR - cannot find MbdVertexMap" << std::endl;
0085     exit(-1);
0086   }
0087
0088   //---------------------
0089   // Smear Truth Vertexes (only one per crossing right now)
0090   //---------------------
0091
0092   PHG4VtxPoint *point = truthinfo->GetPrimaryVtx(truthinfo->GetPrimaryVertexIndex());
0093   if (!point)
0094   {
0095     return Fun4AllReturnCodes::EVENT_OK;
0096   }
0097
0098   MbdVertex *vertex = new MbdVertexv2();
0099
0100   if (m_T_Smear >= 0.0)
0101   {
0102     vertex->set_t(point->get_t() + gsl_ran_gaussian(RandomGenerator, m_T_Smear));
0103     vertex->set_t_err(m_T_Smear);
0104   }
0105   else
0106   {
0107     vertex->set_t(point->get_t() + 2.0 * gsl_rng_uniform_pos(RandomGenerator) * m_T_Smear);
0108     vertex->set_t_err(std::abs(m_T_Smear) / sqrt(12.));
0109   }
0110
0111   if (m_Z_Smear >= 0.0)
0112   {
0113     vertex->set_z(point->get_z() + gsl_ran_gaussian(RandomGenerator, m_Z_Smear));
0114     vertex->set_z_err(m_Z_Smear);
0115   }
0116   else
0117   {
0118     vertex->set_z(point->get_z() + 2.0 * gsl_rng_uniform_pos(RandomGenerator) * m_Z_Smear);
0119     vertex->set_z_err(std::abs(m_Z_Smear) / sqrt(12.));
0120   }
0121
0122   vertexes->insert(vertex);
0123
0124   return Fun4AllReturnCodes::EVENT_OK;
0125 }
0126
0127 int MbdVertexFastSimReco::CreateNodes(PHCompositeNode *topNode)
0128 {
0129   PHNodeIterator iter(topNode);
0130
0131   // Looking for the DST node
0132   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0133   if (!dstNode)
0134   {
0135     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0136     return Fun4AllReturnCodes::ABORTRUN;
0137   }
0138
0139   // store the BBC stuff under a sub-node directory
0140   PHCompositeNode *bbcNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "BBC"));
0141   if (!bbcNode)
0142   {
0143     bbcNode = new PHCompositeNode("BBC");
0144     dstNode->addNode(bbcNode);
0145   }
0146
0147   // create the MbdVertexMap
0148   MbdVertexMap *vertexes = findNode::getClass<MbdVertexMap>(topNode, "MbdVertexMap");
0149   if (!vertexes)
0150   {
0151     vertexes = new MbdVertexMapv1();
0152     PHIODataNode<PHObject> *VertexMapNode = new PHIODataNode<PHObject>(vertexes, "MbdVertexMap", "PHObject");
0153     bbcNode->addNode(VertexMapNode);
0154   }
0155   else
0156   {
0157     std::cout << PHWHERE << "::ERROR - MbdVertexMap pre-exists, but should not if running FastSim" << std::endl;
0158     exit(-1);
0159   }
0160
0161   return Fun4AllReturnCodes::EVENT_OK;
0162 }