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

In This Package:

Public Types | Public Member Functions | Private Member Functions | Private Attributes
McTruthFigMuon Class Reference

#include <McTruthFigMuon.h>

Collaboration diagram for McTruthFigMuon:
Collaboration graph
[legend]

List of all members.

Public Types

enum  GenTypeId_t { kUnknown = 0, kIBD = 1, kMuon = 2, kRad = 3 }

Public Member Functions

 McTruthFigMuon (const std::string &name, ISvcLocator *pSvcLocator)
virtual ~McTruthFigMuon ()
virtual StatusCode initialize ()
virtual StatusCode execute ()
virtual StatusCode finalize ()

Private Member Functions

int processGenHeader (int run, const DayaBay::Detector &detector, const DayaBay::GenHeader *pGenHdr)
int processSimHeader (int run, const DayaBay::Detector &detector, const DayaBay::SimHeader *pSimHdr)
void FillXY (const DayaBay::Detector &detector, const Gaudi::XYZPoint &point, TH2 *hist)
void FillXZ (const DayaBay::Detector &detector, const Gaudi::XYZPoint &point, TH2 *hist)
void FillYZ (const DayaBay::Detector &detector, const Gaudi::XYZPoint &point, TH2 *hist)
TH1 * getOrMakeHist (int run, const DayaBay::Detector &detector, int histogram)
std::string getPath (int run, const char *histName)

Private Attributes

IStatisticsSvcm_statsSvc
ICoordSysSvcm_coordSvc
TimeStampm_firstTriggerTime
TimeStamp m_lastTriggerTime
std::map< int, TH1 ** > m_shortCuts
std::vector< TH1 * > m_normalize

Detailed Description

Definition at line 61 of file McTruthFigMuon.h.


Member Enumeration Documentation

Enumerator:
kUnknown 
kIBD 
kMuon 
kRad 

Definition at line 66 of file McTruthFigMuon.h.

                   { 
    kUnknown = 0,
    kIBD = 1,                                                                                                                             
    kMuon = 2, 
    kRad = 3
  }; 

Constructor & Destructor Documentation

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

Definition at line 29 of file McTruthFigMuon.cc.

  : GaudiAlgorithm(name,pSvcLocator),
    m_statsSvc(0),
    m_coordSvc(0),
    m_firstTriggerTime(0)
{
}
McTruthFigMuon::~McTruthFigMuon ( ) [virtual]

Definition at line 38 of file McTruthFigMuon.cc.

{
}

Member Function Documentation

StatusCode McTruthFigMuon::initialize ( ) [virtual]

Definition at line 42 of file McTruthFigMuon.cc.

{
  // Initialize the necessary services
  StatusCode sc = this->service("StatisticsSvc",m_statsSvc,true);
  if(sc.isFailure()){
    error() << "Failed to get StatisticsSvc" << endreq;
    return sc;
  }

  sc = service("CoordSysSvc", m_coordSvc, true);
  if (sc.isFailure()) {
    error() << "Failed to get CoordSysSvc" << endreq;
    return sc;
  }
  
  return sc;
}
StatusCode McTruthFigMuon::execute ( ) [virtual]

Definition at line 60 of file McTruthFigMuon.cc.

{
  // Add the current event into histograms
  DayaBay::ReadoutHeader* readoutHeader = 
    get<DayaBay::ReadoutHeader>("/Event/Readout/ReadoutHeader");
  if(!readoutHeader){
    error() << "Failed to get readout header." << endreq;
    return StatusCode::FAILURE;
  }

  const DayaBay::DaqCrate* readout = readoutHeader->daqCrate();
  if(!readout){
    error() << "Failed to get readout from header" << endreq;
    return StatusCode::FAILURE;
  }
  
  int runNumber = readout->runNumber();
  const DayaBay::Detector& detector = readout->detector();
  
  const DayaBay::GenHeader *genHdr = 0;
  const DayaBay::SimHeader *simHdr = 0;

  std::vector<const DayaBay::IHeader*> req_headers;
  std::vector<const DayaBay::IHeader*>::const_iterator header_itr;

  const DayaBay::RegistrationSequence *regseq = get<DayaBay::RegistrationSequence>
    (DayaBay::RegistrationSequence::defaultLocation() );

  const DayaBay::RegistrationSequence::Registrations& reglist = regseq->registrations();
  DayaBay::RegistrationSequence::Registrations::const_iterator regIt, regend = reglist.end();

  for (regIt=reglist.begin(); regIt!=regend; regIt++) {
    const DataObject *objreg = (*regIt).object();
    
    if (objreg->clID() == DayaBay::SimHeader::classID()) {
      simHdr = dynamic_cast<const DayaBay::SimHeader*>(objreg);
    
      bool isMuon = false;

      req_headers = simHdr->inputHeaders();
      if (req_headers.size()==0) warning() << "No related genHeaders found." << endreq;

      for (header_itr=req_headers.begin(); header_itr!=req_headers.end(); ++header_itr) {
        genHdr = dynamic_cast<const DayaBay::GenHeader*> (*header_itr);
        
        if (processGenHeader(runNumber, detector, genHdr) == kMuon) isMuon = true;        
      }
      
      if (isMuon) processSimHeader(runNumber, detector, simHdr);
    }   
  }

  return StatusCode::SUCCESS;
}
StatusCode McTruthFigMuon::finalize ( ) [virtual]

Definition at line 115 of file McTruthFigMuon.cc.

