/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 | Public Attributes | Private Attributes
AdExamplePlots Class Reference
Collaboration diagram for AdExamplePlots:
Collaboration graph
[legend]

List of all members.

Public Member Functions

 AdExamplePlots (DayaBay::Detector det, MsgStream &msg)
StatusCode book (ITHistSvc *hsvc, string filepath)
StatusCode fill (ICableSvc *, const ServiceMode &, ICoordSysSvc *coord, const ReadoutHeader *roh)

Public Attributes

TH2F * genVertexZRho
TH1F * nPhotoElectrons
TH1F * hitPmts
TH1F * dtReadout
TH1F * dtReadoutShort
TH1F * adcSpectrum
TH1F * tdcSpectrum
TH1F * chargeSum
TH2F * hitPmtVsChargeSum
bool findCoincidence
TH1F * dtNeutronCapture
TH1F * promptEnergy
TH1F * delayedEnergy
TH2F * delayedVsPromptEnergy

Private Attributes

DayaBay::Detector detector
MsgStream & log
bool firstReadout
TimeStamp lastTriggerTime
vector< AdSummaryadHistory

Detailed Description

Definition at line 36 of file AdPlotterAlg.cc.


Constructor & Destructor Documentation

AdExamplePlots::AdExamplePlots ( DayaBay::Detector  det,
MsgStream &  msg 
) [inline]

Definition at line 65 of file AdPlotterAlg.cc.

    : detector(det), 
      log(msg),
      firstReadout(true),
      findCoincidence(false)
  { }

Member Function Documentation

StatusCode AdExamplePlots::book ( ITHistSvc *  hsvc,
string  filepath 
) [inline]

