In This Package:

GenHeaderCnv.cc

Go to the documentation of this file.

#include "GenHeaderCnv.h"
#include "PerBaseEvent/HeaderObjectCnv.h"

using namespace DayaBay;
using namespace std;

GenHeaderCnv::GenHeaderCnv(ISvcLocator* svc)
    : RootIOTypedCnv<PerGenHeader,GenHeader>("PerGenHeader",
                                             classID(),svc)
{
  // Set threshold for maximum number of vertices.
  // GenHeaders with more vertices will not be written out 
  // This is temporary for MDC10a
  // wangzhe
  // Raise it to 50 for safe. When adding muon spallation background, it is ofter above 10.
  m_maxVertices = 50;  
  // Wz 
  MsgStream log(msgSvc(), "GenHeaderCnv");

  log << MSG::INFO 
      << "Only minimal GenHeader info will be persisted for GenHeaders with > " 
      << m_maxVertices 
      << " vertices. THIS IS TEMPORARY FOR MDC10a."<<endreq;
}

GenHeaderCnv::~GenHeaderCnv()
{
}


// From Persistent to Transient
StatusCode GenHeaderCnv::PerToTran(const PerGenHeader& perobj,
                                   DayaBay::GenHeader& tranobj)
{
  StatusCode sc = HeaderObjectCnv::toTran(perobj,tranobj);
  if (sc.isFailure()) return sc;

  MsgStream log(msgSvc(), "GenHeaderCnv::PerToTran");

  tranobj.setGeneratorName(perobj.name);

  PerGenEvent* pgevt = perobj.event;
  if (!pgevt) {
    tranobj.setEvent(0);
    log << MSG::WARNING
        << "Null PerGenEvent pointer"
        << endreq;
    return StatusCode::SUCCESS;
  }

  HepMC::GenEvent* gevt = new HepMC::GenEvent(pgevt->signalProcessId, 
                                              pgevt->eventNumber,
                                              0,
                                              pgevt->weights);

  // Map from persistent to transient particles
  std::map<PerGenParticle*,HepMC::GenParticle*> particleMap;
  for (size_t ind=0; ind<pgevt->particles.size(); ++ind) {
    PerGenParticle* pgpart = pgevt->particles[ind];
    HepMC::GenParticle* gpart = 
      new HepMC::GenParticle(pgpart->momentum,
                             pgpart->pdgId,
                             pgpart->status,
                             HepMC::Flow(),
                             HepMC::Polarization(pgpart->polTheta,
                                                 pgpart->polPhi));

    particleMap[pgpart] = gpart;
  }

  // Convert vertices and add in/out particles
  std::map<PerGenVertex*,HepMC::GenVertex*> vertexMap;
  for (size_t ind=0; ind<pgevt->vertices.size(); ++ind) {
    PerGenVertex* pgvtx = pgevt->vertices[ind];
    HepMC::GenVertex* gvtx = new HepMC::GenVertex(pgvtx->position,
                                                  pgvtx->id,
                                                  pgvtx->weights);
    gevt->add_vertex(gvtx);
    for (size_t pi_ind=0; pi_ind < pgvtx->particlesIn.size(); ++pi_ind) {
      PerGenParticle* pgpart = pgevt->particles[pgvtx->particlesIn[pi_ind]];
      HepMC::GenParticle* gpart = particleMap[pgpart];
      gvtx->add_particle_in(gpart);
    }
    for (size_t po_ind=0; po_ind < pgvtx->particlesOut.size(); ++po_ind) {
      PerGenParticle* pgpart = pgevt->particles[pgvtx->particlesOut[po_ind]];
      HepMC::GenParticle* gpart = particleMap[pgpart];
      gvtx->add_particle_out(gpart);
    }

    if ((int)ind == pgevt->signalProcessIndex) {
      gevt->set_signal_process_vertex(gvtx);
    }
  }

  tranobj.setEvent(gevt);

  return StatusCode::SUCCESS;
}

