RndmEmbed.cxx

Go to the documentation of this file.

/*
 * RndmEmbed.C
 *
 *  Created on: Oct 29, 2017
 *      Author: yuhw@nmsu.edu
 */
 
 
#include "RndmEmbed.h"
 
#include <interface_main/SQHit.h>
#include <interface_main/SQHit_v1.h>
#include <interface_main/SQMCHit_v1.h>
#include <interface_main/SQHitMap_v1.h>
#include <interface_main/SQHitVector_v1.h>
#include <interface_main/SQEvent_v1.h>
#include <interface_main/SQRun_v1.h>
#include <interface_main/SQSpill_v1.h>
#include <interface_main/SQSpillMap_v1.h>
 
#include <ktracker/SRecEvent.h>
 
#include <geom_svc/GeomSvc.h>
 
#include <fun4all/Fun4AllReturnCodes.h>
#include <fun4all/PHTFileServer.h>
#include <phool/PHNodeIterator.h>
#include <phool/PHIODataNode.h>
#include <phool/getClass.h>
#include <phool/PHRandomSeed.h>
 
#include <g4main/PHG4TruthInfoContainer.h>
#include <g4main/PHG4HitContainer.h>
#include <g4main/PHG4Hit.h>
#include <g4main/PHG4Particle.h>
#include <g4main/PHG4HitDefs.h>
#include <g4main/PHG4VtxPoint.h>
 
#include <TFile.h>
#include <TTree.h>
 
#include <cstring>
#include <cmath>
#include <cfloat>
#include <stdexcept>
#include <limits>
#include <tuple>
#include <memory>
 
// gsl
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
 
#include <boost/lexical_cast.hpp>
 
using namespace std;
 
RndmEmbed::RndmEmbed(const std::string& name) :
SubsysReco(name),
_hit_container_type("Vector"),
_noise_rate(),
_event(0),
_hit_map(nullptr),
_hit_vector(nullptr),
_out_name("RndmEmbedEval.root")
{
}
 
int RndmEmbed::Init(PHCompositeNode* topNode) {
        return Fun4AllReturnCodes::EVENT_OK;
}
 
int RndmEmbed::InitRun(PHCompositeNode* topNode) {
 
        ResetEvalVars();
        InitEvalTree();
 
        p_geomSvc = GeomSvc::instance();
 
        int ret = GetNodes(topNode);
        if(ret != Fun4AllReturnCodes::EVENT_OK) return ret;
 
        return Fun4AllReturnCodes::EVENT_OK;
}
 