Definition at line 72 of file AdPlotterAlg.cc.

                                                    {
    // Initialize and register all the histograms
    gStyle->SetHistLineWidth(2);
    gStyle->SetHistLineColor(4);
    {
      string name = detector.detName() + "_HitPmts";
      string title = "Number of Hit Pmts for " + detector.detName();
      hitPmts = new TH1F(name.c_str(),title.c_str(),200,0,200);
      hitPmts->GetXaxis()->SetTitle("Number of Hit PMTs");
      hitPmts->GetYaxis()->SetTitle("Readouts");
      name = filepath+"/"+detector.detName()+"/hitPmts";
      if (hsvc->regHist(name,hitPmts).isFailure()) {
        log << MSG::ERROR << "Could not register " << name << endreq;
        delete hitPmts; hitPmts = 0;
        return StatusCode::FAILURE;
      }
    }
    {
      string name = detector.detName() + "_DtReadout";
      string title = "Time between readouts for " + detector.detName();
      dtReadout = new TH1F(name.c_str(),title.c_str(),200,0,0.02);
      dtReadout->GetXaxis()->SetTitle("#Deltat [s]");
      dtReadout->GetYaxis()->SetTitle("Readouts / bin");
      name = filepath+"/"+detector.detName()+"/dtReadout";
      if (hsvc->regHist(name,dtReadout).isFailure()) {
        log << MSG::ERROR << "Could not register " << name << endreq;
        delete dtReadout; dtReadout = 0;
        return StatusCode::FAILURE;
      }
    }
    {
      string name = detector.detName() + "_DtReadoutShort";
      string title = "Short time between readouts for " + detector.detName();
      dtReadoutShort = new TH1F(name.c_str(),title.c_str(),200,0,100.e-6);
      dtReadoutShort->GetXaxis()->SetTitle("#Deltat [s]");
      dtReadoutShort->GetYaxis()->SetTitle("Readouts / bin");
      name = filepath+"/"+detector.detName()+"/dtReadoutShort";
      if (hsvc->regHist(name,dtReadoutShort).isFailure()) {
        log << MSG::ERROR << "Could not register " << name << endreq;
        delete dtReadoutShort; dtReadoutShort = 0;
        return StatusCode::FAILURE;
      }
    }
    {
      string name = detector.detName() + "_AdcSpectrum";
      string title = "ADC spectrum for all channels in " + detector.detName();
      adcSpectrum = new TH1F(name.c_str(),title.c_str(),4096,0,4096);
      adcSpectrum->GetXaxis()->SetTitle("ADC value");
      adcSpectrum->GetYaxis()->SetTitle("Entries");
      name = filepath+"/"+detector.detName()+"/adcSpectrum";
      if (hsvc->regHist(name,adcSpectrum).isFailure()) {
        log << MSG::ERROR << "Could not register " << name << endreq;
        delete adcSpectrum; adcSpectrum = 0;
        return StatusCode::FAILURE;
      }
    }
    {
      string name = detector.detName() + "_TdcSpectrum";
      string title = "TDC spectrum for all channels in " + detector.detName();
      tdcSpectrum = new TH1F(name.c_str(),title.c_str(),300,0,300);
      tdcSpectrum->GetXaxis()->SetTitle("TDC value");
      tdcSpectrum->GetYaxis()->SetTitle("Entries");
      name = filepath+"/"+detector.detName()+"/tdcSpectrum";
      if (hsvc->regHist(name,tdcSpectrum).isFailure()) {
        log << MSG::ERROR << "Could not register " << name << endreq;
        delete tdcSpectrum; tdcSpectrum = 0;
        return StatusCode::FAILURE;
      }
    }
    {
      string name = detector.detName() + "_ChargeSum";
      string title = "ADC Charge sum for " + detector.detName();
      chargeSum = new TH1F(name.c_str(),title.c_str(),500,0,20000);
      chargeSum->GetXaxis()->SetTitle("Sum(ADC - Baseline)");
      chargeSum->GetYaxis()->SetTitle("Readouts");
      name = filepath+"/"+detector.detName()+"/chargeSum";
      if (hsvc->regHist(name,chargeSum).isFailure()) {
        log << MSG::ERROR << "Could not register " << name << endreq;
        delete chargeSum; chargeSum = 0;
        return StatusCode::FAILURE;
      }
    }
    {
      string name = detector.detName() + "_HitPmtVsChargeSum";
      string title = "Number of Hit PMTs versus the Charge Sum for" 
        + detector.detName();
      hitPmtVsChargeSum = new TH2F(name.c_str(),title.c_str(),
                                   200,0,20000,
                                   200,0,200);
      hitPmtVsChargeSum->GetXaxis()->SetTitle("Sum(ADC - Baseline)");
      hitPmtVsChargeSum->GetYaxis()->SetTitle("Number of Hit PMTs");
      name = filepath+"/"+detector.detName()+"/hitPmtVsChargeSum";
      if (hsvc->regHist(name,hitPmtVsChargeSum).isFailure()) {
        log << MSG::ERROR << "Could not register " << name << endreq;
        delete hitPmtVsChargeSum; hitPmtVsChargeSum = 0;
        return StatusCode::FAILURE;
      }
    }
    // Simulated Hits
    {
      string name = detector.detName() + "_NPhotoElectrons";
      string title = "Number of simulated photoelectrons for " 
        + detector.detName();
      nPhotoElectrons = new TH1F(name.c_str(),title.c_str(),1000,0,1000);
      nPhotoElectrons->GetXaxis()->SetTitle("Number of Photoelectrons");
      nPhotoElectrons->GetYaxis()->SetTitle("Readouts");
      name = filepath+"/"+detector.detName()+"/nPhotoElectrons";
      if (hsvc->regHist(name,nPhotoElectrons).isFailure()) {
        log << MSG::ERROR << "Could not register " << name << endreq;
        delete nPhotoElectrons; nPhotoElectrons = 0;
        return StatusCode::FAILURE;
      }
    }
    // Generated Vertices
    {
      string name = detector.detName() + "_GenVertexZRho";
      string title = "Position of generated vertices for " + detector.detName();
      genVertexZRho = new TH2F(name.c_str(),title.c_str(),
                               200,0,3,
                               200,-5,5);
      genVertexZRho->GetXaxis()->SetTitle("Particle (#rho / 3 m)^{2}");
      genVertexZRho->GetYaxis()->SetTitle("Particle z [m]");
      name = filepath+"/"+detector.detName()+"/genVertexZRho";
      if (hsvc->regHist(name,genVertexZRho).isFailure()) {
        log << MSG::ERROR << "Could not register " << name << endreq;
        delete genVertexZRho; genVertexZRho = 0;
        return StatusCode::FAILURE;
      }
    }
    if(findCoincidence){
      // add coincidence histograms
      {
        // IBD Neutron Capture Time
        string name = detector.detName() + "_DtNeutronCapture";
        string title = "IBD Neutron capture spectrum for " + detector.detName();
        dtNeutronCapture = new TH1F(name.c_str(),title.c_str(),200,0,100.e-6);
        dtNeutronCapture->GetXaxis()->SetTitle("#Deltat [s]");
        dtNeutronCapture->GetYaxis()->SetTitle("Readouts");
        name = filepath+"/"+detector.detName()+"/dtNeutronCapture";
        if (hsvc->regHist(name,dtNeutronCapture).isFailure()) {
          log << MSG::ERROR << "Could not register " << name << endreq;
          delete dtNeutronCapture; dtNeutronCapture = 0;
          return StatusCode::FAILURE;
        }       
      }
      {
        // IBD Prompt Energy
        string name = detector.detName() + "_PromptEnergy";
        string title = "IBD Prompt Energy for " + detector.detName();
        promptEnergy = new TH1F(name.c_str(),title.c_str(),200,0,12);
        promptEnergy->GetXaxis()->SetTitle("Prompt Energy [MeV]");
        promptEnergy->GetYaxis()->SetTitle("Readouts");
        name = filepath+"/"+detector.detName()+"/promptEnergy";
        if (hsvc->regHist(name,promptEnergy).isFailure()) {
          log << MSG::ERROR << "Could not register " << name << endreq;
          delete promptEnergy; promptEnergy = 0;
          return StatusCode::FAILURE;
        }       
      }
      {
        // IBD Delayed Energy
        string name = detector.detName() + "_DelayedEnergy";
        string title = "IBD Delayed Energy for " + detector.detName();
        delayedEnergy = new TH1F(name.c_str(),title.c_str(),200,0,12);
        delayedEnergy->GetXaxis()->SetTitle("Delayed Energy [MeV]");
        delayedEnergy->GetYaxis()->SetTitle("Readouts");
        name = filepath+"/"+detector.detName()+"/delayedEnergy";
        if (hsvc->regHist(name,delayedEnergy).isFailure()) {
          log << MSG::ERROR << "Could not register " << name << endreq;
          delete delayedEnergy; delayedEnergy = 0;
          return StatusCode::FAILURE;
        }       
      }
      {
        // IBD Delayed vs. Prompt Energy
        string name = detector.detName() + "_DelayedVsPromptEnergy";
        string title = "IBD Delayed vs. Prompt Energy for " + detector.detName();
        delayedVsPromptEnergy = new TH2F(name.c_str(),title.c_str(),
                                         200,0,12,
                                         200,0,12);
        delayedVsPromptEnergy->GetXaxis()->SetTitle("Prompt Energy [MeV]");
        delayedVsPromptEnergy->GetYaxis()->SetTitle("Delayed Energy [MeV]");
        name = filepath+"/"+detector.detName()+"/delayedVsPromptEnergy";
        if (hsvc->regHist(name,delayedVsPromptEnergy).isFailure()) {
          log << MSG::ERROR << "Could not register " << name << endreq;
          delete delayedVsPromptEnergy; delayedVsPromptEnergy = 0;
          return StatusCode::FAILURE;
        }       
      }
    }
    return StatusCode::SUCCESS;
  }