{
  // Handle Normalized histograms
  for(unsigned int normIdx=0; normIdx<m_normalize.size(); normIdx++){
    TH1* hist = m_normalize[normIdx];
    hist->Sumw2();
    if(hist->GetNormFactor()>0){
      double scale = 1./hist->GetNormFactor();
      hist->SetNormFactor(1);
      hist->Scale( scale );
    }
    hist->SetBit(TH1::kIsAverage);
  }
  m_normalize.clear();
  // Clean-up histogram shortcuts
  std::map<int,TH1**>::iterator histIter, histEnd = m_shortCuts.end();
  for(histIter = m_shortCuts.begin(); histIter!=histEnd; histIter++){
    delete [] (histIter->second);
    histIter->second = 0;
  }
  if( m_statsSvc ) m_statsSvc->release();
  if(m_firstTriggerTime){
    delete m_firstTriggerTime;
    m_firstTriggerTime = 0;
  }
  return StatusCode::SUCCESS;
}
int McTruthFigMuon::processGenHeader ( int  run,
const DayaBay::Detector detector,
const DayaBay::GenHeader pGenHdr 
) [private]

Definition at line 143 of file McTruthFigMuon.cc.

                                                                                                                    {

  debug() << "processing GenHeader" << endreq;

  if (!pGenHdr) {
    warning() << "GenHeader is empty, return" << endreq;
    return -1;
  }

  const std::string generatorName = pGenHdr->generatorName();
  const std::string Unknown = "Unknown";
  const std::string Muon = "Muon";

  int genType = kUnknown;

  size_t found = generatorName.find(Muon);

  if (found != std::string::npos) {
    genType = kMuon;
    debug() << "found Muon events" << endreq;
  }

  if (genType!=kMuon) return genType;
  
  TH1F* mcDeltaT = dynamic_cast<TH1F*>(this->getOrMakeHist(runNumber, detector, DELTATMU));
  TH1F* mcLogDeltaT = dynamic_cast<TH1F*>(this->getOrMakeHist(runNumber, detector, LOGDELTATMU));

  // Time between triggers      

  debug() << "GenHeader TimeStamp: " << pGenHdr->timeStamp().AsString() << endreq;

  if (m_firstTriggerTime) {
    TimeStamp dtTriggerTime = pGenHdr->timeStamp();
    dtTriggerTime.Subtract(m_lastTriggerTime);
    
    mcDeltaT->Fill(dtTriggerTime.GetSeconds());
    mcLogDeltaT->Fill(TMath::Log10(dtTriggerTime.GetSeconds()));
  }

  if (!m_firstTriggerTime) {
    m_firstTriggerTime =  new TimeStamp(pGenHdr->timeStamp()); 
  }
  
  m_lastTriggerTime = pGenHdr->timeStamp(); 

  TH2F* mcVtxyVSx = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber, detector, VTXYVSXMU));
  TH2F* mcVtxzVSx = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber, detector, VTXZVSXMU));
  TH2F* mcVtxzVSy = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber, detector, VTXZVSYMU));

  TH1F* mcEMu = dynamic_cast<TH1F*>(this->getOrMakeHist(runNumber, detector, EMU));
  TH1F* mcLogEMu = dynamic_cast<TH1F*>(this->getOrMakeHist(runNumber, detector, LOGEMU));
  TH1F* mcThetaMu = dynamic_cast<TH1F*>(this->getOrMakeHist(runNumber, detector, THETAMU));
  TH1F* mcPhiMu = dynamic_cast<TH1F*>(this->getOrMakeHist(runNumber, detector, PHIMU));

  TH1F* mcClosestLMu = dynamic_cast<TH1F*>(this->getOrMakeHist(runNumber, detector, CLOSESTLMU));

  const HepMC::GenEvent* event = pGenHdr->event();
  HepMC::GenEvent::particle_const_iterator p_it;

  Gaudi::XYZPoint  starting_vertex;
  Gaudi::XYZVector momentum_vector;
   
  for ( p_it = event->particles_begin(); p_it!=event->particles_end(); ++p_it ) {
    
    int pdgid = (*p_it)->pdg_id();
    int status = (*p_it)->status();

    Gaudi::XYZVector global_vector((*p_it)->momentum().px(),
                                   (*p_it)->momentum().py(),
                                   (*p_it)->momentum().pz());

    double energy = (*p_it)->momentum().e();

    if ((pdgid==13 || pdgid==-13) && status==1) { // mu+/-
      
      Gaudi::XYZPoint global_point( (*p_it)->production_vertex()->position().x(),
                                    (*p_it)->production_vertex()->position().y(),
                                    (*p_it)->production_vertex()->position().z());
      
      IDetectorElement *de = m_coordSvc->coordSysDE(global_point, 1); // site coordinate system is used
                        
      if (de) {
        debug() << "detector element: " << de->name() << endreq;
        Gaudi::XYZPoint local_point = de->geometry()->toLocal(global_point);
        Gaudi::XYZVector local_vector = de->geometry()->toLocal(global_vector);

        starting_vertex = local_point;
        momentum_vector = local_vector;
      }
    
      debug() << "muon starting vertex: ["
              << starting_vertex.X() << ", "
              << starting_vertex.Y() << ", "
              << starting_vertex.Z() << "]" << endreq;
      
      mcVtxyVSx->Fill(starting_vertex.X()/Gaudi::Units::m, starting_vertex.Y()/Gaudi::Units::m);
      mcVtxzVSx->Fill(starting_vertex.X()/Gaudi::Units::m, starting_vertex.Z()/Gaudi::Units::m);
      mcVtxzVSy->Fill(starting_vertex.Y()/Gaudi::Units::m, starting_vertex.Z()/Gaudi::Units::m);
      
      double theta = 180 - momentum_vector.theta() * 180 / TMath::Pi();
      double phi = momentum_vector.phi() * 180 / TMath::Pi();

      mcEMu->Fill(energy / Gaudi::Units::GeV);
      mcLogEMu->Fill(TMath::Log10(energy / Gaudi::Units::GeV));

      mcThetaMu->Fill(TMath::Cos(TMath::Pi() - momentum_vector.theta()));
      mcPhiMu->Fill(phi);

      debug() << "muon momentum: ["
              << momentum_vector.X() << ", "
              << momentum_vector.Y() << ", "
              << momentum_vector.Z() << "]" 
              << " energy: ["
              << energy << ", "
              << theta  << ", "
              << phi    << "]"
              << endreq;

      IDetectorElement *ows = 0;

      if (detector.site() == Site::kDayaBay) {
        ows = getDet<IDetectorElement>("/dd/Structure/Pool/db-ows");
      } 

      if (detector.site() == Site::kLingAo) {
        ows = getDet<IDetectorElement>("/dd/Structure/Pool/la-ows");
      }

      if (detector.site() == Site::kFar) {
        ows = getDet<IDetectorElement>("/dd/Structure/Pool/far-ows");
      }
      
      if (ows) {

        Gaudi::XYZPoint local_origin(0, 0, 0);
        Gaudi::XYZPoint global_origin = ows->geometry()->toGlobal(local_origin);
        
        if (de) {
          Gaudi::XYZPoint site_origin = de->geometry()->toLocal(global_origin);

          // calculate the closest point and distance to the ows center
          Gaudi::XYZVector displacement_vector(site_origin.X()-starting_vertex.X(), 
                                               site_origin.Y()-starting_vertex.Y(), 
                                               site_origin.Z()-starting_vertex.Z());
          
          double cos_theta = (momentum_vector.X()*displacement_vector.X() + 
                              momentum_vector.Y()*displacement_vector.Y() + 
                              momentum_vector.Z()*displacement_vector.Z()) /
            TMath::Sqrt(momentum_vector.Mag2() * displacement_vector.Mag2());
          
          double sin_theta = TMath::Sqrt(1-cos_theta*cos_theta);

          double closest_dist = TMath::Sqrt(displacement_vector.Mag2()) * sin_theta;
          double longitude_dist = TMath::Sqrt(displacement_vector.Mag2()) * cos_theta;
          
          Gaudi::XYZPoint closest_point(starting_vertex.X() + momentum_vector.unit().X() * longitude_dist,
                                        starting_vertex.Y() + momentum_vector.unit().Y() * longitude_dist,
                                        starting_vertex.Z() + momentum_vector.unit().Z() * longitude_dist);

          mcClosestLMu->Fill(closest_dist / Gaudi::Units::m);

          debug() << " site_origin: [" 
                  << site_origin.X() << ", "
                  << site_origin.Y() << ", "
                  << site_origin.Z() << "] "
                  << " closest_point: ["
                  << closest_point.X() << ", "
                  << closest_point.Y() << ", "
                  << closest_point.Z() << "] "
                  << " closest_dist: " << closest_dist << endreq;
        }
      }
    }
  }

  return genType;
}
int McTruthFigMuon::processSimHeader ( int  run,
const DayaBay::Detector detector,
const DayaBay::SimHeader pSimHdr 
) [private]