int RndmEmbed::process_event(PHCompositeNode* topNode) {
 
        if(Verbosity() >= Fun4AllBase::VERBOSITY_SOME)
                std::cout << "Entering RndmEmbed::process_event: " << _event << std::endl;
 
        ResetEvalVars();
 
  if(!_hit_vector and _hit_container_type.find("Vector") == std::string::npos) {
        LogInfo("!_hit_vector and _hit_container_type.find(\"Vector\") == std::string::npos");
        return Fun4AllReturnCodes::ABORTRUN;
  }
 
#ifndef __CINT__
 
  gsl_rng* RandomGenerator = gsl_rng_alloc(gsl_rng_mt19937);
  unsigned int seed = PHRandomSeed();  // fixed seed is handled in this funtcion
  //seed = 128;
        //std::cout << Name() << " random seed: " << seed << std::endl;
  gsl_rng_set(RandomGenerator, seed);
 
  for (auto i : _noise_rate) {
        auto detector_name = i.first;
        auto noise_rate = i.second;
 
        auto detector_id = p_geomSvc->getDetectorID(detector_name);
        if(detector_id==0) {
                LogInfo("detector_id==0");
                continue;
        }
 
        Plane plane = p_geomSvc->getPlane(detector_id);
 
        for( int element_id=1; element_id <= plane.nElements; ++element_id) {
                bool acc = gsl_ran_flat(RandomGenerator, 0, 1) < noise_rate ? true : false;
                if(!acc) continue;
 
                {
                                double drift = plane.spacing * gsl_ran_flat(RandomGenerator, -1, 1);
 
                                auto up_digiHit = std::unique_ptr<SQMCHit_v1> (new SQMCHit_v1());
                                auto digiHit = up_digiHit.get();
 
                                digiHit->set_hit_id(_hit_vector->size());
 
                                digiHit->set_detector_id(detector_id);
                                digiHit->set_element_id(element_id);
                                digiHit->set_drift_distance(drift);
 
                                digiHit->set_pos(p_geomSvc->getMeasurement(digiHit->get_detector_id(), digiHit->get_element_id()));
 
                                // FIXME figure this out
                                digiHit->set_in_time(1);
                                digiHit->set_hodo_mask(1);
 
                                digiHit->set_track_id(std::numeric_limits<int>::max());
                                digiHit->set_g4hit_id(std::numeric_limits<int>::max());
 
                                digiHit->set_truth_x(std::numeric_limits<float>::max());
                                digiHit->set_truth_y(std::numeric_limits<float>::max());
                                digiHit->set_truth_z(std::numeric_limits<float>::max());
 
                                _hit_vector->push_back(digiHit);
                }
 
                {
                                double drift = plane.spacing * gsl_ran_flat(RandomGenerator, -1, 1);
 
                                auto up_digiHit = std::unique_ptr<SQMCHit_v1> (new SQMCHit_v1());
                                auto digiHit = up_digiHit.get();
 
                                digiHit->set_hit_id(_hit_vector->size());
 
                                digiHit->set_detector_id(detector_id+1);
                                digiHit->set_element_id(element_id);
                                digiHit->set_drift_distance(drift);
 
                                digiHit->set_pos(p_geomSvc->getMeasurement(digiHit->get_detector_id(), digiHit->get_element_id()));
 
                                // FIXME figure this out
                                digiHit->set_in_time(1);
                                digiHit->set_hodo_mask(1);
 
                                digiHit->set_track_id(std::numeric_limits<int>::max());
                                digiHit->set_g4hit_id(std::numeric_limits<int>::max());
 
                                digiHit->set_truth_x(std::numeric_limits<float>::max());
                                digiHit->set_truth_y(std::numeric_limits<float>::max());
                                digiHit->set_truth_z(std::numeric_limits<float>::max());
 
                                _hit_vector->push_back(digiHit);
                }
 
        }
  }
 
  gsl_rng_free(RandomGenerator);
 
#endif //__CINT__
 
  _tout->Fill();
 
  if(Verbosity() >= Fun4AllBase::VERBOSITY_SOME)
    std::cout << "Leaving RndmEmbed::process_event: " << _event << std::endl;
  ++_event;
 
  return Fun4AllReturnCodes::EVENT_OK;
}
 
int RndmEmbed::End(PHCompositeNode* topNode) {
  if(Verbosity() >= Fun4AllBase::VERBOSITY_SOME)
    std::cout << "RndmEmbed::End" << std::endl;
 
  PHTFileServer::get().cd(_out_name.c_str());
  _tout->Write();
 
  return Fun4AllReturnCodes::EVENT_OK;
}
 
int RndmEmbed::InitEvalTree() {
  PHTFileServer::get().open(_out_name.c_str(), "RECREATE");
 
  _tout = new TTree("T", "RndmEmbed");
 
  return 0;
}
 
int RndmEmbed::ResetEvalVars() {
  return 0;
}
 
int RndmEmbed::GetNodes(PHCompositeNode* topNode) {
 
  if(_hit_container_type.find("Map") != std::string::npos) {
    _hit_map = findNode::getClass<SQHitMap>(topNode, "SQHitMap");
    if (!_hit_map) {
      LogInfo("!_hit_map");
      return Fun4AllReturnCodes::ABORTEVENT;
    }
  }
 
  if(_hit_container_type.find("Vector") != std::string::npos) {
    _hit_vector = findNode::getClass<SQHitVector>(topNode, "SQHitVector");
    if (!_hit_vector) {
        LogInfo("!_hit_vector");
      return Fun4AllReturnCodes::ABORTEVENT;
    }
  }
 
  return Fun4AllReturnCodes::EVENT_OK;
}