---

mMuiSlowSim.cxx

Go to the documentation of this file.

//
// Utility class: mMuiSlowSim:
// Author: J.Nagle, 
// major changes made 9Jul03 by Chun Zhang 
// comments added by Donald Hornback
// Date: 6/23/03
// Description: Convert PISA to detector hits in MUID
//              

// MUIOO headers
#include <mMuiSlowSim.h>
#include <mMuiSlowSimPar.h>
#include <TMutNode.h>
#include <MUIOO.h>

// PISA hit table,
#include <root_ptrk.h>
#include <rootAnc.h>

// Geometry modules
#include <MuiGeomClasses.hh>

// PHENIX headers
#include <PHGeometryObject.h>
#include <PHVector.h>
#include <PHException.h>
#include <PHTimer.h>
//#include <gufld.h>

// STL/BOOST
//
#include <iostream>
#include <string>

mMuiSlowSim::mMuiSlowSim() : 
  _timer( PHTimeServer::get()->insert_new("mMuiSlowSim") )
{
  name = "mMuiSlowSim";
  MUTOO::TRACE("initializing module " + std::string(name.getString()));  
}

// Destructor

mMuiSlowSim::~mMuiSlowSim() { }

// Event method.

PHBoolean mMuiSlowSim::event(PHCompositeNode* top_node)
{

  _timer.get()->restart(); 
  
  try { 
    // Reset IOC pointers
    //
    set_interface_ptrs(top_node);    

    // Do the conversion of PISA hits to TMuiMCHit objects
    //
    digitize();

    // Do the association between MCMuiHit and MCTrack,
    // MuTr slowsim has to been run first.
    //
    associate_mctrk();

  } catch(std::exception& e) {
    // MODULE CODE GOES HERE
    //
    MUTOO::TRACE(e.what());
    return False;
  }  
  
  // If verbose dump the contents of the cluster map
  // Verbosity methods declared in TMuiParBase
  
  _timer.get()->stop();
  if(_mod_par->get_verbosity() >= MUIOO::ALOT) _mc_hit_map->print();
  if(_mod_par->get_verbosity() >= MUIOO::SOME) _timer.get()->print();     
  return True;
}

void 
mMuiSlowSim::set_interface_ptrs(PHCompositeNode* top_node)
{  
  // module runtime parameters
  //
  _mod_par = TMutNode<mMuiSlowSimPar>::find_node(top_node,"mMuiSlowSimPar");

  // Interface Object Container pointers for muid MC hits.
  //
  _mc_hit_map = TMutNode<TMuiMCHitMapO>::find_node(top_node,"TMuiMCHitMapO");

  // Interface Object Containter pointers for MC tracks.
  // 
  _mc_trk_map = TMutNode<TMutMCTrkMap>::find_node(top_node,"TMutMCTrkMap");

  // PISA input in wrapped STAF table form here.  Use of wrapped STAF tables
  // has been depreciated -- this interface will eventually go away.
  //
  munhitsWrapper* munhits_ptr= TMutNode<munhitsWrapper>::find_io_node(top_node,"munhits");
  _munhits_h = munhits_ptr->TableHeader();
  _munhits = munhits_ptr->TableData();
} 