Definition at line 322 of file McTruthFigMuon.cc.

                                                                                                                    {
  
  debug() << "processing SimHeader" << endreq;

  if (!pSimHdr) {
    warning() << "SimHeader is empty, return" << endreq;
    return -1;
  }

  // unobserver statistics

  double AllEDepInLS1 = 0;
  double AllEDepInGdLS1 = 0;
  double AllEDepInMO1 = 0;
    
  double AllEDepInLS2 = 0;
  double AllEDepInGdLS2 = 0;
  double AllEDepInMO2 = 0;

  double AllEDepInIWS = 0;
  double AllEDepInOWS = 0;

  double AllQEDepInLS1 = 0;
  double AllQEDepInGdLS1 = 0;

  double AllQEDepInLS2 = 0;
  double AllQEDepInGdLS2 = 0;
  
  double MuTrackInLS1 = 0;
  double MuTrackInGdLS1 = 0;
  double MuTrackInMO1 = 0;
  
  double MuTrackInLS2 = 0;
  double MuTrackInGdLS2 = 0;
  double MuTrackInMO2 = 0;

  double MuTrackInIWS = 0;
  double MuTrackInOWS = 0;

  double TotalEnergyinAD = 0;
  double TotalEnergyinWS = 0;

  double DepositeEnergyinLS = 0;
  double QuenchedEnergyinLS = 0;

  double TrackLengthinAD = 0;
  double TrackLengthinWS = 0;

  double SumHitsEnergyinAD = 0;
  double SumHitsEnergyinWS = 0;

  TH1F* mcNSpallationN = dynamic_cast<TH1F*>(this->getOrMakeHist(runNumber, detector, NSPALLATIONN));
  
  TH1F* mcTNCapMU = dynamic_cast<TH1F*>(this->getOrMakeHist(runNumber, detector, TNCAPMU));
  TH2F* mcNCapyVSx = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber, detector, YVSXNCAPMU));
  TH2F* mcNCapzVSx = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber, detector, ZVSXNCAPMU));
  TH2F* mcNCapzVSy = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber, detector, ZVSYNCAPMU));

  TH2F* mcDepEVSTrkLinAD = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber, detector, DEPOVSLTRKAD));
  TH2F* mcDepEVSTrkLinWS = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber, detector, DEPOVSLTRKWS));

  TH2F* mcQuenVSDepEinAD = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber, detector, QUENVSDEPOAD));

  TH2F* mcSumhitsVSTrkLinAD = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber, detector, HITSVSLTRKAD));
  TH2F* mcSumhitsVSTrkLinWS = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber, detector, HITSVSLTRKWS));

  const DayaBay::SimUnobservableStatisticsHeader* unobservableheader = pSimHdr->unobservableStatistics();

  if (!unobservableheader) {
    warning() << "No SimUnobservableStatisticsHeader" << endreq;
  } else {

    const DayaBay::SimUnobservableStatisticsHeader::stat_map& statmap = unobservableheader->stats();
    DayaBay::SimUnobservableStatisticsHeader::stat_map::const_iterator stIt, stend = statmap.end();
    debug() << " Unobservable Statistics: " << statmap.size() << endreq;

    for (stIt = statmap.begin(); stIt != stend; ++stIt) {

      if (stIt->first == "AllEDepInLS1")    AllEDepInLS1 = stIt->second.sum();
      if (stIt->first == "AllEDepInGdLS1")  AllEDepInGdLS1 = stIt->second.sum();
      if (stIt->first == "AllEDepInMO1")    AllEDepInMO1 = stIt->second.sum();
      if (stIt->first == "AllEDepInLS2")    AllEDepInLS2 = stIt->second.sum();
      if (stIt->first == "AllEDepInGdLS2")  AllEDepInGdLS2 = stIt->second.sum();
      if (stIt->first == "AllEDepInMO2")    AllEDepInMO2 = stIt->second.sum();
      if (stIt->first == "AllEDepInIWS")    AllEDepInIWS = stIt->second.sum();
      if (stIt->first == "AllEDepInOWS")    AllEDepInOWS = stIt->second.sum();

      if (stIt->first == "AllQEDepInLS1")    AllQEDepInLS1 = stIt->second.sum();
      if (stIt->first == "AllQEDepInGdLS1")  AllQEDepInGdLS1 = stIt->second.sum();
      if (stIt->first == "AllQEDepInLS2")    AllQEDepInLS2 = stIt->second.sum();
      if (stIt->first == "AllQEDepInGdLS2")  AllQEDepInGdLS2 = stIt->second.sum();
      
      if (stIt->first == "MuTrackInLS1")     MuTrackInLS1 = stIt->second.sum();
      if (stIt->first == "MuTrackInGdLS1")   MuTrackInGdLS1 = stIt->second.sum();
      if (stIt->first == "MuTrackInMO1")     MuTrackInMO1 = stIt->second.sum();
      if (stIt->first == "MuTrackInLS2")     MuTrackInLS2 = stIt->second.sum();
      if (stIt->first == "MuTrackInGdLS2")   MuTrackInGdLS2 = stIt->second.sum();
      if (stIt->first == "MuTrackInMO2")     MuTrackInMO2 = stIt->second.sum();
      if (stIt->first == "MuTrackInIWS")     MuTrackInIWS = stIt->second.sum();
      if (stIt->first == "MuTrackInOWS")     MuTrackInOWS = stIt->second.sum();
    }
  }
      
  TotalEnergyinAD = (AllEDepInLS1 + AllEDepInGdLS1 + AllEDepInMO1
                     + AllEDepInLS2 + AllEDepInGdLS2 + AllEDepInMO2) / Gaudi::Units::GeV;

  TotalEnergyinWS = (AllEDepInIWS + AllEDepInOWS) / Gaudi::Units::GeV;

  DepositeEnergyinLS = (AllEDepInLS1 + AllEDepInGdLS1 
                        + AllEDepInLS2 + AllEDepInGdLS2) / Gaudi::Units::GeV;
  
  QuenchedEnergyinLS = (AllQEDepInLS1 + AllQEDepInGdLS1 
                        + AllQEDepInLS2 + AllQEDepInGdLS2) / Gaudi::Units::GeV;
  
  TrackLengthinAD = (MuTrackInLS1 + MuTrackInGdLS1 + MuTrackInMO1 
                     + MuTrackInLS2 + MuTrackInGdLS2 + MuTrackInMO2) / Gaudi::Units::m;

  TrackLengthinWS = (MuTrackInIWS + MuTrackInOWS) / Gaudi::Units::m;

  debug() << " AllEDepInLS1: "  << AllEDepInLS1
          << " AllEDepInGdLS1: " << AllEDepInGdLS1
          << " AllEDepInMO1: " << AllEDepInMO1
          << " AllEDepInLS2: "  << AllEDepInLS2
          << " AllEDepInGdLS2: " << AllEDepInGdLS2
          << " AllEDepInMO2: " << AllEDepInMO2
          << " AllEDepInIWS: " << AllEDepInIWS
          << " AllEDepInOWS: " << AllEDepInOWS
          
          << " AllQEDepInLS1: " << AllQEDepInLS1
          << " AllQEDepInGdLS1: " << AllQEDepInGdLS1
          << " AllQEDepInLS2: " << AllQEDepInLS2
          << " AllQEDepInGdLS2: " << AllQEDepInGdLS2
    
          << " MuTrackInLS1: " << MuTrackInLS1
          << " MuTrackInGdLS1: " << MuTrackInGdLS1
          << " MuTrackInMO1: " << MuTrackInMO1
          << " MuTrackInLS2: " << MuTrackInLS2
          << " MuTrackInGdLS2: " << MuTrackInGdLS2
          << " MuTrackInMO2: " << MuTrackInMO2
          << " MuTrackInIWS: " << MuTrackInIWS
          << " MuTrackInOWS: " << MuTrackInOWS
          << endreq;

  int nSpallationNeutron = 0;
  
  // historian information

  const DayaBay::SimParticleHistory* history = pSimHdr->particleHistory();

  if (!history) {
    warning() << "No SimParticleHistory" << endreq;
  } else {
    
    const std::list<DayaBay::SimTrack*>& tracks = history->tracks();
    debug() << "total " << tracks.size() << " tracks in the history" << endreq;

    std::list<DayaBay::SimTrack*>::const_iterator trkIt, trkEnd = tracks.end();
    
    for (trkIt = tracks.begin(); trkIt!=trkEnd; ++trkIt) {
      debug() << " particle: " << (*trkIt)->particle() 
              << " parentId: " << (*trkIt)->parentParticle()
              << endreq;
      
      /*
        if ( (*trkIt)->particle() == 13 || (*trkIt)->particle() == -13 ) { // muon 

        const std::vector<DayaBay::SimVertex*>& vertex = (*trkIt)->vertices();
        
        if (vertex.size()>0) {
          
          Gaudi::XYZPoint global_point( vertex[0]->position().x(),
                                        vertex[0]->position().y(),
                                        vertex[0]->position().z());

          IDetectorElement *de = m_coordSvc->coordSysDE(global_point, 1); // site coordinate system is used
          
          if (de) {
            debug() << "detector element: " << de->name() << endreq;
            Gaudi::XYZPoint local_point = de->geometry()->toLocal(global_point);
            
            debug() << "muon starting point: [" 
                    << local_point.X() << ", "
                    << local_point.Y() << ", "
                    << local_point.Z() << "] "
                    << endreq;
          }
        }

        std::vector<DayaBay::SimVertex*>::const_iterator vtxIt, vtxend = vertex.end();
        for (vtxIt = vertex.begin(); vtxIt!=vtxend; ++vtxIt) {
          debug() << " muon status: " << (*vtxIt)->process().name() 
                  << " kE: " << (*vtxIt)->kineticEnergy()
                  << " time: " << (*vtxIt)->time()
                  << " position: [" 
                  << (*vtxIt)->position().X() << ", "
                  << (*vtxIt)->position().Y() << ", "
                  << (*vtxIt)->position().Z() << "] "
                  << endreq;
        }
      }
      */

      if ( (*trkIt)->particle() == 2112 ) {

        const std::vector<DayaBay::SimVertex*>& vertex = (*trkIt)->vertices();
        std::vector<DayaBay::SimVertex*>::const_iterator vtxIt, vtxend = vertex.end();

        for (vtxIt = vertex.begin(); vtxIt!=vtxend; ++vtxIt) {

          std::string processname = (*vtxIt)->process().name();
          DayaBay::SimProcess::Type type = (*vtxIt)->process().type();
            
          if (type == DayaBay::SimProcess::kHadronic) {
              
            double capture_time = (*vtxIt)->time();
              
            if (processname == "nCapture")  {
            
              Gaudi::XYZPoint global_point((*vtxIt)->position().X(),
                                           (*vtxIt)->position().Y(),
                                           (*vtxIt)->position().Z());
            
              IDetectorElement *de = m_coordSvc->coordSysDE(global_point, 1); 
              
              if (de) {
                //debug() << "detector element: " << de->name() << endreq;
                Gaudi::XYZPoint end_point = de->geometry()->toLocal(global_point);

                mcTNCapMU->Fill(capture_time/Gaudi::Units::microsecond);
                mcNCapyVSx->Fill(end_point.X()/Gaudi::Units::m, end_point.Y()/Gaudi::Units::m);
                mcNCapzVSx->Fill(end_point.X()/Gaudi::Units::m, end_point.Z()/Gaudi::Units::m);
                mcNCapzVSy->Fill(end_point.Y()/Gaudi::Units::m, end_point.Z()/Gaudi::Units::m);
              }
              nSpallationNeutron++;
            }
          }

          debug() << " neutron status: " << processname
                  << " type: " << type
                  << " kE: " << (*vtxIt)->kineticEnergy()
                  << endreq;
        }
      }
    }
  }

  mcNSpallationN->Fill(nSpallationNeutron);
  
  // sim hit information

  const DayaBay::SimHitHeader* simhitheader = pSimHdr->hits();

  if (!simhitheader) {
    warning() << "No SimHitHeader" << endreq;
    
  } else {

    const DayaBay::SimHitHeader::hc_map& hcmap = simhitheader->hitCollection();
    DayaBay::SimHitHeader::hc_map::const_iterator hc_it, hc_end = hcmap.end();
    debug() << "size of hit collection " << hcmap.size() << endreq;
    
    int sumPEinAD = 0;
    int sumPEinWS = 0;
    
    for (hc_it=hcmap.begin(); hc_it!=hc_end; ++hc_it) { // loop over hit collections

      DayaBay::Detector detector(hc_it->first);

      const std::vector<DayaBay::SimHit*>& hitvec = hc_it->second->collection();
      std::vector<DayaBay::SimHit*>::const_iterator hvec_it, hvec_end = hitvec.end();
      debug() << detector.detName() << " has " << hitvec.size() << " hits." << endreq;

      if (detector.isAD()) {

        int ring, column = 0;

        for (hvec_it=hitvec.begin(); hvec_it!=hvec_end; ++hvec_it) { 
        
          DayaBay::AdPmtSensor pmtId = DayaBay::AdPmtSensor((*hvec_it)->sensDetId());

          ring = pmtId.ring()-1;
          column = pmtId.column()-1;
          if (ring<0||ring>7||column<0||column>23) continue;
        
          //double hittime = (*hvec_it)->hitTime();
          sumPEinAD++;
        }
      }

      if (detector.isWaterShield()) {
        sumPEinWS += hitvec.size();
      }
    }
    
    SumHitsEnergyinAD = sumPEinAD;
    SumHitsEnergyinWS = sumPEinWS;
  }

  if (TrackLengthinAD>0) {
    mcDepEVSTrkLinAD->Fill(TrackLengthinAD, TotalEnergyinAD/TrackLengthinAD);
    mcSumhitsVSTrkLinAD->Fill(TrackLengthinAD, SumHitsEnergyinAD/TrackLengthinAD);
  }

  if (TrackLengthinWS>0) {
    mcDepEVSTrkLinWS->Fill(TrackLengthinWS, TotalEnergyinWS/TrackLengthinWS);
    mcSumhitsVSTrkLinWS->Fill(TrackLengthinWS, SumHitsEnergyinWS/TrackLengthinWS);
  }

  if (QuenchedEnergyinLS>0) {
    mcQuenVSDepEinAD->Fill(QuenchedEnergyinLS, DepositeEnergyinLS/QuenchedEnergyinLS);
  }
    
  return 0;
}
void McTruthFigMuon::FillXY ( const DayaBay::Detector detector,
const Gaudi::XYZPoint &  point,
TH2 *  hist 
) [private]

