In This Package:

Ana15.cc

Go to the documentation of this file.

#include "Ana15.h"
#include "Stage/HeaderStageData.h"

// for GenHeader
#include "Event/GenHeader.h"

// for GenEvent
#include "HepMC/GenEvent.h"

// detector element
#include "DetDesc/IGeometryInfo.h"
#include "DetDesc/DetectorElement.h"
#include "DetDesc/SolidBase.h"

// for SimHit
#include "Event/SimHeader.h"
#include "Event/SimHitHeader.h"
#include "Event/SimHitCollection.h"
#include "Event/SimPmtHit.h"
#include "Event/SimHit.h"

// for ElecEvent
#include "Event/ElecHeader.h"

// SimTrig
#include "Event/SimTrigHeader.h"

// SimReadout
#include "Event/SimReadoutHeader.h"

// for Readout
#include "Event/ReadoutHeader.h"
#include "Event/ReadoutPmtCrate.h"
#include "Event/ReadoutRpcCrate.h"

// for history info

//#include "Event/SimParticleHistoryHeader.h"
#include "Event/SimParticleHistory.h"
#include "Event/SimUnobservableStatisticsHeader.h"
#include "Event/SimStatistic.h"

// for Detector ID etc.
#include "Conventions/Detectors.h"

#include "GaudiKernel/ITHistSvc.h"

#include <Event/HeaderObject.h>

using namespace std;
using namespace DayaBay;

Ana15::Ana15(const std::string& name, ISvcLocator* pSvcLocator):
  GaudiAlgorithm(name, pSvcLocator)
{
  declareProperty("TopStage",m_TopStage="SingleLoader","Name of top stage");  
}

Ana15::~Ana15()
{
}

StatusCode Ana15::initialize()
{
  this->GaudiAlgorithm::initialize();

  // init path
  if(m_TopStage=="SingleLoader") m_path = ReadoutHeader::defaultLocation();
  if(m_TopStage=="TrigRead")     m_path = SimReadoutHeader::defaultLocation();
  if(m_TopStage=="Electronic")   m_path = ElecHeader::defaultLocation();
  if(m_TopStage=="Detector")     m_path = SimHeader::defaultLocation();
  if(m_TopStage=="Kinematic")    m_path = GenHeader::defaultLocation();
  
  // init root tree
  m_valiTree=new ValidationTree;
  m_valiTree->create();

  StatusCode status;

  return StatusCode::SUCCESS;
}

