/search.css" rel="stylesheet" type="text/css"/> /search.js">
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

In This Package:

Public Member Functions | Private Member Functions | Private Attributes
ChkMuon Class Reference

#include <ChkMuon.h>

Collaboration diagram for ChkMuon:
Collaboration graph
[legend]

List of all members.

Public Member Functions

 ChkMuon (const std::string &name, ISvcLocator *pSvcLocator)
 Check muon sample.
virtual ~ChkMuon ()
virtual StatusCode initialize ()
virtual StatusCode execute ()
virtual StatusCode finalize ()

Private Member Functions

StatusCode processGenHeader (const GenHeader *pGenHeader)
StatusCode processSimHeader (const SimHeader *pSimHeader)
StatusCode processHistory (const DayaBay::SimParticleHistory *pHist)

Private Attributes

MuonTreem_valiTree
 create a root tree
TimeStamp m_t0
 t0 of this set of data
SimTrackm_muon
 Primary muon track.
vector< SimTrack * > m_neutrons
IDetectorElement * m_dbOws
 A few key detector elements, near site only -- db.
IDetectorElement * m_dbOil1
IDetectorElement * m_dbOil2
IDetectorElement * m_dbLso1
IDetectorElement * m_dbLso2

Detailed Description

Definition at line 27 of file ChkMuon.h.


Constructor & Destructor Documentation

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

Check muon sample.

Oct 5, 2009 created by Zhe Wang

Definition at line 30 of file ChkMuon.cc.

                                                               :
  GaudiAlgorithm(name, pSvcLocator)
{
}
ChkMuon::~ChkMuon ( ) [virtual]

Definition at line 35 of file ChkMuon.cc.

{
}

Member Function Documentation

StatusCode ChkMuon::initialize ( ) [virtual]

Get geometry information

Definition at line 39 of file ChkMuon.cc.

{
  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 ( ) [virtual]

reset tree

Definition at line 57 of file ChkMuon.cc.

{
  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 ( ) [virtual]

Definition at line 105 of file ChkMuon.cc.

{
  m_valiTree->close();
  delete m_valiTree;

  return this->GaudiAlgorithm::finalize();
}
StatusCode ChkMuon::processGenHeader ( const GenHeader pGenHeader) [private]

in CLHEP/Units/SystemOfUnits.h the default energy unit is MeV.

Definition at line 117 of file ChkMuon.cc.

{
  // 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) [private]

fill sim hit SimHitHeader

typedef std::map<short int,SimHitCollection*> hc_map

loop over every SimHitCollection

loop over every SimHit

Definition at line 194 of file ChkMuon.cc.

{
  // 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) [private]

analyze muon track 1. track length in each material

analyze each neutron track

Analyze each muon induced neutron

Get where it is generated

Definition at line 309 of file ChkMuon.cc.

{
  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;

}

Member Data Documentation

create a root tree

Definition at line 45 of file ChkMuon.h.

t0 of this set of data

Definition at line 48 of file ChkMuon.h.

Primary muon track.

Definition at line 51 of file ChkMuon.h.

vector<SimTrack*> ChkMuon::m_neutrons [private]

Definition at line 52 of file ChkMuon.h.

IDetectorElement* ChkMuon::m_dbOws [private]

A few key detector elements, near site only -- db.

Definition at line 55 of file ChkMuon.h.

IDetectorElement* ChkMuon::m_dbOil1 [private]

Definition at line 56 of file ChkMuon.h.

IDetectorElement* ChkMuon::m_dbOil2 [private]

Definition at line 57 of file ChkMuon.h.

IDetectorElement* ChkMuon::m_dbLso1 [private]

Definition at line 58 of file ChkMuon.h.

IDetectorElement* ChkMuon::m_dbLso2 [private]

Definition at line 59 of file ChkMuon.h.


The documentation for this class was generated from the following files:
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 10:07:04 for MDC09b by doxygen 1.7.4