Definition at line 935 of file McTruthFigMuon.cc.

{
  
  if (detector.site() == Site::kDayaBay || detector.site() == Site::kLingAo) { 
   
    Gaudi::XYZPoint shifted_point( point.X(), 
                                   point.Y() + 3000,
                                   point.Z());
    if (hist)  hist->Fill(shifted_point.X()/Gaudi::Units::m, shifted_point.Y()/Gaudi::Units::m);
  }

  if (detector.site() == Site::kFar) {
    if (hist) hist->Fill(point.X()/Gaudi::Units::m, point.Y()/Gaudi::Units::m);
  }
  
  return;
}
void McTruthFigMuon::FillXZ ( const DayaBay::Detector detector,
const Gaudi::XYZPoint &  point,
TH2 *  hist 
) [private]

Definition at line 953 of file McTruthFigMuon.cc.

{

  if (detector.site() == Site::kDayaBay || detector.site() == Site::kLingAo) {
    
    Gaudi::XYZPoint shifted_point(point.X(),
                                  point.Y(),
                                  point.Z() + 8000);
    if (hist) hist->Fill(shifted_point.X()/Gaudi::Units::m, shifted_point.Z()/Gaudi::Units::m);
  }

  if (detector.site() == Site::kFar) {
    
    Gaudi::XYZPoint shifted_point;

    if (point.Y()>0) { // shift AD1 & AD2 into upper Z
      shifted_point.SetXYZ(point.X(), 
                           point.Y(),
                           point.Z() + 8000);
    } else {
      shifted_point.SetXYZ(point.X(),
                           point.Y(),
                           point.Z() + 2000);
    }

    if (hist) hist->Fill(shifted_point.X()/Gaudi::Units::m, shifted_point.Z()/Gaudi::Units::m);
  }
  
  return;
}
void McTruthFigMuon::FillYZ ( const DayaBay::Detector detector,
const Gaudi::XYZPoint &  point,
TH2 *  hist 
) [private]

