In This Package:

AdExamplePlots Class Reference

Collaboration diagram for AdExamplePlots:[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< AdSummary >	adHistory

---

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

DayaBay::Detector AdExamplePlots::detector [private]

Definition at line 37 of file AdPlotterAlg.cc.

MsgStream& AdExamplePlots::log [private]

Definition at line 38 of file AdPlotterAlg.cc.

bool AdExamplePlots::firstReadout [private]

Definition at line 39 of file AdPlotterAlg.cc.

TimeStamp AdExamplePlots::lastTriggerTime [private]

Definition at line 40 of file AdPlotterAlg.cc.

vector<AdSummary> AdExamplePlots::adHistory [private]

Definition at line 41 of file AdPlotterAlg.cc.

TH2F* AdExamplePlots::genVertexZRho

Definition at line 46 of file AdPlotterAlg.cc.

TH1F* AdExamplePlots::nPhotoElectrons

Definition at line 48 of file AdPlotterAlg.cc.

TH1F* AdExamplePlots::hitPmts

Definition at line 50 of file AdPlotterAlg.cc.

TH1F* AdExamplePlots::dtReadout

Definition at line 51 of file AdPlotterAlg.cc.

TH1F* AdExamplePlots::dtReadoutShort

Definition at line 52 of file AdPlotterAlg.cc.

TH1F* AdExamplePlots::adcSpectrum

Definition at line 53 of file AdPlotterAlg.cc.

TH1F* AdExamplePlots::tdcSpectrum

Definition at line 54 of file AdPlotterAlg.cc.

TH1F* AdExamplePlots::chargeSum

Definition at line 55 of file AdPlotterAlg.cc.

TH2F* AdExamplePlots::hitPmtVsChargeSum

Definition at line 56 of file AdPlotterAlg.cc.

bool AdExamplePlots::findCoincidence

Definition at line 59 of file AdPlotterAlg.cc.

TH1F* AdExamplePlots::dtNeutronCapture

Definition at line 60 of file AdPlotterAlg.cc.

TH1F* AdExamplePlots::promptEnergy

Definition at line 61 of file AdPlotterAlg.cc.

TH1F* AdExamplePlots::delayedEnergy

Definition at line 62 of file AdPlotterAlg.cc.

TH2F* AdExamplePlots::delayedVsPromptEnergy

Definition at line 63 of file AdPlotterAlg.cc.

---

The documentation for this class was generated from the following file:

• /home/dayabay/installs/trunk_2011_04_11_opt/NuWa-trunk/dybgaudi/Validation/DataQuality/src/AdPlotterAlg.cc

---