In This Package:

DetSimVali.cc

Go to the documentation of this file.

/*
 * DetSimVali -- DetSimValidation
 *
 * May, 2008 Zhe Wang, Kevin Zhang
 * Sep, 2008 Wei Wang, make following changes:
 * 1. Adding in the neutron capture info and relevant tree branch changes
 *    in corresponding header files;
 * 2. Adding in the center-of-gravity of deposited energy and relavant
 *    tree branch changes in cooresponding header files;;
 * 3. Restructuring the loop of the program to make it faster.
 */

#include "DetSimVali.h"

// for GenHeader
#include "Event/GenHeader.h"
#include "HepMC/GenEvent.h"

// detector element
#include "DetDesc/IGeometryInfo.h"
#include "DetDesc/DetectorElement.h"
#include "DetDesc/SolidBase.h"

// for SimHit
#include "Event/SimHeader.h"
#include "Event/SimHitHeader.h"
#include "Event/SimHitCollection.h"
#include "Event/SimPmtHit.h"
#include "Event/SimHit.h"

// for history info

#include "Event/SimParticleHistory.h"
#include "Event/SimUnobservableStatisticsHeader.h"
#include "Event/SimStatistic.h"

// for Detector ID etc.
#include "Conventions/Detectors.h"

#include "GaudiKernel/ITHistSvc.h"

#include "DetHelpers/ICoordSysSvc.h"

#include "CLHEP/Units/SystemOfUnits.h"

using namespace std;

DetSimVali::DetSimVali(const std::string& name, ISvcLocator* pSvcLocator):
  GaudiAlgorithm(name, pSvcLocator), m_tsvc(0),m_csvc(0)
{
  declareProperty("Volume",m_volume="",
                  "TDS path of volume used to generated vertices");  
}

DetSimVali::~DetSimVali()
{
}

StatusCode DetSimVali::initialize()
{
  this->GaudiAlgorithm::initialize();
  
  if ( service("THistSvc", m_tsvc).isFailure()) {
    error() << " No THistSvc available." << endreq;
    return StatusCode::FAILURE;
  } 
  
  if ( service("CoordSysSvc", m_csvc).isFailure()) {
    error() << " No CoordSysSvc available." << endreq;
    return StatusCode::FAILURE;
  }
  
  m_valiTree=new ValidationTree;
  m_valiTree->create();
  
  return StatusCode::SUCCESS;
}

