/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
McTruthFigIBD Class Reference

#include <McTruthFigIBD.h>

Collaboration diagram for McTruthFigIBD:
Collaboration graph
[legend]

List of all members.

Public Types

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

Public Member Functions

 McTruthFigIBD (const std::string &name, ISvcLocator *pSvcLocator)
virtual ~McTruthFigIBD ()
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 79 of file McTruthFigIBD.h.


Member Enumeration Documentation

Enumerator:
kUnknown 
kIBD 
kMuon 
kRad 

Definition at line 84 of file McTruthFigIBD.h.

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

Constructor & Destructor Documentation

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

Definition at line 29 of file McTruthFigIBD.cc.

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

Definition at line 38 of file McTruthFigIBD.cc.

{
}

Member Function Documentation

StatusCode McTruthFigIBD::initialize ( ) [virtual]

Definition at line 42 of file McTruthFigIBD.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 McTruthFigIBD::execute ( ) [virtual]

Definition at line 60 of file McTruthFigIBD.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 isIBD = 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) == kIBD) isIBD = true;        
      }
      
      if (isIBD) processSimHeader(runNumber, detector, simHdr);
    }   
  }

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

Definition at line 115 of file McTruthFigIBD.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 McTruthFigIBD::processGenHeader ( int  run,
const DayaBay::Detector detector,
const DayaBay::GenHeader pGenHdr 
) [private]

Definition at line 143 of file McTruthFigIBD.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 IBD = "IBD";

  int genType = kUnknown;

  size_t found = generatorName.find(IBD);

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

  if (genType!=kIBD) return genType;
  
  TH1F* mcDeltaT = dynamic_cast<TH1F*>(this->getOrMakeHist(runNumber, detector, DELTAT));
  TH1F* mcLogDeltaT = dynamic_cast<TH1F*>(this->getOrMakeHist(runNumber, detector, LOGDELTAT));

  // Time between triggers      

  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(); 

  TH1F* mcThetaNu = dynamic_cast<TH1F*>(this->getOrMakeHist(runNumber, detector, THETANU));
  TH1F* mcPhiNu = dynamic_cast<TH1F*>(this->getOrMakeHist(runNumber, detector, PHINU));

  TH1F* mcEanue = dynamic_cast<TH1F*>(this->getOrMakeHist(runNumber, detector, EANUE));
  TH1F* mcKEp = dynamic_cast<TH1F*>(this->getOrMakeHist(runNumber, detector, KEP));
  TH1F* mcThetaP = dynamic_cast<TH1F*>(this->getOrMakeHist(runNumber, detector, THETAP));
  TH1F* mcKEn = dynamic_cast<TH1F*>(this->getOrMakeHist(runNumber, detector, KEN));
  TH1F* mcThetaN = dynamic_cast<TH1F*>(this->getOrMakeHist(runNumber, detector, THETAN));

  TH2F* mcVtxyVSx = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber, detector, VTXYVSX));
  TH2F* mcVtxzVSx = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber, detector, VTXZVSX));
  TH2F* mcVtxzVSy = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber, detector, VTXZVSY));

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

  Gaudi::XYZPoint  interaction_vertex;
  Gaudi::XYZVector antinue_vector;
  Gaudi::XYZVector positron_vector;
  Gaudi::XYZVector neutron_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();
    double kinermatic = energy - (*p_it)->momentum().m();

    if (pdgid==-12 && status==2) { // anti-nue
      
      Gaudi::XYZPoint global_point( (*p_it)->end_vertex()->position().x(),
                                    (*p_it)->end_vertex()->position().y(),
                                    (*p_it)->end_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);
        interaction_vertex = local_point;

        // get the anti-nue direction in the site coordinates
        Gaudi::XYZVector local_vector = de->geometry()->toLocal(global_vector);
        mcThetaNu->Fill(local_vector.Theta() * 180 / TMath::Pi());
        mcPhiNu->Fill(local_vector.Phi() * 180 / TMath::Pi());
      }
    
      debug() << "ane end vertex: ["
             << interaction_vertex.X() << ", "
             << interaction_vertex.Y() << ", "
             << interaction_vertex.Z() << "]" << endreq;

      this->FillXY(detector, interaction_vertex, mcVtxyVSx);
      this->FillXZ(detector, interaction_vertex, mcVtxzVSx);
      this->FillYZ(detector, interaction_vertex, mcVtxzVSy);

      antinue_vector = global_vector;
           
      debug() << "ane momentum: ["
              << antinue_vector.X() << ", "
              << antinue_vector.Y() << ", "
              << antinue_vector.Z() << "]" << endreq;
      
      mcEanue->Fill(energy);
    }
      
    if (pdgid == 2112 && status==1) { // neutron
        
      neutron_vector = global_vector;
        
      debug() << "neutron momentum: ["
             << neutron_vector.X() << ", "
             << neutron_vector.Y() << ", "
             << neutron_vector.Z() << "]" << endreq;

      mcKEn->Fill(kinermatic);

      if (antinue_vector.R()>0 && neutron_vector.R()>0) {
        double costheta = antinue_vector.Dot(neutron_vector) / antinue_vector.R() / neutron_vector.R();
        debug() << "costheta : " << costheta << endreq;
        mcThetaN->Fill(costheta);
      }
    }

    if (pdgid == -11 && status==1) { // positron
        
      positron_vector = global_vector;

      debug() << "positron momentum: ["
             << positron_vector.X() << ", "
             << positron_vector.Y() << ", "
             << positron_vector.Z() << "]" << endreq;

      mcKEp->Fill(kinermatic);
      
      if (antinue_vector.R()>0 && positron_vector.R()>0) {
        double costheta = antinue_vector.Dot(positron_vector) / antinue_vector.R() / positron_vector.R();
        debug() << "costheta: " << costheta << endreq;
        mcThetaP->Fill(costheta);
      }
    }
  }

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