StatusCode Ana15::execute()
{
  double m_e=5.109989e-1;  // MeV
  double m_alpha=3727.379; // MeV

  m_valiTree->reset();

  StatusCode   sc;
  DataObject   *pObject;
  sc=eventSvc()->retrieveObject(m_path,pObject);     // get a new data
  if(sc.isFailure()) {
    error()<<"Failed to read data"<<endreq;
    return sc;
  }

  const HeaderObject     *pHeader=0;
  pHeader=dynamic_cast<HeaderObject*>(pObject);

  info()<<"to grep: new data arrived at time: "<<pHeader->earliest()<<endreq;

  const GenHeader        *pGenHeader=0;
  const SimHeader        *pSimHeader=0;
  const ElecHeader       *pElecHeader=0;
  const SimTrigHeader    *pSimTrigHeader=0;
  const SimReadoutHeader *pSimReadoutHeader=0;
  const ReadoutHeader    *pReadoutHeader=0;

  const IHeader* firstHeader=0;

  // get the header object pointer of the top stage
  if(m_TopStage=="SingleLoader") {
    // ReadoutHeader
    pReadoutHeader=dynamic_cast<ReadoutHeader*>(pObject);
  } else if(m_TopStage=="TrigRead") { 
    // SimReadoutHeader
    pSimReadoutHeader=dynamic_cast<SimReadoutHeader*>(pObject);
  } else if(m_TopStage=="Electronic") {
    // ElecHeader
    pElecHeader=dynamic_cast<ElecHeader*>(pObject);
  } else if (m_TopStage=="Detector") {
    // SimHeader
    pSimHeader=dynamic_cast<SimHeader*>(pObject);
  } else if (m_TopStage=="Kinematic") {
    // GenHeader
    pGenHeader=dynamic_cast<GenHeader*>(pObject);
  }

  // get the header object pointer of it lower stages
  // only the first one, for low event rate, it is enough
  if(pReadoutHeader!=0) {
    // SimReadoutHeader
    if(pReadoutHeader->inputHeaders().size() != 0) {
      firstHeader=*((pReadoutHeader->inputHeaders()).begin());
      pSimReadoutHeader= dynamic_cast<const SimReadoutHeader*>(firstHeader);
    }
  }
  if(pSimReadoutHeader!=0) {
    debug() << "Number of input SimTrigHeaders:" << pSimReadoutHeader->inputHeaders().size() << endreq;
    // SimTrigHeader
    if(pSimReadoutHeader->inputHeaders().size() != 0) {
      firstHeader=*((pSimReadoutHeader->inputHeaders()).begin());
      pSimTrigHeader= dynamic_cast<const SimTrigHeader*>(firstHeader);
    }
  }
  if(pSimTrigHeader!=0) {
    // ElecHeader
    if(pSimTrigHeader->inputHeaders().size() != 0) {
      firstHeader=*((pSimTrigHeader->inputHeaders()).begin());
      pElecHeader= dynamic_cast<const ElecHeader*>(firstHeader);
    }
  }
  if(pElecHeader!=0) {
    // SimHeader
    if(pElecHeader->inputHeaders().size() != 0) {
      firstHeader=*((pElecHeader->inputHeaders()).begin());
      pSimHeader= dynamic_cast<const SimHeader*>(firstHeader);
    }
  }
  if(pSimHeader!=0) {
    // GenHeader
    if(pSimHeader->inputHeaders().size() != 0) {
      firstHeader=*((pSimHeader->inputHeaders()).begin());
      pGenHeader= dynamic_cast<const GenHeader*>(firstHeader);
    }
  }


  // event information
  m_valiTree->time=pHeader->earliest().GetSeconds();

  // check ReadoutHeader information
  if(pReadoutHeader!=0) {
    m_valiTree->execSing=pReadoutHeader->execNumber();
    
    // adc tdc
    const Readout* pReadout = pReadoutHeader->readout();
    Detector det = pReadout->detector();
    
    if(det.detectorId() == DetectorId::kAD1 ||  // AD information
       det.detectorId() == DetectorId::kAD2 ||
       det.detectorId() == DetectorId::kAD3 ||
       det.detectorId() == DetectorId::kAD4) {
      
      // convert it to PMT Crate Readout
      const ReadoutPmtCrate* pmtcrate = 0;
      pmtcrate = dynamic_cast<const ReadoutPmtCrate*>(pReadout);
      
      // channnel map
      const ReadoutPmtCrate::PmtChannelReadouts& channelMap
          = pmtcrate->channelReadout();
      
      ReadoutPmtCrate::PmtChannelReadouts::const_iterator cmci;

      for(cmci=channelMap.begin();cmci!=channelMap.end();++cmci) {
        const FeeChannelId&      channelId = cmci->first;
        const ReadoutPmtChannel& aChannel  = cmci->second;

        const vector<int>&          tdc=aChannel.tdc();
        const vector<int>&          adc=aChannel.adc();

        int raw_tdc=tdc[0];
        int raw_adc=adc[0];
        
        m_valiTree->adc_sum+=raw_adc;

        verbose()<<channelId<<" tdc "<<raw_tdc<<" adc "<<raw_adc<<endreq;
      }
    }  // end of AD crate
  } // end of pReadoutHeader

  // check SimReadoutHeader information
  if(pSimReadoutHeader!=0) {
    m_valiTree->execTrig=pSimReadoutHeader->execNumber();

    const SimReadoutHeader::SimReadoutContainer 
      simReadout=pSimReadoutHeader->readouts();
    debug()<<"No. of SimReadout: "<<simReadout.size()<<endreq;
    debug()<<pSimReadoutHeader->execNumber()<<endreq;
  } // end of pSimReadoutHeader

  // check electrontic simulation information if available
  if(pElecHeader!=0) {
    m_valiTree->execElec=pElecHeader->execNumber();
  } // end of pElecHeader

  // check trigger header information if available
  if(pSimTrigHeader!=0) {
    debug()<<pSimTrigHeader->commandHeader()->collections().size();
  }

  // check detector simulation information if available
  if(pSimHeader!=0) {
    m_valiTree->execDets=pSimHeader->execNumber();
    const SimHitHeader* simhitheader = pSimHeader->hits();
    const SimHitHeader::hc_map& hcmap = simhitheader->hitCollection();
    SimHitHeader::hc_map::const_iterator it;

    m_valiTree->SimHitsEntries=0;
    int index(0);

    debug()<<"size of hit collection "<< hcmap.size()<<endreq;
    for (it=hcmap.begin(); it != hcmap.end(); ++it) {

      Detector det(it->first);
      debug() << "Got hit collection from " << det.detName()<<" id= "
              << det.siteDetPackedData()<<endreq;

      const std::vector<SimHit*>& hitvec=it->second->collection();
      std::vector<SimHit*>::const_iterator it_hvec;
      
      debug() <<"size of hit vector "<< hitvec.size()<<endreq;
      for (it_hvec=hitvec.begin(); it_hvec!=hitvec.end(); ++it_hvec) {
        
        index=m_valiTree->SimHitsEntries;
        if(index<MaxSimHitEntries){

          m_valiTree->hitTime[index]  = (*it_hvec)->hitTime();
          m_valiTree->hitx[index]     = (*it_hvec)->localPos().x();
          m_valiTree->hity[index]     = (*it_hvec)->localPos().y();
          m_valiTree->hitz[index]     = (*it_hvec)->localPos().z();
          m_valiTree->sensID[index]   = (*it_hvec)->sensDetId();
          m_valiTree->weight[index]   = (*it_hvec)->weight();
          
          // optical photon's ancestor's pdg
          if((*it_hvec)->ancestor().track()) {
            m_valiTree->ancestorPdg[index] = (*it_hvec)->ancestor().track()->particle();
          }
          else {
            m_valiTree->ancestorPdg[index] = 0;
          }
          //optical photon's wavelength
          const SimPmtHit* pmthit = static_cast<const SimPmtHit*>(*it_hvec);
          m_valiTree->wavelength[index]   = (*pmthit).wavelength();
          m_valiTree->polx[index]   = (*pmthit).pol().x();
          m_valiTree->poly[index]   = (*pmthit).pol().y();
          m_valiTree->polz[index]   = (*pmthit).pol().z();
          m_valiTree->px[index]   = (*pmthit).dir().x();
          m_valiTree->py[index]   = (*pmthit).dir().y();
          m_valiTree->pz[index]   = (*pmthit).dir().z();
          
          m_valiTree->SimHitsEntries = index+1;
          
        }
        else {
          debug()<< " Reached the max hit limit!!!!!" <<endl;
        }
      }
    }
  } // end of pSimHeader

  //  if(0==1) {
  if(pGenHeader!=0) {
    m_valiTree->execKine=pGenHeader->execNumber();
    debug()<<"Generator name: "<<pGenHeader->generatorName()<<endreq;
    debug()<<"execNum: "<<pGenHeader->execNumber()<<endreq;

    const HepMC::GenEvent* event=pGenHeader->event();
    //  event->print();

    HepMC::GenEvent::particle_const_iterator p;
    int count=0;
    for ( p = event->particles_begin(); p!=event->particles_end(); ++p ) {      
      (*p)->print();
      if((*p)->production_vertex()) {   // only outgoing particle has a production vertex
        count++;

        if(count==1) {
          
          m_valiTree->trk1_pdgId=(*p)->pdg_id();
          
          m_valiTree->trk1_px=(*p)->momentum().px();
          m_valiTree->trk1_py=(*p)->momentum().py();
          m_valiTree->trk1_pz=(*p)->momentum().pz();
          m_valiTree->trk1_e =(*p)->momentum().e();
          
          m_valiTree->trk1_x=(*p)->production_vertex()->position().x();
          m_valiTree->trk1_y=(*p)->production_vertex()->position().y();
          m_valiTree->trk1_z=(*p)->production_vertex()->position().z();
          m_valiTree->trk1_t=(*p)->production_vertex()->position().t();

          debug()<<"trk1_pdgId: "<<m_valiTree->trk1_pdgId<<endreq;
          
          if(m_valiTree->trk1_pdgId==22) { // gamma
            m_valiTree->trk1_ke=m_valiTree->trk1_e;
            m_valiTree->trk1_ve=m_valiTree->trk1_e;
          }else if (m_valiTree->trk1_pdgId==11) { // e-
            m_valiTree->trk1_ke=m_valiTree->trk1_e-m_e;
            m_valiTree->trk1_ve=m_valiTree->trk1_e-m_e;
          }else if (m_valiTree->trk1_pdgId==-11) { // e+
            m_valiTree->trk1_ke=m_valiTree->trk1_e-m_e;
            m_valiTree->trk1_ve=m_valiTree->trk1_e+m_e;
          }else if (m_valiTree->trk1_pdgId==1000020040) { // alpha
            m_valiTree->trk1_ke=m_valiTree->trk1_e-m_alpha;
            m_valiTree->trk1_ve=m_valiTree->trk1_e-m_alpha;
          }
        }  // outgoing particle 1
        if(count==2) {
          m_valiTree->trk2_pdgId=(*p)->pdg_id();
          
          m_valiTree->trk2_px=(*p)->momentum().px();
          m_valiTree->trk2_py=(*p)->momentum().py();
          m_valiTree->trk2_pz=(*p)->momentum().pz();
          m_valiTree->trk2_e =(*p)->momentum().e();
          
          m_valiTree->trk2_x=(*p)->production_vertex()->position().x();
          m_valiTree->trk2_y=(*p)->production_vertex()->position().y();
          m_valiTree->trk2_z=(*p)->production_vertex()->position().z();
          m_valiTree->trk2_t=(*p)->production_vertex()->position().t();
          if(m_valiTree->trk2_pdgId==22) { // gamma
            m_valiTree->trk2_ke=m_valiTree->trk2_e;
            m_valiTree->trk2_ve=m_valiTree->trk2_e;
          }else if (m_valiTree->trk2_pdgId==11) { // e-
            m_valiTree->trk2_ke=m_valiTree->trk2_e-m_e;
            m_valiTree->trk2_ve=m_valiTree->trk2_e-m_e;
          }else if (m_valiTree->trk2_pdgId==-11) { // e+
            m_valiTree->trk2_ke=m_valiTree->trk2_e-m_e;
            m_valiTree->trk2_ve=m_valiTree->trk2_e+m_e;
          }else if (m_valiTree->trk2_pdgId==1000020040) { // alpha
            m_valiTree->trk2_ke=m_valiTree->trk2_e-m_alpha;
            m_valiTree->trk2_ve=m_valiTree->trk2_e-m_alpha;
          }
        }  // outgoing particle 2
      } // end of outgoing particles
    }
  } // end of pGenHeader

  m_valiTree->fill();
  return StatusCode::SUCCESS;
}

StatusCode Ana15::finalize()
{
  m_valiTree->close();
  delete m_valiTree;

  return this->GaudiAlgorithm::finalize();
}

---