StatusCode DetSimVali::execute()
{
  debug()<<"executing ... "<<endreq;
  
  m_valiTree->reset();

  
  DayaBay::GenHeader* genheader = 
    get<DayaBay::GenHeader>(DayaBay::GenHeaderLocation::Default);
    
  HepMC::GenEvent* genevt = genheader->event();
  if (!genevt) { 
    warning() << "No gen event!!" << endreq;
    //        return StatusCode::SUCCESS; 
  }

  HepMC::GenEvent::vertex_const_iterator vtxIter, vtxEnd = genevt->vertices_end();
  for(vtxIter = genevt->vertices_begin(); vtxIter != vtxEnd; vtxIter++){  
    HepMC::GenVertex* genvtx = *vtxIter;
    HepMC::GenVertex::particles_out_const_iterator pt = genvtx->particles_out_const_begin(), done = genvtx->particles_out_const_end();
  
    double energy=0;
    for (; pt != done; ++pt) {
      HepMC::GenParticle* part = *pt;
      if (part->pdg_id()==13 || part->pdg_id()==-13){
        HepMC::FourVector momentum = part->momentum();
        
        m_valiTree->MuonE=momentum.e();
        m_valiTree->MuonPx=momentum.x();
        m_valiTree->MuonPy=momentum.y();
        m_valiTree->MuonPz=momentum.z();
      }
    }
    
    //get the Detector Element:
    
    DetectorElement* detelem=0;
    
    if ("" == m_volume) {
      warning() << "No Volume property set, guessing at volume bounds." 
                << endreq;
    }
    
    if (!existDet<DetectorElement>(m_volume)) {
      fatal() << "Failed to get detector element " << m_volume 
              << ", will guess at volume bounds" << endreq;
    } else {
      detelem = getDet<DetectorElement>(m_volume);
    }
    
    if (detelem) {
      m_geo = detelem->geometry();
      
      HepMC::FourVector r4 = genvtx->position();
    
      Gaudi::XYZPoint gp(r4.x(),r4.y(),r4.z());
      Gaudi::XYZPoint lp = m_geo->toLocal(gp);
      
      m_valiTree->MuonTime=r4.t()/CLHEP::second;
      m_valiTree->MuonX=lp.x()/CLHEP::meter;
      m_valiTree->MuonY=lp.y()/CLHEP::meter;
      m_valiTree->MuonZ=lp.z()/CLHEP::meter;
    }
  }

  DayaBay::SimHeader* simheader = 
    get<DayaBay::SimHeader>(DayaBay::SimHeaderLocation::Default);
  
  
  const DayaBay::SimHitHeader*         simhitheader = simheader->hits();
  
  const DayaBay::SimHitHeader::hc_map& hcmap = simhitheader->hitCollection();
  
  DayaBay::SimHitHeader::hc_map::const_iterator it;
  
  m_valiTree->SimHitsEntries=0;
  int index(0);
  
  for (it=hcmap.begin(); it != hcmap.end(); ++it) {
    
    DayaBay::Detector det(it->first);
    debug() << "Got hit collection from " << det.detName()<<" id= "
            << det.siteDetPackedData()<<endreq;
    
    const std::vector<DayaBay::SimHit*>& hitvec=it->second->collection();
    std::vector<DayaBay::SimHit*>::const_iterator it_hvec;

    for (it_hvec=hitvec.begin(); it_hvec!=hitvec.end(); ++it_hvec) {
      
      index=m_valiTree->SimHitsEntries;
      if(index<MaxSimHitEntries){
        
        m_valiTree->hitTime[index]  = (*it_hvec)->hitTime();
        m_valiTree->hitx[index]     = (*it_hvec)->localPos().x();
        m_valiTree->hity[index]     = (*it_hvec)->localPos().y();
        m_valiTree->hitz[index]     = (*it_hvec)->localPos().z();
        m_valiTree->sensID[index]   = (*it_hvec)->sensDetId();
        m_valiTree->weight[index]   = (*it_hvec)->weight();
        
        // optical photon's ancestor's pdg
        if((*it_hvec)->ancestor().track()) {
          m_valiTree->ancestorPdg[index] = 
            (*it_hvec)->ancestor().track()->particle();
        } else {
          m_valiTree->ancestorPdg[index] = 0;
        }
        //optical photon's wavelength
        const DayaBay::SimPmtHit* pmthit = 
          static_cast<const DayaBay::SimPmtHit*>(*it_hvec);
        m_valiTree->wavelength[index]   = (*pmthit).wavelength();
        m_valiTree->polx[index]   = (*pmthit).pol().x();
        m_valiTree->poly[index]   = (*pmthit).pol().y();
        m_valiTree->polz[index]   = (*pmthit).pol().z();
        m_valiTree->px[index]   = (*pmthit).dir().x();
        m_valiTree->py[index]   = (*pmthit).dir().y();
        m_valiTree->pz[index]   = (*pmthit).dir().z();
        
        m_valiTree->SimHitsEntries = index+1;
        
      } else {
        warning()<< " Reached the max hit limit!!!!!" <<endl;
      }
    }
  }

  
  // Get Primary Track Info. <<<<<<<<<<<<<--------------
  const DayaBay::SimParticleHistory* shist = simheader->particleHistory();
  const std::list<const DayaBay::SimTrack*>& priTrks = shist->primaryTracks();
  std::list<const DayaBay::SimTrack*>::const_iterator vt,priIt, priEnd=priTrks.end();
  
  /*
  if(priTrks.size()!=0){
    
    for (vt=priTrks.begin(); vt != priEnd; ++vt) {
      cout<<"IDIDID:"<<(*vt)->parentParticle()<<endl;;
      if((*vt)->parentParticle()==0 ){
        const std::vector<DayaBay::SimVertex*> pvtx=(*vt)->vertices();
        std::vector<DayaBay::SimVertex*>::const_iterator tt, tEt=pvtx.end();
        cout<<""<<pvtx.back()->volume()->
        
        for(tt=pvtx.begin();tt!=tEt;tt++){
          cout<<"!!!!!!!"<<(*tt)->process().name()<<endl;
        }
                
                  std::string pname=(*vt)->vertices().back()->process().name();
                  if(pname=="Decay")
                  m_valiTree->StopFlag = 1;
                  if(pname=="muMinusCaptureAtRest")
                  m_valiTree->StopFlag = 2;
                  cout<<"!!!!!!!!!!"<<pname<<endl;
        
      }
    }
  }
  else
    debug()<<" Empty Track List" <<endreq;
  */
  // end of track information ------------>>>>>>>>>>>

  
  //Get UnObservable Statistics info <<<<<<<<<<<<<----------
  const DayaBay::SimUnobservableStatisticsHeader* unobSt = 
    simheader->unobservableStatistics();
  const DayaBay::SimUnobservableStatisticsHeader::stat_map& 
    statmap = unobSt->stats();
  DayaBay::SimUnobservableStatisticsHeader::stat_map::const_iterator 
    st, stdone = statmap.end();
  
  debug()<<"UnObs size "<<statmap.size()<<endreq;
  
  if(unobSt) {
    
    double ptrk1,vxtrk1,vytrk1,vztrk1;
    double ptrk2,vxtrk2,vytrk2,vztrk2;
    double QESum;
    double xQESumGdLS, yQESumGdLS, zQESumGdLS;
    double xQESumLS, yQESumLS, zQESumLS;
    double xQESumMO, yQESumMO, zQESumMO;
    
    for (st=statmap.begin(); st != stdone; ++st) {
      
      // Muon information
      
      if(st->first=="MuonTrkLengthInOws")  
        m_valiTree->MuonTrkLengthInOws = st->second.sum();
      if(st->first=="MuonTrkLengthInIws")  
        m_valiTree->MuonTrkLengthInIws = st->second.sum();
      if(st->first=="MuonTrkLengthInLS")  
        m_valiTree->MuonTrkLengthInLS = st->second.sum();
      if(st->first=="MuonTrkLengthInGdLS")  
        m_valiTree->MuonTrkLengthInGdLS = st->second.sum();

      if(st->first=="dEInn")  m_valiTree->dEInn = st->second.sum();
      if(st->first=="dEOut")  m_valiTree->dEOut = st->second.sum();
      if(st->first=="MuonStop" && st->second.sum()>0) {
        m_valiTree->StopFlag=int(st->second.sum());
        cout<< "EEEE: "<<st->second.sum();
      }
      
      // Energy information
      
      if(st->first=="EDepInGdLS") 
        m_valiTree->EDepInGdLS = st->second.sum();
      if(st->first=="EDepInLS") 
          m_valiTree->EDepInLS = st->second.sum();
      if(st->first=="EDepInIAV") 
         m_valiTree->EDepInIAV = st->second.sum();
      if(st->first=="EDepInOAV") 
         m_valiTree->EDepInOAV = st->second.sum();
      if(st->first=="EDepInOIL") 
         m_valiTree->EDepInOIL = st->second.sum();
      
      if(st->first=="QEDepInGdLS") 
        m_valiTree->QEDepInGdLS = st->second.sum();
      if(st->first=="QEDepInLS")                  
          m_valiTree->QEDepInLS = st->second.sum();
      if(st->first=="QEDepInIAV")                 
         m_valiTree->QEDepInIAV = st->second.sum();
      if(st->first=="QEDepInOAV")                 
         m_valiTree->QEDepInOAV = st->second.sum();
      if(st->first=="QEDepInOIL")                 
         m_valiTree->QEDepInOIL = st->second.sum();
      
      if(st->first == "tQESumGdLS")
        m_valiTree->tQESumGdLS = st->second.sum();
      if(st->first == "xQESumGdLS")
        xQESumGdLS = st->second.sum();
      if(st->first == "yQESumGdLS")
        yQESumGdLS = st->second.sum();
      if(st->first == "zQESumGdLS")
        zQESumGdLS = st->second.sum();
      
      if(st->first == "tQESumLS")
        m_valiTree->tQESumLS = st->second.sum();
      if(st->first == "xQESumLS")
        xQESumLS = st->second.sum();
      if(st->first == "yQESumLS")
        yQESumLS = st->second.sum();
      if(st->first == "zQESumLS")
        zQESumLS = st->second.sum();
      
      if(st->first == "tQESumMO") 
        m_valiTree->tQESumMO = st->second.sum();
      if(st->first == "xQESumMO")                                 
        xQESumMO = st->second.sum();
      if(st->first == "yQESumMO")                                 
        yQESumMO = st->second.sum();
      if(st->first == "zQESumMO")                                 
        zQESumMO = st->second.sum();
      
      // Neutron Capture

      if(st->first == "tGen") 
        m_valiTree->t_gen = st->second.sum()/CLHEP::microsecond;
      if(st->first == "xGen")                                  
        m_valiTree->x_gen = st->second.sum();
      if(st->first == "yGen")                                  
        m_valiTree->y_gen = st->second.sum();
      if(st->first == "zGen")                                  
        m_valiTree->z_gen = st->second.sum();
      
      if(st->first == "tCap") 
        m_valiTree->t_cap = st->second.sum()/CLHEP::microsecond;
      if(st->first == "xCap")                                  
        m_valiTree->x_cap = st->second.sum();
      if(st->first == "yCap")                                  
        m_valiTree->y_cap = st->second.sum();
      if(st->first == "zCap")                                  
        m_valiTree->z_cap = st->second.sum();
      
      if(st->first=="capTarget") m_valiTree->capTarget = (int)st->second.sum();
      
      
      if(st->first=="pdgId_Trk1")  
        m_valiTree->pdgId_Trk1 = (int)st->second.sum();

      if(st->first=="t_Trk1")      
        m_valiTree->t_Trk1 = st->second.sum()/CLHEP::microsecond;
      if(st->first=="x_Trk1")      
        m_valiTree->x_Trk1 = st->second.sum();
      if(st->first=="y_Trk1")      
        m_valiTree->y_Trk1 = st->second.sum();
      if(st->first=="z_Trk1")      
        m_valiTree->z_Trk1 = st->second.sum();

      if(st->first=="tEnd_Trk1") 
        m_valiTree->tEnd_Trk1 = st->second.sum()/CLHEP::microsecond;
      if(st->first=="xEnd_Trk1")                                   
        m_valiTree->xEnd_Trk1 = st->second.sum();
      if(st->first=="yEnd_Trk1")                                   
        m_valiTree->yEnd_Trk1 = st->second.sum();
      if(st->first=="zEnd_Trk1")                                   
        m_valiTree->zEnd_Trk1 = st->second.sum();
      
      if(st->first=="e_Trk1")
        m_valiTree->e_Trk1 = st->second.sum() / CLHEP::MeV;

      if(st->first=="p_Trk1")  ptrk1 = st->second.sum() / CLHEP::MeV;
      if(st->first=="vx_Trk1") vxtrk1 = st->second.sum();
      if(st->first=="vy_Trk1") vytrk1 = st->second.sum();
      if(st->first=="vz_Trk1") vztrk1 = st->second.sum();
      
      if(st->first=="ke_Trk1") 
        m_valiTree->ke_Trk1 = st->second.sum() / CLHEP::MeV;
      
      if(st->first=="TrkLength_GD_Trk1")  
        m_valiTree->TrkLength_GD_Trk1 = st->second.sum()/CLHEP::centimeter;
      if(st->first=="TrkLength_iAV_Trk1") 
        m_valiTree->TrkLength_iAV_Trk1 = st->second.sum()/CLHEP::centimeter;
      if(st->first=="TrkLength_LS_Trk1")  
        m_valiTree->TrkLength_LS_Trk1 = st->second.sum()/CLHEP::centimeter;
      if(st->first=="TrkLength_oAV_Trk1") 
        m_valiTree->TrkLength_oAV_Trk1 = st->second.sum()/CLHEP::centimeter;
      if(st->first=="TrkLength_Oil_Trk1") 
        m_valiTree->TrkLength_Oil_Trk1 = st->second.sum()/CLHEP::centimeter;


      if(st->first=="pdgId_Trk2")  
        m_valiTree->pdgId_Trk2 = (int)st->second.sum();
      if(st->first=="t_Trk2")      
        m_valiTree->t_Trk2 = st->second.sum()/CLHEP::microsecond;
      if(st->first=="x_Trk2")      
        m_valiTree->x_Trk2 = st->second.sum();
      if(st->first=="y_Trk2")      
        m_valiTree->y_Trk2 = st->second.sum();
      if(st->first=="z_Trk2")      
        m_valiTree->z_Trk2 = st->second.sum();
      
      if(st->first=="tEnd_Trk2") 
        m_valiTree->tEnd_Trk2 = st->second.sum()/CLHEP::microsecond;
      if(st->first=="xEnd_Trk2")                                   
        m_valiTree->xEnd_Trk2 = st->second.sum();
      if(st->first=="yEnd_Trk2")                                   
        m_valiTree->yEnd_Trk2 = st->second.sum();
      if(st->first=="zEnd_Trk2")                                   
        m_valiTree->zEnd_Trk2 = st->second.sum();
      
      if(st->first=="e_Trk2") 
        m_valiTree->e_Trk2 = st->second.sum() / CLHEP::MeV;

      if(st->first=="p_Trk2")  ptrk2 = st->second.sum() / CLHEP::MeV;
      if(st->first=="vx_Trk2") vxtrk2 = st->second.sum();
      if(st->first=="vy_Trk2") vytrk2 = st->second.sum();
      if(st->first=="vz_Trk2") vztrk2 = st->second.sum();
                                                            
      if(st->first=="ke_Trk2")                              
        m_valiTree->ke_Trk2 = st->second.sum() / CLHEP::MeV;
      
      if(st->first=="TrkLength_GD_Trk2")  
        m_valiTree->TrkLength_GD_Trk2 = st->second.sum()/CLHEP::centimeter;
      if(st->first=="TrkLength_iAV_Trk2") 
        m_valiTree->TrkLength_iAV_Trk2 = st->second.sum()/CLHEP::centimeter;
      if(st->first=="TrkLength_LS_Trk2")  
        m_valiTree->TrkLength_LS_Trk2 = st->second.sum()/CLHEP::centimeter;
      if(st->first=="TrkLength_oAV_Trk2") 
        m_valiTree->TrkLength_oAV_Trk2 = st->second.sum()/CLHEP::centimeter;
      if(st->first=="TrkLength_Oil_Trk2") 
        m_valiTree->TrkLength_Oil_Trk2 = st->second.sum()/CLHEP::centimeter;
      
    }
    
    Gaudi::XYZPoint gp;
    Gaudi::XYZPoint lp;
    IDetectorElement *de;
    
    if(m_valiTree->pdgId_Trk1==2112 || m_valiTree->pdgId_Trk2==2112)
    {
      gp.SetCoordinates(m_valiTree->x_gen,m_valiTree->y_gen,m_valiTree->z_gen);
      if(gp.x()!=0 || gp.y()!=0 || gp.z()!=0)
      {
        de = m_csvc->coordSysDE(gp);
        if(de!=0)
        {
          lp = de->geometry()->toLocal(gp);
          m_valiTree->x_gen = lp.x()/CLHEP::centimeter;
          m_valiTree->y_gen = lp.y()/CLHEP::centimeter;
          m_valiTree->z_gen = lp.z()/CLHEP::centimeter;
        }
        else
        {
          m_valiTree->x_gen = 0;
          m_valiTree->y_gen = 0;
          m_valiTree->z_gen = 0;
        }
      }
      
      gp.SetCoordinates(m_valiTree->x_cap,m_valiTree->y_cap,m_valiTree->z_cap);
      if(gp.x()!=0 || gp.y()!=0 || gp.z()!=0)
      {
        de = m_csvc->coordSysDE(gp);
        if(de!=0)
        {
          lp = de->geometry()->toLocal(gp);
          m_valiTree->x_cap = lp.x()/CLHEP::centimeter;
          m_valiTree->y_cap = lp.y()/CLHEP::centimeter;
          m_valiTree->z_cap = lp.z()/CLHEP::centimeter;
        }
        else
        {
          m_valiTree->x_cap = 0;
          m_valiTree->y_cap = 0;
          m_valiTree->z_cap = 0;
        }
      }
    }
    
    gp.SetCoordinates(m_valiTree->x_Trk1,m_valiTree->y_Trk1,m_valiTree->z_Trk1);
    if(gp.x()!=0 || gp.y()!=0 || gp.z()!=0)
    {
      de = m_csvc->coordSysDE(gp);
      if(de!=0)
      {
        lp = de->geometry()->toLocal(gp);
        m_valiTree->x_Trk1 = lp.x()/CLHEP::centimeter;
        m_valiTree->y_Trk1 = lp.y()/CLHEP::centimeter;
        m_valiTree->z_Trk1 = lp.z()/CLHEP::centimeter;
      }
      else
      {
        m_valiTree->x_Trk1 = 0;
        m_valiTree->y_Trk1 = 0;
        m_valiTree->z_Trk1 = 0;
      }
    }
    
    gp.SetCoordinates(m_valiTree->xEnd_Trk1,m_valiTree->yEnd_Trk1,m_valiTree->zEnd_Trk1);
    if(gp.x()!=0 || gp.y()!=0 || gp.z()!=0)
    {
      de = m_csvc->coordSysDE(gp);
      if(de!=0)
      {
        lp = de->geometry()->toLocal(gp);
        m_valiTree->xEnd_Trk1 = lp.x()/CLHEP::centimeter;
        m_valiTree->yEnd_Trk1 = lp.y()/CLHEP::centimeter;
        m_valiTree->zEnd_Trk1 = lp.z()/CLHEP::centimeter;
      }
      else
      {
        m_valiTree->xEnd_Trk1 = 0;
        m_valiTree->yEnd_Trk1 = 0;
        m_valiTree->zEnd_Trk1 = 0;
      }
    }
    
    gp.SetCoordinates(m_valiTree->x_Trk2,m_valiTree->y_Trk2,m_valiTree->z_Trk2);
    if(gp.x()!=0 || gp.y()!=0 || gp.z()!=0)
    {
      de = m_csvc->coordSysDE(gp);
      if(de!=0)
      {
        lp = de->geometry()->toLocal(gp);
        m_valiTree->x_Trk2 = lp.x()/CLHEP::centimeter;
        m_valiTree->y_Trk2 = lp.y()/CLHEP::centimeter;
        m_valiTree->z_Trk2 = lp.z()/CLHEP::centimeter;
      }
      else
      {
        m_valiTree->x_Trk2 = 0;
        m_valiTree->y_Trk2 = 0;
        m_valiTree->z_Trk2 = 0;
      }
    }
    
    gp.SetCoordinates(m_valiTree->xEnd_Trk2,m_valiTree->yEnd_Trk2,m_valiTree->zEnd_Trk2);
    if(gp.x()!=0 || gp.y()!=0 || gp.z()!=0)
    {
      de = m_csvc->coordSysDE(gp);
      if(de!=0)
      {
        lp = de->geometry()->toLocal(gp);
        m_valiTree->xEnd_Trk2 = lp.x()/CLHEP::centimeter;
        m_valiTree->yEnd_Trk2 = lp.y()/CLHEP::centimeter;
        m_valiTree->zEnd_Trk2 = lp.z()/CLHEP::centimeter;
      }
      else
      {
        m_valiTree->xEnd_Trk2 = 0;
        m_valiTree->yEnd_Trk2 = 0;
        m_valiTree->zEnd_Trk2 = 0;
      }
    }
    
    m_valiTree->px_Trk1 = ptrk1 * vxtrk1;
    m_valiTree->py_Trk1 = ptrk1 * vytrk1;
    m_valiTree->pz_Trk1 = ptrk1 * vztrk1;
    
    m_valiTree->px_Trk2 = ptrk2 * vxtrk2;
    m_valiTree->py_Trk2 = ptrk2 * vytrk2;
    m_valiTree->pz_Trk2 = ptrk2 * vztrk2;
    
    if (m_valiTree->QEDepInGdLS > 0.) {
      m_valiTree->tQESumGdLS /= (m_valiTree->QEDepInGdLS * CLHEP::microsecond);
      Gaudi::XYZPoint grQESumGdLS(xQESumGdLS / m_valiTree->QEDepInGdLS, 
                                  yQESumGdLS / m_valiTree->QEDepInGdLS, 
                                  zQESumGdLS / m_valiTree->QEDepInGdLS);
      Gaudi::XYZPoint lrQESumGdLS = m_geo->toLocal(grQESumGdLS);
    
      m_valiTree->xQESumGdLS = lrQESumGdLS.x()/CLHEP::centimeter;
      m_valiTree->yQESumGdLS = lrQESumGdLS.y()/CLHEP::centimeter;
      m_valiTree->zQESumGdLS = lrQESumGdLS.z()/CLHEP::centimeter;
    }

    if (m_valiTree->QEDepInLS > 0.) {
      m_valiTree->tQESumLS /= (m_valiTree->QEDepInLS * CLHEP::microsecond);
      Gaudi::XYZPoint grQESumLS(xQESumLS / m_valiTree->QEDepInGdLS, 
                                yQESumLS / m_valiTree->QEDepInGdLS, 
                                zQESumLS / m_valiTree->QEDepInGdLS);
      
      Gaudi::XYZPoint lrQESumLS = m_geo->toLocal(grQESumLS);

      m_valiTree->xQESumLS = lrQESumLS.x()/CLHEP::centimeter;
      m_valiTree->yQESumLS = lrQESumLS.y()/CLHEP::centimeter;
      m_valiTree->zQESumLS = lrQESumLS.z()/CLHEP::centimeter;
    }

    QESum = m_valiTree->QEDepInGdLS + m_valiTree->QEDepInLS;
    if (QESum > 0.) {
      m_valiTree->tCtrQE = 
        (m_valiTree->tQESumGdLS + m_valiTree->tQESumLS) / QESum
        / CLHEP::microsecond;
      Gaudi::XYZPoint grCtrQE( (xQESumGdLS + xQESumLS) / QESum, 
                               (yQESumGdLS + yQESumLS) / QESum,
                               (zQESumGdLS + zQESumLS) / QESum );
      Gaudi::XYZPoint lrCtrQE = m_geo->toLocal(grCtrQE);
      m_valiTree->xCtrQE = lrCtrQE.x()/CLHEP::centimeter;
      m_valiTree->yCtrQE = lrCtrQE.y()/CLHEP::centimeter;
      m_valiTree->zCtrQE = lrCtrQE.z()/CLHEP::centimeter;
    }
    
    debug() << "Capture target: " << m_valiTree->capTarget << endreq;
    
  } else {

    m_valiTree->MuonTrkLengthInOws = 0;
    m_valiTree->MuonTrkLengthInIws = 0;
    m_valiTree->MuonTrkLengthInLS = 0;
    m_valiTree->MuonTrkLengthInGdLS = 0;
    m_valiTree->StopFlag=0;
    m_valiTree->dEInn = 0;
    m_valiTree->dEOut =0;

  }
  
  const DayaBay::SimParticleHistory* history = simheader->particleHistory();
  
  if(history) {
    const std::list<DayaBay::SimTrack*>& TrackList = history->tracks();
    std::list<DayaBay::SimTrack*>::const_iterator cit;

    debug() << "Track list size: " << TrackList.size() << endreq;

    double t0, x0, y0, z0;
    double t1, x1, y1, z1;
    
    /* the center of AD1 in the global frame */
    /* Hard-coding is bad. But for now, it is only for test purpose.
       the information is available in the *_gen and *_cap branches.
       --- Wei */
    const Double_t AD1ctr_x = -1496.055;
    const Double_t AD1ctr_y = -80452.055;
    const Double_t AD1ctr_z = -710.0;
    
    for(cit=TrackList.begin(); cit!=TrackList.end(); ++cit) {

      debug()<<"pdg: "<<(*cit)->particle()<<endreq;
      debug()<<"parent:"<<(*cit)->parentParticle()<<endreq;
      
      const std::vector<DayaBay::SimVertex*>& trkVtxList = (*cit)->vertices();
      std::vector<DayaBay::SimVertex*>::const_iterator vtxi;
      
      debug() << "# of vertices: " << trkVtxList.size() << endreq;
      
      int trkVtxSize = trkVtxList.size();
      int count = 0;
      
      for (vtxi=trkVtxList.begin(); vtxi != trkVtxList.end(); vtxi++) {
        
        Gaudi::XYZPoint gp = (*vtxi)->position();
        Gaudi::XYZPoint lp = m_geo->toLocal(gp);
        
        if (count == 0) {
          
          m_valiTree->tBuff1 = (*vtxi)->time()/CLHEP::microsecond;
          
          m_valiTree->xBuff1 = lp.x()/CLHEP::centimeter;
          m_valiTree->yBuff1 = lp.y()/CLHEP::centimeter;
          m_valiTree->zBuff1 = lp.z()/CLHEP::centimeter;
          
          m_valiTree->xBuff2 = 
            (*vtxi)->position().x()/CLHEP::centimeter - AD1ctr_x;
          m_valiTree->yBuff2 = 
            (*vtxi)->position().y()/CLHEP::centimeter - AD1ctr_y;
          m_valiTree->zBuff2 = 
            (*vtxi)->position().z()/CLHEP::centimeter - AD1ctr_z;
          
        } 
        
        if (count == trkVtxSize-1) {
          t1 = (*vtxi)->time()/CLHEP::microsecond;
          x1 = (*vtxi)->position().x()/CLHEP::centimeter - AD1ctr_x;
          y1 = (*vtxi)->position().y()/CLHEP::centimeter - AD1ctr_y;
          z1 = (*vtxi)->position().z()/CLHEP::centimeter - AD1ctr_z;
        }
        
        debug() << "Neutron verties of track: "
                << (*vtxi)->process() << endreq;
        debug() << "Time: "
                << (*vtxi)->time()/CLHEP::microsecond << endreq;
        count++;
      }
      
    }
    
    const std::list<DayaBay::SimVertex*>& VertexList = history->vertices();
    std::list<DayaBay::SimVertex*>::const_iterator ci;
    
    int count=0;
    debug()<<"vertex list size: "<<VertexList.size()<<endreq;

    for(ci=VertexList.begin(); ci!=VertexList.end(); ++ci) {
      count++;
      if(count<30) {
        debug()<<"process: "<<(*ci)->process()<<endreq;
        debug()<<"energy: "<<(*ci)->totalEnergy()<<endreq;      
        debug()<< "secondaries: " << (*ci)->secondaries().size() << endreq;
        debug()<<"mass: "<<(*ci)->mass()<<endreq;
        debug()<<"totalEnergy: "<<(*ci)->totalEnergy()<<endreq;
        debug()<<"kineticEnergy: "<<(*ci)->kineticEnergy()<<endreq;
      }
    }

  }
  
  m_valiTree->fill();  
  
  return StatusCode::SUCCESS;
}

StatusCode DetSimVali::finalize()
{
  m_valiTree->close();
  delete m_valiTree;

  return this->GaudiAlgorithm::finalize();
}

---