StatusCode AdExamplePlots::fill ( ICableSvc ,
const ServiceMode ,
ICoordSysSvc coord,
const ReadoutHeader roh 
) [inline]

Definition at line 265 of file AdPlotterAlg.cc.

  {
    
    const DayaBay::ReadoutPmtCrate* ro = 
      dynamic_cast<const DayaBay::ReadoutPmtCrate*>(roh->readout());
    if (!ro) {
      log << MSG::DEBUG << "Got a non-PmtCrate readout, skipping..." << endreq;
      return StatusCode::SUCCESS;
    }
  
    double adcChargeSum = 0;
    int nHitPmts = 0;
    const ReadoutPmtCrate::PmtChannelReadouts& 
      channelMap = ro->channelReadout();
    ReadoutPmtCrate::PmtChannelReadouts::const_iterator
      it, done = channelMap.end();
    for (it=channelMap.begin(); it != done; ++it) {
      // Ignore calibration PMTs
      if(it->first.board() > 12) continue;
      nHitPmts++;
    
      // ADC data
      double baseline = 500.;
      float charge = 0;
      const vector<int> &adc = it->second.adc();
      vector<int>::const_iterator cit, cdone=adc.end();
      for (cit=adc.begin(); cit != cdone; ++cit) {
        adcSpectrum->Fill(*cit);
        charge += *cit - baseline;
      }
      adcChargeSum += charge;      

      // TDC data
      const vector<int> &tdc = it->second.tdc();
      vector<int>::const_iterator tdcIter, tdcDone = tdc.end();
      for (tdcIter=tdc.begin(); tdcIter != tdcDone; ++tdcIter) {
        tdcSpectrum->Fill(*tdcIter);
      }

    }
    hitPmts->Fill(nHitPmts);
    chargeSum->Fill(adcChargeSum);
    hitPmtVsChargeSum->Fill(adcChargeSum, nHitPmts);
  
    // Time between readouts
    if( !firstReadout ){
      TimeStamp dt = ro->triggerTime() - lastTriggerTime;
      dtReadout->Fill( dt.GetSeconds() );
      dtReadoutShort->Fill( dt.GetSeconds() );
    }else{
      firstReadout = false;
    }
    lastTriggerTime = ro->triggerTime();

    // If Simulation headers are available, fill additional figures
    const SimReadoutHeader* simReadoutHeader = 0;
    const ElecHeader* elecHeader = 0;
    vector<const SimHeader*> simHeaders;
    vector<const GenHeader*> genHeaders;
    vector<const IHeader*> inputHeaders = roh->inputHeaders();
    vector<const IHeader*>::const_iterator headerIter, 
      headerDone = inputHeaders.end();
    log << MSG::INFO << "Processing Input Headers: " << inputHeaders.size() 
        << endreq;
    for (headerIter=inputHeaders.begin(); headerIter != headerDone; 
         ++headerIter) {
      simReadoutHeader = dynamic_cast<const SimReadoutHeader*>(*headerIter);
      if(simReadoutHeader){
        log << MSG::INFO << " Found SimReadoutHeader" << endreq;
        break;
      }
    }
    if(simReadoutHeader){
      inputHeaders = simReadoutHeader->inputHeaders();
      headerDone = inputHeaders.end();
      for (headerIter=inputHeaders.begin(); headerIter != headerDone; 
           ++headerIter) {
        elecHeader = dynamic_cast<const ElecHeader*>(*headerIter);
        if(elecHeader){
          log << MSG::INFO << " Found ElecHeader" << endreq;
          break;
        }
      }
    }
    if(elecHeader){
      inputHeaders = elecHeader->inputHeaders();
      headerDone = inputHeaders.end();
      for (headerIter=inputHeaders.begin(); headerIter != headerDone; 
         ++headerIter) {
        const SimHeader* simHeader 
          = dynamic_cast<const SimHeader*>(*headerIter);
        if(simHeader){
          log << MSG::INFO << " Found SimHeader" << endreq;
          simHeaders.push_back(simHeader);
        }
      }
    }
    if(simHeaders.size() > 0){
      for(unsigned int simHeaderIdx=0; simHeaderIdx<simHeaders.size(); 
          simHeaderIdx++){
        inputHeaders = (simHeaders[simHeaderIdx])->inputHeaders();
        headerDone = inputHeaders.end();
        for (headerIter=inputHeaders.begin(); headerIter != headerDone; 
             ++headerIter) {
          const GenHeader* genHeader 
            = dynamic_cast<const GenHeader*>(*headerIter);
          if(genHeader){
            log << MSG::INFO << " Found GenHeader" << endreq;
            genHeaders.push_back(genHeader);
          }
        }
      }
    }
    
    // Fill additional histograms
    // Loop over simulated hits
    for(unsigned int simHeaderIdx=0; simHeaderIdx<simHeaders.size(); 
        simHeaderIdx++){
      log << MSG::INFO << " simHeaderIdx = " << simHeaderIdx << endreq;
      const SimHeader* simHeader = simHeaders[simHeaderIdx];
      const SimHitHeader* shh = simHeader->hits();
      if(!shh){
        log << MSG::INFO << "Can't find hit header." << endreq;
        continue;
      }
      short int detId = detector.siteDetPackedData();
      SimHitHeader::hc_map::const_iterator hcmapIter 
        = shh->hitCollection().find(detId);
      if(hcmapIter == shh->hitCollection().end()){
        log << MSG::INFO << "Can't find hit collection for detector " 
            << detId << endreq;
        log << MSG::INFO << "Have hit collections for " 
            << shh->hitCollection().size() << " detectors: "; 
        SimHitHeader::hc_map::const_iterator hcmapDone 
          = shh->hitCollection().end();
        for(hcmapIter = shh->hitCollection().begin(); hcmapIter != hcmapDone;
            hcmapIter++){
          log << MSG::INFO << hcmapIter->first << "\t";
        }
        continue;
      }
      const SimHitCollection* hc = hcmapIter->second;
      if(!hc) continue;
      int nPeHits = 0;
      const SimHitCollection::hit_container& hits = hc->collection();
      SimHitCollection::hit_container::const_iterator hitIter, 
        hitEnd = hits.end();
      for(hitIter = hits.begin(); hitIter != hitEnd; hitIter++){
        const SimHit* hit = *hitIter;
        if(!hit) continue;
        // Ignore calibration pmts
        if(AdPmtSensor(hit->sensDetId()).ring()<1) continue;
        nPeHits++;
      }
      log << MSG::INFO << " Npe = " << nPeHits << endreq;
      nPhotoElectrons->Fill(nPeHits);
    }
    
    for(unsigned int genHeaderIdx=0; genHeaderIdx<genHeaders.size(); 
        genHeaderIdx++){
      const GenHeader* genHeader = genHeaders[genHeaderIdx];
      const HepMC::GenEvent* event = genHeader->event();
      if(!event) continue;
      HepMC::GenEvent::vertex_const_iterator vtxIter, 
        vtxDone = event->vertices_end();
      for(vtxIter = event->vertices_begin(); vtxIter != vtxDone; vtxIter++){
        const HepMC::GenVertex* vtx = *vtxIter;
        if(!vtx) continue;
        double x = vtx->position().x();
        double y = vtx->position().y();
        double z = vtx->position().z();
        log << MSG::INFO << "\n" << vtx->position().t() << endreq;
        log << MSG::INFO << " Vertex time: " << vtx->position().t() << endreq;
        log << MSG::INFO << "\n" << vtx->position().t() << endreq;
        Gaudi::XYZPoint globalPoint(x,y,z);
        IDetectorElement* coordDE = coord->coordSysDE(globalPoint);
        if(!coordDE){
          log << MSG::ERROR << "Failed to find detector element for point: "
              << x << "\t" << y << "\t" << z 
              << endreq;
          continue;
        }
        Gaudi::XYZPoint localPoint = coordDE->geometry()->toLocal(globalPoint);
        double rho_3meters = localPoint.Rho()/(3.*Gaudi::Units::meter);
        double rhoSq = rho_3meters * rho_3meters;
        double z_meters = localPoint.Z()/Gaudi::Units::meter;
        genVertexZRho->Fill(rhoSq, z_meters);
        if( rhoSq < 0 || rhoSq > 3. || z_meters>5. || z_meters<-5. ){
          log << MSG::VERBOSE << "vtx position: " << localPoint.x() << "\t" 
              << localPoint.y() << "\t" << localPoint.z() 
              << endreq;
        }
        // Example loop over particles in vertex
        /*
        HepMC::GenVertex::particles_out_const_iterator pIter, 
          pDone = vtx->particles_out_const_end();
        for(pIter = vtx->particles_out_const_begin(); pIter != pDone; pIter++){
          const HepMC::GenParticle* particle = *pIter;
          if(!particle) continue;
          log << MSG::VERBOSE << "particle energy: " 
              << particle->momentum().e() << endreq;    
        }
        */
        // End example
      }
    }

    if(findCoincidence){
      double dtCut_s = 100.0e-6; 
      double epCutLow = 0.7;
      double epCutHigh = 8.0;
      double edCutLow = 6.0;
      double edCutHigh = 10.0;
      double chargeScale = 14000./8.;
      double ed = adcChargeSum / chargeScale;

      vector<AdSummary> newHistory;
      vector<AdSummary>::iterator histIter, histEnd = adHistory.end();
      for(histIter = adHistory.begin(); histIter != histEnd; histIter++){
        double dt = ro->triggerTime() - histIter->triggerTime;
        if(dt < dtCut_s){
          AdSummary promptSummary = *histIter;
          newHistory.push_back(promptSummary);
          double ep = promptSummary.adcSum / chargeScale;
          dtNeutronCapture->Fill(dt);
          if( ed > edCutLow && ed < edCutHigh ) promptEnergy->Fill(ep);
          if( ep > epCutLow && ep < epCutHigh ) delayedEnergy->Fill(ed);
          delayedVsPromptEnergy->Fill(ep,ed);
        }
      }
      AdSummary currentSummary;
      currentSummary.triggerTime = ro->triggerTime();
      currentSummary.adcSum = adcChargeSum;
      newHistory.push_back(currentSummary);
      adHistory = newHistory;
    }

    return StatusCode::SUCCESS;
  }