Definition at line 298 of file McTruthFigIBD.cc.

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

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

  TH1F* mcNCAPZ = dynamic_cast<TH1F*>(this->getOrMakeHist(runNumber, detector, NCAPZ));

  TH1F* mcTNCapGd = dynamic_cast<TH1F*>(this->getOrMakeHist(runNumber, detector, TNCAPGD));
  TH2F* mcyVSxNCapGd = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber, detector, YVSXNCAPGD)); 
  TH2F* mczVSxNCapGd = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber, detector, ZVSXNCAPGD)); 
  TH2F* mczVSyNCapGd = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber, detector, ZVSYNCAPGD)); 
  
  TH1F* mcTNCapH = dynamic_cast<TH1F*>(this->getOrMakeHist(runNumber, detector, TNCAPH));
  TH2F* mcyVSxNCapH = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber, detector, YVSXNCAPH));
  TH2F* mczVSxNCapH = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber, detector, ZVSXNCAPH));
  TH2F* mczVSyNCapH = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber, detector, ZVSYNCAPH));

  TH1F* mcTNCapC = dynamic_cast<TH1F*>(this->getOrMakeHist(runNumber, detector, TNCAPC));
  TH2F* mcyVSxNCapC = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber, detector, YVSXNCAPC));
  TH2F* mczVSxNCapC = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber, detector, ZVSXNCAPC));
  TH2F* mczVSyNCapC = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber, detector, ZVSYNCAPC));

  TH1F* mcTNCapO = dynamic_cast<TH1F*>(this->getOrMakeHist(runNumber, detector, TNCAPO));
  TH2F* mcyVSxNCapO = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber, detector, YVSXNCAPO));
  TH2F* mczVSxNCapO = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber, detector, ZVSXNCAPO));
  TH2F* mczVSyNCapO = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber, detector, ZVSYNCAPO));

  TH1F* mcTNAbs = dynamic_cast<TH1F*>(this->getOrMakeHist(runNumber, detector, TNABS));
  TH2F* mcyVSxNAbs = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber, detector, YVSXNABS));
  TH2F* mczVSxNAbs = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber, detector, ZVSXNABS));
  TH2F* mczVSyNAbs = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber, detector, ZVSYNABS));

  TH2F* mcSumhitsVSqen = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber, detector, HITSVSQUEN));
  TH2F* mcQuenVSDepo = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber, detector, QUENVSDEPO));

  // unobserver statistics

  double EDepInGdLS = 0;
  double EDepInLS = 0;

  double QEDepInGdLS = 0;
  double QEDepInLS = 0;

  int capTarget = 0;
     
  double tGen = 0;
  double tCap = 0;
  double xCap = 0;
  double yCap = 0;
  double zCap = 0;
  
  double t_Trk2 = 0;
  double tEnd_Trk2 = 0;
  double xEnd_Trk2 = 0;
  double yEnd_Trk2 = 0;
  double zEnd_Trk2 = 0;

  double DepositedEnergy = 0;
  double QuenchedEnergy = 0;
  double SumHitsEnergy = 0;
 
  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 == "EDepInGdLS") EDepInGdLS = stIt->second.sum();
      if (stIt->first == "EDepInLS")   EDepInLS = stIt->second.sum();
 
      if (stIt->first == "QEDepInGdLS") QEDepInGdLS = stIt->second.sum();
      if (stIt->first == "QEDepInLS")   QEDepInLS = stIt->second.sum();

      if (stIt->first == "capTarget")  capTarget = stIt->second.sum();

      if (stIt->first == "tGen")  tGen = stIt->second.sum();
      if (stIt->first == "tCap")  tCap = stIt->second.sum();
      if (stIt->first == "xCap")  xCap = stIt->second.sum();
      if (stIt->first == "yCap")  yCap = stIt->second.sum();
      if (stIt->first == "zCap")  zCap = stIt->second.sum();

      if (stIt->first == "t_Trk2") t_Trk2 = stIt->second.sum();
      if (stIt->first == "tEnd_Trk2") tEnd_Trk2 = stIt->second.sum();
      if (stIt->first == "xEnd_Trk2") xEnd_Trk2 = stIt->second.sum();
      if (stIt->first == "yEnd_Trk2") yEnd_Trk2 = stIt->second.sum();
      if (stIt->first == "zEnd_Trk2") zEnd_Trk2 = stIt->second.sum();
    }
  }

  debug() << " capTarget: "  << capTarget 
          << " EDepInGdLS: " << EDepInGdLS
          << " EDepInLS:   " << EDepInLS
          << " QEDepInGdLS: " << QEDepInGdLS
          << " QEDepInLS: " << QEDepInLS
          << " tGen: " << tGen
          << " tCap: " << tCap
          << " xCap: " << xCap
          << " yCap: " << yCap
          << " zCap: " << zCap
          << " t_Trk2: " << t_Trk2
          << " tEnd_Trk2: " << tEnd_Trk2
          << " xEnd_Trk2: " << xEnd_Trk2
          << " yEnd_Trk2: " << yEnd_Trk2
          << " zEnd_Trk2: " << zEnd_Trk2 
          << endreq;
  
  Gaudi::XYZPoint capture_point;
  Gaudi::XYZPoint absorption_point;

  if (capTarget>0) { // neutron capture
    Gaudi::XYZPoint global_point( xCap, yCap, zCap);
    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);
      capture_point = local_point;
    }
  } else if (capTarget == 0){ // neutron absorption
    Gaudi::XYZPoint global_point( xEnd_Trk2, yEnd_Trk2, zEnd_Trk2);
    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);
      absorption_point = local_point;
    }
  } else {
    warning() << "Unknown Target : " << capTarget << endreq;
  }

  mcNCAPZ->Fill(capTarget);

  if (capTarget == 64) {  // Gd
    mcTNCapGd->Fill((tCap - tGen)/Gaudi::Units::microsecond);
    this->FillXY(detector, capture_point, mcyVSxNCapGd);
    this->FillXZ(detector, capture_point, mczVSxNCapGd);
    this->FillYZ(detector, capture_point, mczVSyNCapGd);
  }

  if (capTarget == 1) { // H
    mcTNCapH->Fill((tCap - tGen)/Gaudi::Units::microsecond);
    this->FillXY(detector, capture_point, mcyVSxNCapH);
    this->FillXZ(detector, capture_point, mczVSxNCapH);
    this->FillYZ(detector, capture_point, mczVSyNCapH);
  }

  if (capTarget == 6) { // C
    mcTNCapC->Fill((tCap - tGen)/Gaudi::Units::microsecond);
    this->FillXY(detector, capture_point, mcyVSxNCapC);
    this->FillXZ(detector, capture_point, mczVSxNCapC);
    this->FillYZ(detector, capture_point, mczVSyNCapC);
  }
  
  if (capTarget>0 && capTarget!=1  && capTarget!=6 && capTarget!=64) {
    mcTNCapO->Fill((tCap - tGen)/Gaudi::Units::microsecond);
    this->FillXY(detector, capture_point, mcyVSxNCapO);
    this->FillXZ(detector, capture_point, mczVSxNCapO);
    this->FillYZ(detector, capture_point, mczVSyNCapO);
  }

  if (capTarget==0) {
    mcTNAbs->Fill((tEnd_Trk2-t_Trk2)/Gaudi::Units::microsecond);
    this->FillXY(detector, absorption_point, mcyVSxNAbs);
    this->FillXZ(detector, absorption_point, mczVSxNAbs);
    this->FillYZ(detector, absorption_point, mczVSyNAbs);
  }

  DepositedEnergy = EDepInGdLS + EDepInLS;
  QuenchedEnergy = QEDepInGdLS + QEDepInLS;
  
  // 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;
    
    for (hc_it=hcmap.begin(); hc_it!=hc_end; ++hc_it) { // loop over hit collections

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

      int sumPE = 0;
      int sumPE_positron = 0;
      int sumPE_neutron = 0;

      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()) continue;

      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();
        sumPE++;

        if (!(*hvec_it)->ancestor().isBad()) {
          const DayaBay::SimTrack* trk = (*hvec_it)->ancestor().track();
                  
          if (trk->particle() == -11)  sumPE_positron++;
          if (trk->particle() == 2112) sumPE_neutron++;
        }
      }
      
      debug() << " totalPE: " << sumPE
              << " sumPE_positron: " << sumPE_positron
              << " sumPE_neutron: " << sumPE_neutron
              << endreq;

      SumHitsEnergy += sumPE;
    }
  }

  debug() << " Reco E: " << SumHitsEnergy 
          << " quen E: " << QuenchedEnergy
          << " Dept E: " << DepositedEnergy
          << endreq; 

  if (QuenchedEnergy>0) {
    mcSumhitsVSqen->Fill(QuenchedEnergy, SumHitsEnergy/QuenchedEnergy); 
  }

  if (DepositedEnergy>0) {
    mcQuenVSDepo->Fill(QuenchedEnergy, QuenchedEnergy/DepositedEnergy);
  }

  return 0;
}
void McTruthFigIBD::FillXY ( const DayaBay::Detector detector,
const Gaudi::XYZPoint &  point,
TH2 *  hist 
) [private]

