g4simulation/g4vertex/GlobalVertexFastSimReco.cc

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