Member Data Documentation

Definition at line 37 of file AdPlotterAlg.cc.

MsgStream& AdExamplePlots::log [private]

Definition at line 38 of file AdPlotterAlg.cc.

Definition at line 39 of file AdPlotterAlg.cc.

Definition at line 40 of file AdPlotterAlg.cc.

Definition at line 41 of file AdPlotterAlg.cc.

Definition at line 46 of file AdPlotterAlg.cc.

Definition at line 48 of file AdPlotterAlg.cc.

Definition at line 50 of file AdPlotterAlg.cc.

Definition at line 51 of file AdPlotterAlg.cc.

Definition at line 52 of file AdPlotterAlg.cc.

Definition at line 53 of file AdPlotterAlg.cc.

Definition at line 54 of file AdPlotterAlg.cc.

Definition at line 55 of file AdPlotterAlg.cc.

Definition at line 56 of file AdPlotterAlg.cc.

Definition at line 59 of file AdPlotterAlg.cc.

Definition at line 60 of file AdPlotterAlg.cc.

Definition at line 61 of file AdPlotterAlg.cc.

Definition at line 62 of file AdPlotterAlg.cc.

Definition at line 63 of file AdPlotterAlg.cc.


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

Generated on Fri May 16 2014 10:06:20 for DataQuality by doxygen 1.7.4