In This Package:

ChkMuon.cc

Go to the documentation of this file.

#include "ChkMuon.h"

// for GenEvent

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

// for history info
//#include "Event/SimParticleHistoryHeader.h"
#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 <Event/HeaderObject.h>


ChkMuon::ChkMuon(const std::string& name, ISvcLocator* pSvcLocator):
  GaudiAlgorithm(name, pSvcLocator)
{
}

ChkMuon::~ChkMuon()
{
}

StatusCode ChkMuon::initialize()
{
  this->GaudiAlgorithm::initialize();
  
  m_dbOws = getDet<IDetectorElement>("/dd/Structure/Pool/db-ows");
  m_dbOil1 = getDet<IDetectorElement>("/dd/Structure/AD/db-oil1");
  m_dbOil2 = getDet<IDetectorElement>("/dd/Structure/AD/db-oil2");
  m_dbLso1 = getDet<IDetectorElement>("/dd/Structure/AD/db-lso1");
  m_dbLso2 = getDet<IDetectorElement>("/dd/Structure/AD/db-lso2");

  // init root tree
  m_valiTree=new MuonTree;
  m_valiTree->create();

  return StatusCode::SUCCESS;
}

StatusCode ChkMuon::execute()
{
  info()<<"Executing ... "<<endreq;

  StatusCode   sc;
  DataObject   *pObject;

  // get ReadoutHeader
  sc=eventSvc()->retrieveObject(SimHeader::defaultLocation(),pObject);
  if(sc.isFailure()) {
    debug()<<"No SimHeader for this event"<<endreq;
  }

  const GenHeader        *pGenHeader=0;
  const SimHeader        *pSimHeader=0;

  const IHeader* firstHeader=0;

  pSimHeader=dynamic_cast<SimHeader*>(pObject);

  // get the header object pointer of it lower stages
  // only the first one, for low event rate, it is enough
  if(pSimHeader!=0) {
    // GenHeader
    if(pSimHeader->inputHeaders().size() != 0) {
      firstHeader=*((pSimHeader->inputHeaders()).begin());
      pGenHeader= dynamic_cast<const GenHeader*>(firstHeader);
    }
  }  

  m_valiTree->reset();

  m_valiTree->evt = pSimHeader->execNumber();

  // then each header
  sc = processGenHeader(pGenHeader);
  if(sc.isFailure()) return StatusCode::FAILURE;

  sc = processSimHeader(pSimHeader);
  if(sc.isFailure()) return StatusCode::FAILURE;

  m_valiTree->fill();
  debug()<<"Leaving Execute()..."<<endreq;

  return StatusCode::SUCCESS;
}

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

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

//  --------------------------------------
//  process each kind of header
//  --------------------------------------