void
mMuiSlowSim::digitize()
{
  // ==================================================================
  // Loop through the PISA hits, and add twopacks to the dMuiRaw table.
  // ==================================================================
  
  for (int hit=0; hit<_munhits_h.nok; hit++) {
    
    // Determine the arm and gap number from the PISA plane number.
    // PISA convention:  2001-2006 are South Arm, 1001-1006 are North.
    short PisaPlane = _munhits[hit].plane_num;   //staf table   
    /*plane_num is data memeber of munhits struct*/
    
    short iArm = -1;
    if (PisaPlane > 2000) {
      iArm = kSOUTH;
    } else {
      iArm = kNORTH;
    }
    
    short  iPlane = (PisaPlane%1000) - 1;
    //divide 1001 through 1006 each by 1000 taken remainder subtract by 1
    
    if ( (iArm<0) || (iArm>=TMuiChannelId::kArmsTotal) ||
         (iPlane<0) || (iPlane>=TMuiChannelId::kPlanesPerArm) ) {
      continue;
    } //sanity check
    
    TMuiMCHitMapO::iterator mc_hit_iter = _mc_hit_map->insert_new(iArm,iPlane);
    //mc_hit_iter is the pointer to the interface object (TMuiHitO the MC Hit Object)
    //insert_new methid is declared in file: TMuiMCHitMapO.h
    /* ONLY insert arm and plane information since there is NOT any panel or orientation
       information in PISA for the MuID */
    /* Set the location of the MCHit object */
    mc_hit_iter->get()->set_arm(iArm);    
    mc_hit_iter->get()->set_plane(iPlane);    
    
    /*the following seven set_* methods are defined in TMuiMCHitO */
    mc_hit_iter->get()->set_track_id(_munhits[hit].track_num);
    mc_hit_iter->get()->set_x(_munhits[hit].rhit[0]);    
    mc_hit_iter->get()->set_y(_munhits[hit].rhit[1]);    
    mc_hit_iter->get()->set_z(_munhits[hit].rhit[2]);    
    mc_hit_iter->get()->set_px(_munhits[hit].phit[0]);    
    mc_hit_iter->get()->set_py(_munhits[hit].phit[1]);    
    mc_hit_iter->get()->set_pz(_munhits[hit].phit[2]); 
    mc_hit_iter->get()->set_pid(_munhits[hit].trk_id);
    
      /* the get()->set_* takes pisa information and puts it into the interface object (here TMuiMCHitO) which has already had its pointer inserted into the container (here the the MCHitMap) */      
    
    // rhit and phit originate in the MuID PISA file: 
    // simulation/pisa2000/src/phnxcore/encodeRootEvntMui.cc
    Hep3Vector X(_munhits[hit].rhit[0],     //rhit is the PISA MuID HIT position
                 _munhits[hit].rhit[1],     //...position GVect
                 _munhits[hit].rhit[2]);
    // make a unite vector from phit.
    //
    float r = sqrt(_munhits[hit].phit[0]*_munhits[hit].phit[0]+_munhits[hit].phit[1]*_munhits[hit].phit[1]+_munhits[hit].phit[2]*_munhits[hit].phit[2]);
    float v1=_munhits[hit].phit[0]/r;
    float v2=_munhits[hit].phit[1]/r;
    float v3=_munhits[hit].phit[2]/r;

    /* need clarification 
    Hep3Vector V(_munhits[hit].phit[0],     //phit is the MuID PISA momentum 
                 _munhits[hit].phit[1],     //...unit vector DirVect
                 _munhits[hit].phit[2]);
    */
    Hep3Vector V(v1,v2,v3);
    std::vector<TMuiChannelId> twoPackList;
    //TMuiChannelId located at offline/packages/muigeom
    
    twoPackList = TMuiGeometry::Geom()->findTwoPacks(iArm,iPlane,X,V);      
    //TMuiGeometry located at offline/packages/muigeom
    
    //for above "findTwoPacks" accessor method declared in TMuiGeometry.hh
    // Find two-packs in the specified plane that lie along a trajectory.
    // The trajectory is given by the position GVect and the unit vector
    // DirVect in global coordinates.
    /* findTwoPacks(const short& Arm, const short& Plane, const Hep3Vector& GVect,
       const Hep3Vector &DirVect); */
    if (twoPackList.size() > 0) { /* size() is method of STL vector class */
            
      for (unsigned int i=0; i<twoPackList.size(); i++) {         
        short iPanel      = twoPackList[i].Panel();
        EOrient_t iOrient = twoPackList[i].Orient();
        short iTwoPack    = twoPackList[i].TwoPack();
        mc_hit_iter->get()->add_twopack(iOrient, iPanel, iTwoPack);
        /*  Panel(), Orient(), TwoPack() are TMuiChannelId accessor methods  */
        /*  add_two_pack is a TMuiMCHitO method */
        /* now the interface object contains the two_pack information */
      }
    }
  }
  //  cout << "mMuiReadout-I8  finished looping through PISA data" << endl;
}