Definition at line 984 of file McTruthFigMuon.cc.

{
  

  if (detector.site() == Site::kDayaBay || detector.site() == Site::kLingAo) {

    Gaudi::XYZPoint shifted_point;

    if (point.X()<0) {  // AD1
      shifted_point.SetXYZ(point.X(),
                           point.Y() - 3000,
                           point.Z() + 8000);
    } else { // AD2
      shifted_point.SetXYZ(point.X(),
                           point.Y() + 3000,
                           point.Z() + 8000);
    }

    if (hist) hist->Fill(shifted_point.Y()/Gaudi::Units::m, shifted_point.Z()/Gaudi::Units::m);
  }

  if (detector.site() == Site::kFar) {
    
    Gaudi::XYZPoint shifted_point;

    if (point.X()<0 && point.Y()>0) { // AD1
      shifted_point.SetXYZ(point.X(),
                           point.Y(),
                           point.Z() + 8000);
    }

    if (point.X()>0 && point.Y()>0) { // AD2
      shifted_point.SetXYZ(point.X(),
                           point.Y() + 6000,
                           point.Z() + 8000);
    }

    if (point.X()<0 && point.Y()<0) { // AD3
      shifted_point.SetXYZ(point.X(),
                           point.Y() - 6000,
                           point.Z() + 2000);
    }

    if (point.X()>0 && point.Y()<0) { // AD4
      shifted_point.SetXYZ(point.X(),
                           point.Y(),
                           point.Z() + 2000);
    }

    if (hist) hist->Fill(shifted_point.Y()/Gaudi::Units::m, shifted_point.Z()/Gaudi::Units::m);
  }

  return;
}
TH1 * McTruthFigMuon::getOrMakeHist ( int  run,
const DayaBay::Detector detector,
int  histogram 
) [private]