StatusCode ChkMuon::processGenHeader(const GenHeader *pGenHeader)
{
  // check GenHeader info if available                                                                           
  if( pGenHeader==0 ) {
    return StatusCode::SUCCESS;
  }
  info()<<"Processing GenHeader"<<endreq;

  // vertex                                                                                                                          
  const HepMC::GenEvent * pgevt = pGenHeader->event();
  const HepMC::GenVertex * pvtx = 0;
  const HepMC::GenParticle * pptl =0;

  debug()<<"vertices_size "<<pgevt->vertices_size()<<endreq;
  // There should be at least one vertex.                                                                                            
  // More than one vertex is also unlikely.                                                                                          
  if ( pgevt->vertices_size() != 1 ) {
    error()<<"More than one primary vertex, out of the scope of this analysing program."<<endreq;
    return StatusCode::FAILURE;
  }
  pvtx=*(pgevt->vertices_begin());
  debug()<<"vertex: "<<*pvtx<<endreq;

  // Only one primary particle                                                                                                       
  if ( pvtx->particles_out_size() != 1 ) {
    warning()<<"More than one primary particle, might be out of the scope of this analysing program."<<endreq;
    //return StatusCode::FAILURE;                                                                                                    
  }
  pptl=*(pvtx->particles_out_const_begin());
  debug()<<"particle: "<<*pptl<<endreq;

  // fill vertex position and time                                                                                                   
  // convert global point to local point with respect to ows
  Gaudi::XYZPoint globalP(pvtx->position().x(),
                          pvtx->position().y(),
                          pvtx->position().z());

  Gaudi::XYZPoint localP = m_dbOws->geometry()->toLocal(globalP);

  m_valiTree->trk1_x=localP.X();
  m_valiTree->trk1_y=localP.Y();
  m_valiTree->trk1_z=localP.Z();
  m_valiTree->trk1_r=sqrt((m_valiTree->trk1_x * m_valiTree->trk1_x) +
                          (m_valiTree->trk1_y * m_valiTree->trk1_y));
  // get t0 of this job                           
  static bool first=true;
  if(first) {
    first=false;
    m_t0=pGenHeader->earliest();  // get t0 of this run
    info()<<"t0 of this sample: "<<m_t0<<endreq;
  }

  // substract m_t0
  TimeStamp now=pGenHeader->earliest();
  TimeStamp dt=now-m_t0;  // get rid of m_t0                                                                                             
  m_valiTree->trk1_t=
    dt.GetSeconds()             // second                                                                                            
    +pvtx->position().t()/1e9;  // nano second -> second                                                                             

  debug()<<"now: "<<now<<"\t"<<"dt: "<<dt<<endreq;
  debug()<<"dt(s): "<<dt.GetSeconds()<<"\t"<<"t_vts(ns): "<<pvtx->position().t()<<endreq;

  // fill particle momentum                                                                                                          
  m_valiTree->trk1_px=pptl->momentum().px();
  m_valiTree->trk1_py=pptl->momentum().py();
  m_valiTree->trk1_pz=pptl->momentum().pz();
  m_valiTree->trk1_e=pptl->momentum().e();


  return StatusCode::SUCCESS;
}

//
//
//