Definition at line 1103 of file McTruthFigIBD.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 McTruthFigIBD::FillXZ ( const DayaBay::Detector detector,
const Gaudi::XYZPoint &  point,
TH2 *  hist 
) [private]

Definition at line 1121 of file McTruthFigIBD.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 McTruthFigIBD::FillYZ ( const DayaBay::Detector detector,
const Gaudi::XYZPoint &  point,
TH2 *  hist 
) [private]

Definition at line 1152 of file McTruthFigIBD.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 * McTruthFigIBD::getOrMakeHist ( int  run,
const DayaBay::Detector detector,
int  histogram 
) [private]

Definition at line 555 of file McTruthFigIBD.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*[MAXMCTRUTHIBDHISTS];
    for(int i=0; i<MAXMCTRUTHIBDHISTS; i++) histograms[i] = 0;
    m_shortCuts[run] = histograms;
  }else{
    // Found run
    histograms = histIter->second;
  }
  
  TH1* hist = 
    histograms[siteInt * NMCTRUTHIBDHISTOGRAMS
               + histogram];
  if(!hist){
    // Make this histogram
    std::string histName;
    if(detector.detectorId()){
      // Make detector histogram
      switch(histogram){

      case THETANU:

        histName = "mc_ThetaNu_IBD";
        hist = new TH1F(histName.c_str(),"MC-IBD Incoming #bar{#nu}_{e} zenith angle",
                        190, -5, 185);
        hist->GetXaxis()->SetTitle("#theta");
        hist->GetYaxis()->SetTitle("Events");
        break;

      case PHINU:

        histName = "mc_PhiNu_IBD";
        hist = new TH1F(histName.c_str(),"MC-IBD Incoming #bar{#nu}_{e} azimuth angle",
                        100, -185, 185);
        hist->GetXaxis()->SetTitle("#phi");
        hist->GetYaxis()->SetTitle("Events");
        break;

      case EANUE:

        histName = "mc_Eanue_IBD";
        hist = new TH1F(histName.c_str(),"MC-IBD Total Energy #bar{#nu}_{e}",
                        60,0,15);
        hist->GetXaxis()->SetTitle("Energy (MeV)");
        hist->GetYaxis()->SetTitle("Events");
        break;

      case KEP:

        histName = "mc_KEP_IBD";
        hist = new TH1F(histName.c_str(), "MC-IBD Kinetic Energy of e^{+}",
                        40, 0, 10);
        hist->GetXaxis()->SetTitle("Kinetic Energy (MeV)");
        hist->GetYaxis()->SetTitle("Events");
        break;

      case THETAP:

        histName = "mc_ThetaP_IBD";
        hist = new TH1F(histName.c_str(), "MC-IBD e^{+} cos#theta w.r.t #bar{#nu}_{e}",
                        44, -1.1, 1.1);
        hist->GetXaxis()->SetTitle("cos#theta");
        hist->GetYaxis()->SetTitle("Events");
        break;

      case KEN:
        
        histName = "mc_KEN_IBD";
        hist = new TH1F(histName.c_str(), "MC-IBD Kinetic Energy of Neutron",
                        30, 0, 0.15);
        hist->GetXaxis()->SetTitle("Kinetic Energy (MeV)");
        hist->GetYaxis()->SetTitle("Events");
        break;

      case THETAN:

        histName = "mc_ThetaN_IBD";
        hist = new TH1F(histName.c_str(), "MC-IBD N cos#theta w.r.t #bar{#nu}_{e}",
                        44, -1.1, 1.1);
        hist->GetXaxis()->SetTitle("Cos#theta");
        hist->GetYaxis()->SetTitle("Events");
        break;

      case VTXYVSX:

        histName = "mc_VtxyVSx_IBD";
        
        if (siteInt!=3) {
          hist = new TH2F(histName.c_str(), "MC-IBD Event Vertex y .vs. x",
                          120, -6, 6, 60, 0, 6);
          hist->GetXaxis()->SetTitle("x (m)");
          hist->GetYaxis()->SetTitle("y [AD1,AD2: y+3](m)");

        } else {
          hist = new TH2F(histName.c_str(), "MC-IBD Event Vertex y .vs. x",
                          120, -6, 6, 120, -6, 6);

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

      case VTXZVSX:
        
        histName = "mc_VtxzVSx_IBD";
        if (siteInt!=3) {
          hist = new TH2F(histName.c_str(), "MC-IBD Event Vertex z .vs. x",
                          120, -6, 6, 70, 0, 7);
          hist->GetXaxis()->SetTitle("x (m)");
          hist->GetYaxis()->SetTitle("z [AD1, AD2: z+8](m)");

        } else {
          hist = new TH2F(histName.c_str(), "MC-IBD Event Vertex z .vs. x",
                          120, -6, 6, 140, -7, 7);
          hist->GetXaxis()->SetTitle("x (m)");
          hist->GetYaxis()->SetTitle("z [AD1, AD2: z+8; AD3,AD4: z+2](m)");
        }
        break;
        
      case VTXZVSY:

        histName = "mc_VtxzVSy_IBD";

        if (siteInt!=3) {
          hist = new TH2F(histName.c_str(), "MC-IBD Event Vertex z .vs. y",
                          120, -6, 6, 70, 0, 7);
          hist->GetXaxis()->SetTitle("y [AD1: y-3; AD2: y+3](m)");
          hist->GetYaxis()->SetTitle("z [AD1,AD2: z+8](m)");
        } else {
          hist = new TH2F(histName.c_str(), "MC-IBD Event Vertex z .vs. y",
                          120, -6, 6, 140, -7, 7);
          hist->GetXaxis()->SetTitle("y [AD2: y+6; AD3: y-6](m)");
          hist->GetYaxis()->SetTitle("z [AD1,AD2: z+8, AD3,AD4: z+2](m)");
        }

        break;

      case DELTAT:

        histName = "mc_DeltaT_IBD";
        hist = new TH1F(histName.c_str(),
                        "MC-IBD Time Between Events",
                        100, 0, 100);

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

      case LOGDELTAT:

        histName = "mc_LogDeltaT_IBD";
        hist = new TH1F(histName.c_str(),
                        "MC-IBD Time Between Events",
                        100, -6, 4);

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

      case NCAPZ:

        histName = "mc_NCapZ_IBD";
        hist = new TH1F(histName.c_str(),
                        "MC-IBD Neutron Capture Target Z", 80, -1, 79);
        hist->GetXaxis()->SetTitle("Target Z");
        hist->GetYaxis()->SetTitle("Events");
        break;

      case TNCAPGD:
        
        histName = "mc_TNCaptureinGd_IBD";
        hist = new TH1F(histName.c_str(), 
                        "MC-IBD Neutron Capture Time in Gd",
                        100, 0, 2000);

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

      case YVSXNCAPGD:
        
        histName = "mc_yVSxNCaptureinGd_IBD";

        if (siteInt!=3) {
          hist = new TH2F(histName.c_str(), "MC-IBD Neutron Capture Vtx y .vs. x in Gd",
                          120, -6, 6, 60, 0, 6);
          hist->GetXaxis()->SetTitle("x (m)");
          hist->GetYaxis()->SetTitle("y [AD1,AD2: y+3](m)");
        } else {
          hist = new TH2F(histName.c_str(), "MC-IBD Neutron Capture Vtx y .vs. x in Gd",
                          120, -6, 6, 120, -6, 6);
          hist->GetXaxis()->SetTitle("x (m)");
          hist->GetYaxis()->SetTitle("y (m)");
        }
        break;

      case ZVSYNCAPGD:

        histName = "mc_zVSyNCaptureinGd_IBD";

        if (siteInt!=3) {
          hist = new TH2F(histName.c_str(), "MC-IBD Neutron Capture Vtx z .vs. y in Gd",
                          120, -6, 6, 70, 0, 7);
          hist->GetXaxis()->SetTitle("y [AD1: y-3; AD2: y+3](m)");
          hist->GetYaxis()->SetTitle("z [AD1,AD2: z+8](m)");
        } else {
          hist = new TH2F(histName.c_str(), "MC-IBD Neutron Caputre Vtx z .vs. y in Gd",
                          120, -6, 6, 140, -7, 7);
          hist->GetXaxis()->SetTitle("y [AD2: y+6; AD3: y-6](m)");
          hist->GetYaxis()->SetTitle("z [AD1,AD2: z+8, AD3,AD4: z+2](m)");
        }
        break;

      case ZVSXNCAPGD:

        histName = "mc_zVSxNCaptureinGd_IBD";
        if (siteInt!=3) {
          hist = new TH2F(histName.c_str(), "MC-IBD Neutron Capture Vtx z .vs. x in Gd",
                          120, -6, 6, 70, 0, 7);
          hist->GetXaxis()->SetTitle("x (m)");
          hist->GetYaxis()->SetTitle("z [AD1, AD2: z+8](m)");
        } else {
          hist = new TH2F(histName.c_str(), "MC-IBD Neutron Capture Vtx z .vs. x in Gd",
                          120, -6, 6, 140, -7, 7);
          hist->GetXaxis()->SetTitle("x (m)");
          hist->GetYaxis()->SetTitle("z [AD1, AD2: z+8; AD3,AD4: z+2](m)");
        }
        break;
        
      case TNCAPH:
        
        histName = "mc_TNCaptureinH_IBD";
        hist = new TH1F(histName.c_str(), 
                        "MC-IBD Neutron Capture Time in H",
                        100, 0, 2000);

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

      case YVSXNCAPH:
       
        histName = "mc_yVSxNCaptureinH_IBD";
        if (siteInt!=3) {
          hist = new TH2F(histName.c_str(), "MC-IBD Neutron Capture Vtx y .vs. x in H",
                          120, -6, 6, 60, 0, 6);
          hist->GetXaxis()->SetTitle("x (m)");
          hist->GetYaxis()->SetTitle("y [AD1,AD2: y+3](m)");
        } else {
          hist = new TH2F(histName.c_str(), "MC-IBD Neutron Capture Vtx y .vs. x in H",
                          120, -6, 6, 120, -6, 6);
          hist->GetXaxis()->SetTitle("x (m)");
          hist->GetYaxis()->SetTitle("y (m)");
        }
        break;

      case ZVSYNCAPH:

        histName = "mc_zVSyNCaptureinH_IBD";
        
        if (siteInt!=3) {
          hist = new TH2F(histName.c_str(), "MC-IBD Neutron Capture Vtx z .vs. y in H",
                          120, -6, 6, 70, 0, 7);
          hist->GetXaxis()->SetTitle("y [AD1: y-3; AD2: y+3](m)");
          hist->GetYaxis()->SetTitle("z [AD1,AD2: z+8](m)");
        } else {
          hist = new TH2F(histName.c_str(), "MC-IBD Neutron Caputre Vtx z .vs. y in H",
                          120, -6, 6, 140, -7, 7);
          hist->GetXaxis()->SetTitle("y [AD2: y+6; AD3: y-6](m)");
          hist->GetYaxis()->SetTitle("z [AD1,AD2: z+8, AD3,AD4: z+2](m)");
        }
        break;

      case ZVSXNCAPH:

        histName = "mc_zVSxNCaptureinH_IBD";
        
        if (siteInt!=3) {
          hist = new TH2F(histName.c_str(), "MC-IBD Neutron Capture Vtx z .vs. x in H",
                          120, -6, 6, 70, 0, 7);
          hist->GetXaxis()->SetTitle("x (m)");
          hist->GetYaxis()->SetTitle("z [AD1, AD2: z+8](m)");
        } else {
          hist = new TH2F(histName.c_str(), "MC-IBD Neutron Capture Vtx z .vs. x in H",
                          120, -6, 6, 140, -7, 7);
          hist->GetXaxis()->SetTitle("x (m)");
          hist->GetYaxis()->SetTitle("z [AD1, AD2: z+8; AD3,AD4: z+2](m)");
        }
        break;

      case TNCAPC:

        histName = "mc_TNCaptureinC_IBD";
        hist = new TH1F(histName.c_str(),
                        "MC-IBD Neutron Capture Time in C",
                        100, 0, 2000);

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

        histName = "mc_yVSxNCaptureinC_IBD";
        if (siteInt!=3) {
          hist = new TH2F(histName.c_str(), "MC-IBD Neutron Capture Vtx y .vs. x in C",
                          120, -6, 6, 60, 0, 6);
          hist->GetXaxis()->SetTitle("x (m)");
          hist->GetYaxis()->SetTitle("y [AD1,AD2: y+3](m)");
        } else {
          hist = new TH2F(histName.c_str(), "MC-IBD Neutron Capture Vtx y .vs. x in C",
                          120, -6, 6, 120, -6, 6);
          hist->GetXaxis()->SetTitle("x (m)");
          hist->GetYaxis()->SetTitle("y (m)");
        }
        break;

      case ZVSYNCAPC:

        histName = "mc_zVSyNCaptureinC_IBD";

        if (siteInt!=3) {
          hist = new TH2F(histName.c_str(), "MC-IBD Neutron Capture Vtx z .vs. y in C",
                          120, -6, 6, 70, 0, 7);
          hist->GetXaxis()->SetTitle("y [AD1: y-3; AD2: y+3](m)");
          hist->GetYaxis()->SetTitle("z [AD1,AD2: z+8](m)");
        } else {
          hist = new TH2F(histName.c_str(), "MC-IBD Neutron Caputre Vtx z .vs. y in C",
                          120, -6, 6, 140, -7, 7);
          hist->GetXaxis()->SetTitle("y [AD2: y+6; AD3: y-6](m)");
          hist->GetYaxis()->SetTitle("z [AD1,AD2: z+8, AD3,AD4: z+2](m)");
        }
        break;

      case ZVSXNCAPC:

        histName = "mc_zVSxNCaptureinC_IBD";

        if (siteInt!=3) {
          hist = new TH2F(histName.c_str(), "MC-IBd Neutron Capture Vtx z .vs. x in C",
                          120, -6, 6, 70, 0, 7);
          hist->GetXaxis()->SetTitle("x (m)");
          hist->GetYaxis()->SetTitle("z [AD1, AD2: z+8](m)");
        } else {
          hist = new TH2F(histName.c_str(), "MC-IBD Neutron Capture Vtx z .vs. x in C",
                          120, -6, 6, 140, -7, 7);
          hist->GetXaxis()->SetTitle("x (m)");
          hist->GetYaxis()->SetTitle("z [AD1, AD2: z+8; AD3,AD4: z+2](m)");
        }
        break;

      case TNCAPO:
        
        histName = "mc_TNCaptureinOthers_IBD";
        hist = new TH1F(histName.c_str(),  "MC-IBD Neutron Capture Time in Others",
                        100, 0, 2000);

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

      case YVSXNCAPO:
        
        histName = "mc_yVSxNCaptureinOthers_IBD";

        if (siteInt!=3) {
          hist = new TH2F(histName.c_str(), "MC-IBD Neutron Capture Vtx y .vs. x in Others",
                          120, -6, 6, 60, 0, 6);
          hist->GetXaxis()->SetTitle("x (m)");
          hist->GetYaxis()->SetTitle("y [AD1,AD2: y+3](m)");
        } else {
          hist = new TH2F(histName.c_str(), "MC-IBD Neutron Capture Vtx y .vs. x in Others",
                          120, -6, 6, 120, -6, 6);
          hist->GetXaxis()->SetTitle("x (m)");
          hist->GetYaxis()->SetTitle("y (m)");
        }
        break;

      case ZVSYNCAPO:

        histName = "mc_zVSyNCaptureinOthers_IBD";
        
        if (siteInt!=3) {
          hist = new TH2F(histName.c_str(), "MC-IBD Neutron Capture Vtx z .vs. y in Others",
                          120, -6, 6, 70, 0, 7);
          hist->GetXaxis()->SetTitle("y [AD1: y-3; AD2: y+3](m)");
          hist->GetYaxis()->SetTitle("z [AD1,AD2: z+8](m)");
        } else {
          hist = new TH2F(histName.c_str(), "MC-IBD Neutron Caputre Vtx z .vs. y in Others",
                          120, -6, 6, 140, -7, 7);
          hist->GetXaxis()->SetTitle("y [AD2: y+6; AD3: y-6](m)");
          hist->GetYaxis()->SetTitle("z [AD1,AD2: z+8, AD3,AD4: z+2](m)");
        }

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

      case ZVSXNCAPO:

        histName = "mc_zVSxNCaptureinOthers_IBD";

        if (siteInt!=3) {
          hist = new TH2F(histName.c_str(), "MC-IBD Neutron Capture Vtx z .vs. x in Others",
                          120, -6, 6, 70, 0, 7);
          hist->GetXaxis()->SetTitle("x (m)");
          hist->GetYaxis()->SetTitle("z [AD1, AD2: z+8](m)");
        } else {
          hist = new TH2F(histName.c_str(), "MC-IBD Neutron Capture Vtx z .vs. x in Others",
                          120, -6, 6, 140, -7, 7);
          hist->GetXaxis()->SetTitle("x (m)");
          hist->GetYaxis()->SetTitle("z [AD1, AD2: z+8; AD3,AD4: z+2](m)");
        }
        break;  

      case TNABS:

        histName = "mc_TNAbsorption_IBD";
        hist = new TH1F(histName.c_str(),  "MC-IBD Neutron Absorption Time",
                        100, 0, 2000);

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

      case YVSXNABS:

        histName = "mc_yVSxNAbsorption_IBD";

        if (siteInt!=3) {
          hist = new TH2F(histName.c_str(), "MC-IBD Neutron Absorption Vtx y .vs. x",
                          120, -6, 6, 60, 0, 6);
          hist->GetXaxis()->SetTitle("x (m)");
          hist->GetYaxis()->SetTitle("y [AD1,AD2: y+3](m)");
        } else {
          hist = new TH2F(histName.c_str(), "MC-IBD Neutron Absorption Vtx y .vs. x",
                          120, -6, 6, 120, -6, 6);
          hist->GetXaxis()->SetTitle("x (m)");
          hist->GetYaxis()->SetTitle("y (m)");
        }
        break;
        
      case ZVSYNABS:

        histName = "mc_zVSyNAbsorption_IBD";

        if (siteInt!=3) {
          hist = new TH2F(histName.c_str(), "MC-IBD Neutron Absorption Vtx z .vs. y in Others",
                          120, -6, 6, 70, 0, 7);
          hist->GetXaxis()->SetTitle("y [AD1: y-3; AD2: y+3](m)");
          hist->GetYaxis()->SetTitle("z [AD1,AD2: z+8](m)");
        } else {
          hist = new TH2F(histName.c_str(), "MC-IBD Neutron Absorption Vtx z .vs. y in Others",
                          120, -6, 6, 140, -7, 7);
          hist->GetXaxis()->SetTitle("y [AD2: y+6; AD3: y-6](m)");
          hist->GetYaxis()->SetTitle("z [AD1,AD2: z+8, AD3,AD4: z+2](m)");
        }
        break;

      case ZVSXNABS:

        histName = "mc_zVSxNAbsorption_IBD";

        if (siteInt!=3) {
          hist = new TH2F(histName.c_str(), "MC-IBD Neutron Absorption Vtx z .vs. x",
                          120, -6, 6, 70, 0, 7);
          hist->GetXaxis()->SetTitle("x (m)");
          hist->GetYaxis()->SetTitle("z [AD1, AD2: z+8](m)");
        } else {
          hist = new TH2F(histName.c_str(), "MC-IBD Neutron Absorption Vtx z .vs. x",
                          120, -6, 6, 140, -7, 7);
          hist->GetXaxis()->SetTitle("x (m)");
          hist->GetYaxis()->SetTitle("z [AD1, AD2: z+8; AD3,AD4: z+2](m)");
        }
        break;
 
      case HITSVSQUEN:
        
        histName = "mc_simhitsVSquen_IBD";
        hist = new TH2F(histName.c_str(),
                        "MC Simhits/Quenched E .vs. Quen E",
                        100, 0,  15, 
                        100, 0, 300);
        hist->GetXaxis()->SetTitle("Quenched E (MeV)");
        hist->GetYaxis()->SetTitle("SumHits PEs / Quenched E");
        break;

      case QUENVSDEPO:

        histName = "mc_quenVSdepo_IBD";
        hist = new TH2F(histName.c_str(),
                        "MC Quenched/Deposit E .vs. Quen E",
                        100, 0, 15, 
                        100, 0, 1.2);
        hist->GetXaxis()->SetTitle("Quenched E (MeV)");
        hist->GetYaxis()->SetTitle("Quenched E / Deposited E");
        break;

      default:
        error() << "Unknown Histogram: " << histogram << endreq;
        return 0;
      }
    }

    debug() << "Making histogram: " << histName << endreq;
    m_statsSvc->put( this->getPath(run, histName.c_str()), 
                     hist );
    histograms[siteInt * NMCTRUTHIBDHISTOGRAMS
               + histogram] = hist;
  }

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

Definition at line 1093 of file McTruthFigIBD.cc.

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

Member Data Documentation

Definition at line 112 of file McTruthFigIBD.h.

Definition at line 115 of file McTruthFigIBD.h.

Definition at line 118 of file McTruthFigIBD.h.

Definition at line 119 of file McTruthFigIBD.h.

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

Definition at line 120 of file McTruthFigIBD.h.

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

Definition at line 121 of file McTruthFigIBD.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