void 
mMuiSlowSim::associate_mctrk() 
{
  // Loop over TMuiMCHitMap [
  //   Loop over TMutMCTrkMap [
  //     if(the trkid from mctrk map matches trkid of muimchit map) [
  //       associate muimchit and mctrk.
  //       set association flag true
  //     ]
  //   ]
  //   if(association flag is not true) [
  //       make a new mctrk and insert it in map.
  //       fill the fields of new mctrk object.
  //   ]
  // ]

  TMuiMCHitMapO::iterator muimchit_iter = _mc_hit_map->range();
  while(TMuiMCHitMapO::pointer muimchit_ptr = muimchit_iter.next()){
    TMutMCTrkMap::iterator mctrk_iter = _mc_trk_map->range();
    bool is_associated = false;
    while(TMutMCTrkMap::pointer mctrk_ptr = mctrk_iter.next()){
      if(mctrk_ptr->get()->get_track_id()!=muimchit_ptr->get()->get_track_id())
        continue;
      is_associated = true;
      PHKey::associate(muimchit_ptr,mctrk_ptr);
    }
    if(is_associated==false) fill_new_mctrk(muimchit_ptr);
  }
}

void 
mMuiSlowSim::fill_new_mctrk(TMuiMCHitMapO::pointer mc_hit_ptr)
{

  // Get the parent MC track information via a pisa call according to 
  //   the track id. 
  // Insert a new mctrk into mctrkmap.
  // Fill in the fields.
  //
  float ptot=0, ptheta=0, pphi=0, r_vertex=0, z_vertex=0, theta_vertex=0, phi_vertex=0;
  int nfile=0, error=0, itparent=0, idparent=0, idpart=0;

  int track_id = mc_hit_ptr->get()->get_track_id();

  // Insert an new TMutMCTrk into map
  //
  TMutMCTrkMap::iterator mctrk_iter = _mc_trk_map->insert_new(mc_hit_ptr->get()->get_arm());

  // Make the PISA call to get track data
  //
  dio_ptrkstack(&track_id, &nfile, &error, &ptot, &ptheta, &pphi,
                &r_vertex, &z_vertex, &theta_vertex, &phi_vertex, 
                &itparent, &idparent, &idpart); 
  
  // Fill in the fields.
  //
  mctrk_iter->get()->set_arm(mc_hit_ptr->get()->get_arm());
  mctrk_iter->get()->set_pid(idpart);
  mctrk_iter->get()->set_track_id(track_id);
  mctrk_iter->get()->set_parent_id(idparent);
  mctrk_iter->get()->set_x_orig(r_vertex*sin(theta_vertex*MUIOO::DEG_TO_RAD)*cos(phi_vertex*MUIOO::DEG_TO_RAD));
  mctrk_iter->get()->set_y_orig(r_vertex*sin(theta_vertex*MUIOO::DEG_TO_RAD)*sin(phi_vertex*MUIOO::DEG_TO_RAD));
  mctrk_iter->get()->set_z_orig(z_vertex);
  mctrk_iter->get()->set_px_orig(ptot*sin(ptheta*MUIOO::DEG_TO_RAD)*cos(pphi*MUIOO::DEG_TO_RAD));
  mctrk_iter->get()->set_py_orig(ptot*sin(ptheta*MUIOO::DEG_TO_RAD)*sin(pphi*MUIOO::DEG_TO_RAD));
  mctrk_iter->get()->set_pz_orig(ptot*cos(ptheta*MUIOO::DEG_TO_RAD));

  // UBER-kludge -- don't know how to get charge info -- set correctly
  // for muon/pion/kaon/proton at least. -- This kludge was made by sean,
  // not by chun, the smarter guy.
  //
  if(idpart==5 || idpart==8 || idpart == 11 || idpart == 14){
    mctrk_iter->get()->set_charge(1);
  } else if(idpart == 6 || idpart == 9 || idpart == 12 || idpart == 15) {
    mctrk_iter->get()->set_charge(-1);
  } else {
    MUIOO::TRACE("mMuiSlowSim got an unrecognized PID -- setting charge to 0");
    mctrk_iter->get()->set_charge(0);
  }
  // Do the association
  //
  PHKey::associate(mc_hit_ptr, mctrk_iter.current());
}

---