StatusCode ChkMuon::processSimHeader(const SimHeader *pSimHeader)
{
  // check SimHeader info if available
  if( pSimHeader==0 ) {
    return StatusCode::SUCCESS;
  }

  info()<<"Processing SimHeader. analyze DayaBay site only!"<<endreq;

  const SimHitHeader* simhitheader = pSimHeader->hits();
  const SimHitHeader::hc_map& hcmap = simhitheader->hitCollection();
  SimHitHeader::hc_map::const_iterator it;
  
  debug()<<"size of hit collection "<< hcmap.size()<<endreq;
  for (it=hcmap.begin(); it != hcmap.end(); ++it) {
    
    Detector det(it->first);
    debug() << "Got hit collection from " << det.detName()<<" id= "
            << det.siteDetPackedData()<<endreq;

    const std::vector<SimHit*>& hitvec=it->second->collection();
    std::vector<SimHit*>::const_iterator it_hvec;

    if( det.detectorId() == DetectorId::kAD1 ) {
      debug() <<"AD1 has "<< hitvec.size()<<" hits"<<endreq;
      m_valiTree->nhits_AD1=hitvec.size();

    } else if ( det.detectorId() == DetectorId::kAD2 ) {
      debug() <<"AD2 has "<< hitvec.size()<<" hits"<<endreq;
      m_valiTree->nhits_AD2=hitvec.size();

    } else if ( det.detectorId() == DetectorId::kIWS ) {
      debug() <<"IWS has "<< hitvec.size()<<" hits"<<endreq;
      m_valiTree->nhits_IWS=hitvec.size();

    } else if ( det.detectorId() == DetectorId::kOWS ) {
      debug() <<"OWS has "<< hitvec.size()<<" hits"<<endreq;
      m_valiTree->nhits_OWS=hitvec.size();

    } else if ( det.detectorId() == DetectorId::kRPC ) {
      debug() <<"RPC has "<< hitvec.size()<<" hits"<<endreq;
      m_valiTree->nhits_RPC=hitvec.size();
    }

  } // loop over each hit collection


  // ------> get geant4 truth                                                                                                          
  // find the truch muon vertex in each AD                                                                                             
  const DayaBay::SimParticleHistory* pHist = pSimHeader->particleHistory();

  if(pHist) {
    processHistory(pHist);
  }

  //Get UnObservable Statistics info <<<<<<<<<<<<<----------
  const DayaBay::SimUnobservableStatisticsHeader* unobSt =
    pSimHeader->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) {

    for (st=statmap.begin(); st != stdone; ++st) {

      // Muon information
      // group 1, tot
      if(st->first=="TotDeInOil1")
        m_valiTree ->TotDeInOil1 = st->second.sum();
      if(st->first=="TotDeInOil2")
        m_valiTree ->TotDeInOil2 = st->second.sum();
      if(st->first=="TotDeInLso1")
        m_valiTree ->TotDeInLso1 = st->second.sum();
      if(st->first=="TotDeInLso2")
        m_valiTree ->TotDeInLso2 = st->second.sum();

      // group 2, muon track only
      if(st->first=="MuDeInOil1")
        m_valiTree ->MuDeInOil1 = st->second.sum();
      if(st->first=="MuDeInOil2")
        m_valiTree ->MuDeInOil2 = st->second.sum();
      if(st->first=="MuDeInLso1")
        m_valiTree ->MuDeInLso1 = st->second.sum();
      if(st->first=="MuDeInLso2")
        m_valiTree ->MuDeInLso2 = st->second.sum();

      // group 3, muon track dx
      if(st->first=="MuDxInOil1")
        m_valiTree ->MuDxInOil1 = st->second.sum();
      if(st->first=="MuDxInOil2")
        m_valiTree ->MuDxInOil2 = st->second.sum();
      if(st->first=="MuDxInLso1")
        m_valiTree ->MuDxInLso1 = st->second.sum();
      if(st->first=="MuDxInLso2")
        m_valiTree ->MuDxInLso2 = st->second.sum();

    }
  }


  return StatusCode::SUCCESS;
} // end of pSimHeader

  
StatusCode ChkMuon::processHistory(const DayaBay::SimParticleHistory* pHist)
{
  m_muon = 0;
  m_neutrons.clear();

  // tracks
  const std::list<DayaBay::SimTrack*>& trk=pHist->tracks();
  std::list<DayaBay::SimTrack*>::const_iterator tkci, tkEnd=trk.end();

  info()<<"History track list size "<<trk.size()<<endreq;
  for(tkci = trk.begin(); tkci != tkEnd; tkci++ ) {
    // Find primary muon
    if( (*tkci)->trackId() == 1 &&
        ( (*tkci)->particle() == 13 || (*tkci)->particle() == -13 ) &&
        (*tkci)->ancestorTrack().track() == 0 ) {
      m_muon = *tkci;
    }

    // Find neutron
    if( (*tkci)->particle() == 2112 ) {
      m_neutrons.push_back( *tkci );
    }
  }
  info()<<m_neutrons.size()<<" spallation neutron(s)"<<endreq;
  
  const std::list<DayaBay::SimVertex*>& vtx=pHist->vertices();
  std::list<DayaBay::SimVertex*>::const_iterator vtci, vtEnd=vtx.end();
  info()<<"History vertex list size "<<vtx.size()<<endreq;

  //  -------------------------------------------------------------------
  const vector<SimVertex*>& MuonVtx = m_muon->vertices();
  for( unsigned int ci=0; ci<MuonVtx.size(); ++ci )  {
    //    info()<<MuonVtx[ci]->position()<<endreq;
  }

  // convert global point of muon vertex to local point with respect to ows
  Gaudi::XYZPoint globalP = MuonVtx[0]->position();
  Gaudi::XYZPoint localP = m_dbOws->geometry()->toLocal(globalP);

  m_valiTree->his_t1 = MuonVtx[0]->time();
  m_valiTree->his_x1 = localP.x();
  m_valiTree->his_y1 = localP.y();
  m_valiTree->his_z1 = localP.z();


  //  -------------------------------------------------------------------
  int nMuonInducedNeutron = 0;
  int nInOws = 0;
  int nInOil1 = 0;
  int nInOil2 = 0;
  int nInLso1 = 0;
  int nInLso2 = 0;

  for( unsigned int nn=0; nn<m_neutrons.size(); ++nn )  {
    SimTrack* pNeu = m_neutrons[nn];
    const vector<SimVertex*>& NeuVtx = pNeu->vertices();
    
    info()<<"spallation neutron: "<<"track id "<<pNeu->trackId()<<" PDG id "<<pNeu->particle()
          <<" ancestor "<<pNeu->ancestorTrack().track()->trackId()<<" has "<<NeuVtx.size()<<" vertices."<<endreq;

    // The first and the last vertex for this neutron track
    SimVertex* FirstVtx = NeuVtx[0];
    SimVertex* LastVtx = NeuVtx[NeuVtx.size()-1];
    info()<<"  "<<FirstVtx->time()<<" "<<FirstVtx->position()<<" "<<FirstVtx->totalEnergy()<<" {"
          <<FirstVtx->process().type()<<" "<<FirstVtx->process().name()<<" }"<<endreq;
    info()<<"  "<<LastVtx->time()<<" "<<LastVtx->position()<<" "<<LastVtx->totalEnergy()<<" {"
          <<LastVtx->process().type()<<" "<<LastVtx->process().name()<<" }"<<endreq;

    if( pNeu->ancestorTrack().track()->trackId() != 1 )  {
      // Not produced by primary muon track, not interesting.
      info()<<"  Not induced by primary muon track"<<endreq;
      continue;
    }

    nMuonInducedNeutron += 1;
    globalP.SetZ( globalP.Z()-5000 );
    if( m_dbOws->geometry()->isInside( globalP ) )   info()<<"inside"<<endreq;
    info()<<"local point "<<m_dbOws->geometry()->toLocal( globalP )<<endreq;

    Gaudi::XYZPoint localP = m_dbOws->geometry()->toLocal( FirstVtx->position() );
    info()<<"local point of neutron vertex "<<localP<<endreq;

    if( m_dbOws->geometry()->isInside( FirstVtx->position() ) )  nInOws++;
    if( m_dbOil1->geometry()->isInside( FirstVtx->position() ) )  nInOil1++;
    if( m_dbOil2->geometry()->isInside( FirstVtx->position() ) )  nInOil2++;
    if( m_dbLso1->geometry()->isInside( FirstVtx->position() ) )  nInLso1++;
    if( m_dbLso2->geometry()->isInside( FirstVtx->position() ) )  nInLso2++;
    
    
    /*
    for( unsigned int ci=0; ci<NeuVtx.size(); ++ci )  {
      info()<<"  "<<NeuVtx[ci]->time()<<" "<<NeuVtx[ci]->position()<<" "<<NeuVtx[ci]->totalEnergy()<<" {"
            <<NeuVtx[ci]->process().type()<<" "<<NeuVtx[ci]->process().name()<<" }"<<endreq;
    }
    */
  }

  info()<<"Number of muon induced spallation neutrons: "<<nMuonInducedNeutron<<endreq;
  m_valiTree->nNeutron = nMuonInducedNeutron;
  m_valiTree->nInOws = nInOws;
  m_valiTree->nInOil1 = nInOil1;
  m_valiTree->nInOil2 = nInOil2;
  m_valiTree->nInLso1 = nInLso1;
  m_valiTree->nInLso2 = nInLso2;


  return StatusCode::SUCCESS;

}

---