Definition at line 639 of file McTruthFigMuon.cc.

{
  int siteInt = 0;
  switch(detector.site()){
  case Site::kUnknown:
    siteInt = 0; break;
  case Site::kDayaBay:
    siteInt = 1; break;
  case Site::kLingAo:
    siteInt = 2; break;
  case Site::kFar:
    siteInt = 3; break;
  case Site::kSAB:
    siteInt = 4; break;
  default:
    error() << "Unknown site: " << detector.detName() << endreq;
    return 0;
  }

  //int detectorInt = int(detector.detectorId());
  std::map<int,TH1**>::iterator histIter = m_shortCuts.find(run);
  TH1** histograms = 0;
  if( histIter == m_shortCuts.end() ){
    // Initialize histogram shortcuts
    histograms = new TH1*[MAXMCTRUTHMUONHISTS];
    for(int i=0; i<MAXMCTRUTHMUONHISTS; i++) histograms[i] = 0;
    m_shortCuts[run] = histograms;
  }else{
    // Found run
    histograms = histIter->second;
  }
  
  TH1* hist = 
    histograms[siteInt * NMCTRUTHMUONHISTOGRAMS
               + histogram];
  if(!hist){
    // Make this histogram
    std::string histName;
    
    if(detector.detectorId()){
      // Make detector histogram
      switch(histogram){

      case EMU:

        histName = "mc_E_MUON";
        hist = new TH1F(histName.c_str(),"MC-Mu Total Energy #mu",
                        500, 0, 500);
        hist->GetXaxis()->SetTitle("Energy (GeV)");
        hist->GetYaxis()->SetTitle("Events");
        break;

      case LOGEMU:

        histName = "mc_LogE_MUON";
        hist = new TH1F(histName.c_str(),"MC-Mu Total Energy #mu",
                        100, -2, 6);
        hist->GetXaxis()->SetTitle("log_{10}^{}(Energy/GeV)");
        hist->GetYaxis()->SetTitle("Events");
        break;

      case THETAMU:

        histName = "mc_Theta_MUON";
        hist = new TH1F(histName.c_str(), "MC-Mu Zenith Angle #mu",
                        100, 0, 1);
        hist->GetXaxis()->SetTitle("cos(#pi-#theta)");
        hist->GetYaxis()->SetTitle("Events");
        break;

      case PHIMU:
        
        histName = "mc_Phi_MUON";
        hist = new TH1F(histName.c_str(), "MC-Mu azimuth Angle #mu",
                        370, -185, 185);
        hist->GetXaxis()->SetTitle("#phi");
        hist->GetYaxis()->SetTitle("Events");
        break;

      case VTXYVSXMU:

        histName = "mc_VtxyVSx_MUON";
        
        if (siteInt!=3) {
          hist = new TH2F(histName.c_str(), "MC-Mu #mu Start Vtx y .vs. x",
                          240, -12, 12, 180, -9, 9);
        } else {
          hist = new TH2F(histName.c_str(), "MC-Mu #mu Start Vtx y .vs. x",
                          240, -12, 12, 240, -12, 12);
        }
        
        hist->GetXaxis()->SetTitle("x (m)");
        hist->GetYaxis()->SetTitle("y (m)");
        break;

      case VTXZVSXMU:
        
        histName = "mc_VtxzVSx_MUON";
        hist = new TH2F(histName.c_str(), "MC-Mu #mu Start Vtx z .vs. x",
                        240, -12, 12, 150, -13, 2);
        hist->GetXaxis()->SetTitle("x (m)");
        hist->GetYaxis()->SetTitle("z (m)");
        break;
        
      case VTXZVSYMU:

        histName = "mc_VtxzVSy_MUON";

        if (siteInt!=3) {
          hist = new TH2F(histName.c_str(), "MC-Mu #mu Start Vtx z .vs. y",
                          180, -9, 9, 150, -13, 2);
        } else {
          hist = new TH2F(histName.c_str(), "MC-Mu #mu Start Vtx z .vs. y",
                          240, -12, 12, 150, -13, 2);
        }

        hist->GetXaxis()->SetTitle("y (m)");
        hist->GetYaxis()->SetTitle("z (m)");
        break;

      case DELTATMU:

        histName = "mc_DeltaT_MUON";
        hist = new TH1F(histName.c_str(),
                        "MC-Mu Time Between Events",
                        100, 0, 0.015);

        hist->GetXaxis()->SetTitle("#Deltat [s]");
        hist->GetYaxis()->SetTitle("Events");
        break;

      case LOGDELTATMU:

        histName = "mc_LogDeltaT_MUON";
        hist = new TH1F(histName.c_str(),
                        "MC-Mu Time Between Events",
                        100, -10, 0);

        hist->GetXaxis()->SetTitle("log_{10}^{}(#Deltat/1s)");
        hist->GetYaxis()->SetTitle("Events");
        break;

      case CLOSESTLMU:

        histName = "mc_ClosestL_MUON";
        hist = new TH1F(histName.c_str(),
                        "MC-Mu Closest Distance to the OWS Center",
                        75, 0, 15);
        hist->GetXaxis()->SetTitle("Distance (m)");
        hist->GetYaxis()->SetTitle("Event");
        break;
        
      case NSPALLATIONN:

        histName = "mc_NSpallationN_MUON";
        hist = new TH1F(histName.c_str(),
                        "MC-MC No. of Spallation Neutrons per muon",
                        30, 0, 30);
        hist->GetXaxis()->SetTitle("NO. of Neutron/#mu");
        hist->GetYaxis()->SetTitle("Event");
        break;

      case TNCAPMU:
        
        histName = "mc_TNCapture_MUON";
        hist = new TH1F(histName.c_str(), 
                        "MC-Mu Neutron Capture Time",
                        100, 0, 2000);

        hist->GetXaxis()->SetTitle("#DeltaT [#mus]");
        hist->GetYaxis()->SetTitle("Events");
        break;

      case YVSXNCAPMU:
        
        histName = "mc_yVSxNCapture_MUON";

        if (siteInt!=3) {
          hist = new TH2F(histName.c_str(), "MC-Mu Neutron Capture Vtx y .vs. x",
                          240, -12, 12, 180, -9, 9);
        } else {
          hist = new TH2F(histName.c_str(), "MC-Mu Neutron Capture Vtx y .vs. x",
                          240, -12, 12, 240, -12, 12);
        }

        hist->GetXaxis()->SetTitle("x (m)");
        hist->GetYaxis()->SetTitle("y (m)");
        break;

      case ZVSYNCAPMU:

        histName = "mc_zVSyNCapture_MUON";

        if (siteInt!=3) {
          hist = new TH2F(histName.c_str(), "MC-Mu Neutron Capture Vtx z .vs. y",
                          180, -9, 9, 150, -13, 2);
        } else {
          hist = new TH2F(histName.c_str(), "MC-Mu Neutron Capture Vtx z .vs. y",
                          240, -12, 12, 150, -13, 2);
        }

        hist->GetXaxis()->SetTitle("y (m)");
        hist->GetYaxis()->SetTitle("z (m)");
        break;

      case ZVSXNCAPMU:

        histName = "mc_zVSxNCapture_MUON";

        hist = new TH2F(histName.c_str(), "MC-Mu Neutron Capture Vtx z .vs. x",
                        240, -12, 12, 150, -13, 2);
        hist->GetXaxis()->SetTitle("x (m)");
        hist->GetYaxis()->SetTitle("z (m)");
        break;

      case DEPOVSLTRKAD:
        
        histName = "mc_EdepositVSTrkLeninAD_MUON";
        hist = new TH2F(histName.c_str(),
                        "MC-Mu Deposited E / Track Length in AD",
                        60, 0,  15,
                        100, 0,  5);
        hist->GetXaxis()->SetTitle("Track Length (m)");
        hist->GetYaxis()->SetTitle("Deposit E / Track Length (GeV/m)");
        break;

      case DEPOVSLTRKWS:
        
        histName = "mc_EdepositVSTrkLeninWS_MUON";
        hist = new TH2F(histName.c_str(),
                        "MC-Mu Deposited E / Track Length in WS",
                        60,  0, 15,
                        100, 0,  5);
        hist->GetXaxis()->SetTitle("Track Length (m)");
        hist->GetYaxis()->SetTitle("Deposit E / Track Length (GeV/m)");
        break;

      case QUENVSDEPOAD:

        histName = "mc_QuenchVSdepositEinLS_MUON";
        hist = new TH2F(histName.c_str(),
                        "MC-Mu Deposited E / Quenched E in LS",
                        100, 0, 5,
                        100, 0, 10);
        hist->GetXaxis()->SetTitle("Quenched Energy (GeV)");
        hist->GetYaxis()->SetTitle("Deposited E / Quenched E");
        break;

      case HITSVSLTRKAD:
        
        histName = "mc_simhitsVSTrkLeninAD_MUON";
        hist = new TH2F(histName.c_str(),
                        "MC-Mu Simhits PE / Trk Length in AD",
                        60, 0,  15, 
                        100, 0, 300);
        hist->GetXaxis()->SetTitle("Track Length (m)");
        hist->GetYaxis()->SetTitle("SimHits PEs / Trk Length (PEs/m)");
        break;

      case HITSVSLTRKWS:

        histName = "mc_simhitsVSTrkLeninWS_MUON";
        hist = new TH2F(histName.c_str(),
                        "MC-Mu Simhits PE / Track Length in WS",
                        60, 0,  15, 
                        100, 0, 300);
        hist->GetXaxis()->SetTitle("Track Length (m))");
        hist->GetYaxis()->SetTitle("SimHits PEs / Trk Length (PEs/m)");
        break;

      default:
        error() << "Unknown Histogram: " << histogram << endreq;
        return 0;
      }
    }
    
    m_statsSvc->put( this->getPath(run, histName.c_str()), 
                     hist );
    histograms[siteInt * NMCTRUTHMUONHISTOGRAMS
               + histogram] = hist;
  }

  return hist;
}
std::string McTruthFigMuon::getPath ( int  run,
const char *  histName 
) [private]

Definition at line 925 of file McTruthFigMuon.cc.

{
  // Construct histogram path in statistics service
  std::stringstream path;
  path << "/file1/diagnostics/run_" << std::setfill('0') << std::setw(7) 
       << run;
  path << "/mctruth/source_muon/" << histName;
  return path.str();
}

Member Data Documentation

Definition at line 94 of file McTruthFigMuon.h.

Definition at line 97 of file McTruthFigMuon.h.

Definition at line 100 of file McTruthFigMuon.h.

Definition at line 101 of file McTruthFigMuon.h.

std::map<int,TH1**> McTruthFigMuon::m_shortCuts [private]

Definition at line 102 of file McTruthFigMuon.h.

std::vector<TH1*> McTruthFigMuon::m_normalize [private]

Definition at line 103 of file McTruthFigMuon.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:06:10 for McTruthFigs by doxygen 1.7.4