---

mMuiResponse.cxx

Go to the documentation of this file.

//
// Utility class: mMuiResponse:
// Author: C.Zhang 
// Date: 7/17/03
// Description: 
//              

// MUTOO/MUIOO headers
//
#include <mMuiResponse.h>
#include <mMuiResponsePar.h>
#include <TMuiMCHitMapO.h>
#include <TMutMCTrkMap.h>
#include <TMuiHitMapO.h>
#include <TMutNode.h>
#include <PHException.h>
#include <MUIOO.h>
#include <PHTimer.h>
#include <TMutGeo.h>
#include <TMuiHVMask.h>

//Mut Geometry


//Mut Calibration

// STL/BOOST/GSL
//
#include <iostream>
#include <string>
#include <gsl/gsl_randist.h>

mMuiResponse::mMuiResponse() :   
  _rng(0),
  _timer( PHTimeServer::get()->insert_new( "mMuiResponse" ) )
{
  name = "mMuiResponse";
  MUIOO::TRACE("initializing module " + std::string(name.getString()));  
}

// Destructor

mMuiResponse::~mMuiResponse() { }

// Event method.
//
PHBoolean mMuiResponse::event(PHCompositeNode* top_node)
{
  _timer.get()->restart();   
  try { 
    // Reset IOC pointers
    //
    set_interface_ptrs(top_node);        
    //   Function scope static gets initialized once
    //   Placed here instead of constructor so initialization is done
    //   only if module is used.
    //
    static bool rng_initialized = initialize_rng();
    if(rng_initialized != true) throw std::runtime_error(DESCRIPTION("mMuiResponse::event: failed to initialize rng"));
    // Loop over TMuiMCHits and generate TMuiHits.
    // Add noise and apply masking where appropriate.
    //
    response_loop();

  } catch(std::exception& e) {
    MUIOO::TRACE(e.what());
    return False;
  }    
  // If verbose dump the contents of the hit map
  //
  _timer.get()->stop();
  if(_mod_par->get_verbosity() >= MUIOO::ALOT) _hit_map->print();
  if(_mod_par->get_verbosity() >= MUIOO::SOME) _timer.get()->print();
  return True;
}

void 
mMuiResponse::set_interface_ptrs(PHCompositeNode* top_node)
{  
  // Module runtime parameters
  //
  _mod_par = TMutNode<mMuiResponsePar>::find_node(top_node,"mMuiResponsePar");

  // IOC 
  //
  _mc_hit_map = TMutNode<TMuiMCHitMapO>::find_node(top_node,"TMuiMCHitMapO");
  _mc_trk_map = TMutNode<TMutMCTrkMap>::find_node(top_node,"TMutMCTrkMap");
  _hit_map = TMutNode<TMuiHitMapO>::find_node(top_node,"TMuiHitMapO");
} 

void
mMuiResponse::response_loop()
{
  // Get an iterator to all the TMuiMCHits
  //
  TMuiMCHitMapO::iterator mc_hit_iter = _mc_hit_map->range();
  // Loop over TMuiMCHits
  //
  while(TMuiMCHitMapO::pointer mc_hit_ptr = mc_hit_iter.next()){    
    // fill hit map and associate it with MC hit.
    //
    fill_hit(mc_hit_ptr);
  }
}

bool 
mMuiResponse::is_alive(TMuiMCHitMapO::const_pointer mc_hit_ptr) 
{
  // TBI, some kind of code to check if the hit is in dead
  // area.
  //  .......................
  // return flase

  // Kill the hit if a random number is larger than the eff.
  //
  double x = gsl_ran_flat(_rng,0.0,100.0);
  if (x > _mod_par->get_twopack_eff()) return false;
  if(_mod_par->get_verbosity() >= MUIOO::ALOT) {
    cout <<" ######## check gsl_ran in mMuiRespose ######### " << endl;
    cout <<" # gsl_ran_uniform = " << x << " in mMuiRespose #" << endl;
    cout <<" # twopack efficiency =  "<<_mod_par->get_twopack_eff() <<"  #" << endl;
    cout <<" #####################################################" <<endl;
  }
  return true;
}

void
mMuiResponse::fill_hit(TMuiMCHitMapO::pointer mc_hit_ptr)
{
  UShort_t arm = mc_hit_ptr->get()->get_arm();
  UShort_t plane = mc_hit_ptr->get()->get_plane();
  
  typedef TMuiMCHitO::twopack_list twopack_list;
  typedef TMuiMCHitO::twopack_list::iterator twopack_iterator;
  
  // Get an iterator for list of TMuiMCTwoPack from TMuiMCHit
  //
  twopack_list* twopacks = mc_hit_ptr->get()->get_mutable_twopack_list();
  twopack_iterator twopack_iter = twopacks->begin();
  
  for(;twopack_iter!=twopacks->end();++twopack_iter){
    
    // If TMuiHit exists then do nothing
    // if not create a new TMuiHit
    //
    UShort_t panel       =  twopack_iter->get_panel();
    UShort_t orientation =  twopack_iter->get_orient();
    UShort_t twopack     =  twopack_iter->get_twopack_index();

    // Check the HV and set the two pack status if this hit
    // is masked
    //
    if(_mod_par->get_use_hv_mask() && 
       !check_efficiency(arm,plane,panel,orientation,twopack)){
      twopack_iter->set_status(TMuiMCTwoPack::MASKED);
      continue;
    }
    TMuiHitMapO::iterator hit_iter = _hit_map->get(arm,
                                                   plane,
                                                   panel,
                                                   orientation,
                                                   twopack);

    // Look for existing hit on this strip, if none then
    // make a new one.
    //
    if(hit_iter.at_end()){
       hit_iter = _hit_map->insert_new(arm,
                                       plane,
                                       panel,
                                       orientation,
                                       twopack);

    }

    associate_mc(hit_iter, mc_hit_ptr);
  }

}

void
mMuiResponse::associate_mc(TMuiHitMapO::iterator hit_iter, 
                           TMuiMCHitMapO::const_pointer mc_hit_ptr)
{
  TMuiHitMapO::pointer hit_ptr = hit_iter.current();
  PHKey::associate(TMuiMCHitMapO::pointer(mc_hit_ptr), hit_ptr); 
}


bool
mMuiResponse::initialize_rng()
{
  const gsl_rng_type *T;
  gsl_rng_env_setup();
  T = gsl_rng_default;
  _rng = gsl_rng_alloc(T);

  return true;
}


bool 
mMuiResponse::check_efficiency(UShort_t arm, 
                               UShort_t plane, 
                               UShort_t panel, 
                               UShort_t orientation, 
                               UShort_t twopack) const
{    
  // Efficiency due to HV Mask
  //
  if(_mod_par->get_use_hv_mask()){
    // Check HV state
    //
    double random = gsl_ran_flat(_rng,0.0,1.0);
    double effic = TMuiHVMask::get_effic_twopack(arm,plane,panel,orientation,twopack);      
    if (random > effic) return false;
  }
  
  // Efficiency due to everything else
  //
  double random = gsl_ran_flat(_rng,0.0,1.0);
  double effic = _mod_par->get_twopack_eff();
  if (random > effic) return false;
  
  return true;
}

---