// From Transient to Persistent
StatusCode GenHeaderCnv::TranToPer(const DayaBay::GenHeader& tranobj,
                                   PerGenHeader& perobj)
{
  StatusCode sc = HeaderObjectCnv::toPer(tranobj,perobj);
  if (sc.isFailure()) return sc;

  MsgStream log(msgSvc(), "GenHeaderCnv::TranToPer");

  perobj.name = tranobj.generatorName();

  const HepMC::GenEvent* genevt = tranobj.event();

  if (!genevt) {
    perobj.event = 0;
    log << MSG::WARNING
        << "Null HepMC::GenEvent pointer"
        << endreq;
    return StatusCode::SUCCESS;
  }

  PerGenEvent* pgevt = new PerGenEvent();

  // Fill the event

  // Basic stuff
  pgevt->eventNumber = genevt->event_number();
  pgevt->signalProcessId = genevt->signal_process_id();


  // Only save vertices when there are a few. 
  // This is temporary for MDC10a
  log << MSG::DEBUG << "# vertices = " << genevt->vertices_size() << endreq;
  if (genevt->vertices_size() <= m_maxVertices ){
    log << MSG::DEBUG << "Saving all vertices " << endreq;

    // Weights
    const HepMC::WeightContainer& weights = genevt->weights();
    pgevt->weights.clear();
    pgevt->weights.insert(pgevt->weights.begin(),
                          weights.begin(), weights.end());
    
    // Particles - first pass
    std::map<HepMC::GenParticle*,int> particleIndexMap;
    pgevt->particles.clear();
    HepMC::GenEvent::particle_const_iterator pit, 
      pdone = genevt->particles_end();
    for (pit = genevt->particles_begin(); pit != pdone; ++pit) {
      HepMC::GenParticle* gpart = *pit;
      particleIndexMap[gpart] = pgevt->particles.size(); // for later lookup
      CLHEP::HepLorentzVector momentum(gpart->momentum().x(),
                                       gpart->momentum().y(),
                                       gpart->momentum().z(),
                                       gpart->momentum().t());
      PerGenParticle* pgpart = 
        new PerGenParticle(momentum,
                           gpart->pdg_id(),
                           gpart->status(),
                           gpart->polarization().theta(),
                           gpart->polarization().phi(),
                           gpart->barcode());
      pgevt->particles.push_back(pgpart);
    }

    // Vertices
    pgevt->vertices.clear();
    HepMC::GenVertex* signal_process_vertex = genevt->signal_process_vertex();
    HepMC::GenEvent::vertex_const_iterator vit,
      vdone = genevt->vertices_end();
    for (vit = genevt->vertices_begin(); vit != vdone; ++vit) {
      HepMC::GenVertex* gvtx = *vit;
      CLHEP::HepLorentzVector position(gvtx->position().x(),
                                       gvtx->position().y(),
                                       gvtx->position().z(),
                                       gvtx->position().t());
      PerGenVertex* pgvtx = 
        new PerGenVertex(position,
                         gvtx->id(),
                         gvtx->barcode(),
                         gvtx->weights());
      
      // Fill in particle info
      HepMC::GenVertex::particles_in_const_iterator piit, 
        pidone = gvtx->particles_in_const_end();
      for (piit = gvtx->particles_in_const_begin(); piit != pidone; ++piit) {
        std::map<HepMC::GenParticle*,int>::iterator pimit = 
          particleIndexMap.find(*piit);
        if (pimit != particleIndexMap.end()) {
          pgvtx->particlesIn.push_back(pimit->second);
        }
      }
      // Fill out particle info
      HepMC::GenVertex::particles_out_const_iterator poit, 
        podone = gvtx->particles_out_const_end();
      for (poit = gvtx->particles_out_const_begin(); poit != podone; ++poit) {
        std::map<HepMC::GenParticle*,int>::iterator pomit = 
          particleIndexMap.find(*poit);
        if (pomit != particleIndexMap.end()) {
          pgvtx->particlesOut.push_back(pomit->second);
        }
      }
      
      if (signal_process_vertex == gvtx) {
        pgevt->signalProcessIndex = pgevt->vertices.size();
      }
      pgevt->vertices.push_back(pgvtx);
    }
  }
  else {
    log << MSG::DEBUG << "Not saving any vertices " << endreq;
  } //(genevt->vertices_size() <= m_maxVertices )
  
  delete perobj.event;
  perobj.event = pgevt;
  return StatusCode::SUCCESS;
}
    
StatusCode GenHeaderCnv::fillRepRefs(IOpaqueAddress*, DataObject* dobj)
{
    MsgStream log(msgSvc(), "GenHeaderCnv::fillRepRefs");
    GenHeader* gh = dynamic_cast<GenHeader*>(dobj);

    log << MSG::DEBUG
        << "Saving links to " << gh->inputHeaders().size() 
        << " input headers" << endreq;

    StatusCode sc = HeaderObjectCnv::fillPer(m_rioSvc,*gh,*m_perOutObj);
    if (sc.isFailure()) {
        log << MSG::ERROR << "Failed to fill HeaderObject part" << endreq;
        return sc;
    }

    // ... fill GenHeader part...
    return sc;
}

StatusCode GenHeaderCnv::fillObjRefs(IOpaqueAddress*, DataObject* dobj)
{
    MsgStream log(msgSvc(), "GenHeaderCnv::fillObjRefs");
    HeaderObject* hobj = dynamic_cast<HeaderObject*>(dobj);
    StatusCode sc = HeaderObjectCnv::fillTran(m_rioSvc,*m_perInObj,*hobj);
    if (sc.isFailure()) {
        log << MSG::ERROR << "Failed to fill HeaderObject part" << endreq;
        return sc;
    }

    log << MSG::DEBUG
        << "Restored links to " << hobj->inputHeaders().size() 
        << " input headers" << endreq;

    // ... fill GenHeader part...
    return sc;
}

---