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