/data/refman/pasoursint/CMSSW_4_1_8_patch9/src/IOMC/NtupleConverter/src/Ntuple2HepMCFiller.cc

Go to the documentation of this file.

#include "IOMC/NtupleConverter/interface/Ntuple2HepMCFiller.h"

#include<iostream>
#include "IOMC/NtupleConverter/interface/NtupleROOTFile.h"  
using namespace std;
using namespace HepMC;

Ntuple2HepMCFiller * Ntuple2HepMCFiller::instance_=0;
//-------------------------------------------------------------------------
Ntuple2HepMCFiller * Ntuple2HepMCFiller::instance(){    
        
        if (instance_== 0) {
                instance_ = new Ntuple2HepMCFiller();
        }
        return instance_;
}


Ntuple2HepMCFiller::Ntuple2HepMCFiller(): 
   initialized_(false), input_(0),
   index_to_particle(3996),evtid(1),ntpl_id(1) { 
        cout << "Constructing a new Ntuple2HepMCFiller" << endl;
        if(instance_ == 0) instance_ = this;    
} 

//-------------------------------------------------------------------------
Ntuple2HepMCFiller::~Ntuple2HepMCFiller(){ 
        cout << "Destructing Ntuple2HepMCFiller" << endl;       
        if(input_) {
                delete input_;
        }
        instance_=0;    
}
//-------------------------------------------------------------------------

void Ntuple2HepMCFiller::setInitialized(bool value){initialized_=value;}

//-------------------------------------------------------------------------
void Ntuple2HepMCFiller::initialize(const string & filename, int id){
        
        if (initialized_) {
                cout << "Ntuple2HepMCFiller was already initialized... reinitializing it " << endl;
                if(input_) {
                        delete input_;
                }               
        }       
        
        cout<<"Ntuple2HepMCFiller::initialize : Opening file "<<filename<<endl;
        ntpl_id=id;
        input_ = new NtupleROOTFile( filename, id);
        initialized_ = true;   
}

//-------------------------------------------------------------------------

bool Ntuple2HepMCFiller::isInitialized(){
        
        return initialized_;
}

//-------------------------------------------------------------------------
bool  Ntuple2HepMCFiller::readCurrentEvent() {
        bool filter=false;
        // 1. create an empty event container
        HepMC::GenEvent* event = new HepMC::GenEvent();
        input_ ->setEvent(evtid);
        // 2. fill the evt container - if the read is successful, set the pointer
        if ( this->toGenEvent( evtid, event ) ) evt=event;
        if (evt){ 
                cout <<"| --- Ntuple2HepMCFiller: Event Nr. "  <<evt->event_number() <<" with " <<evt->particles_size()<<" particles --- !" <<endl;             
                //      printHepMcEvent();      
                filter=true;
        }
        if (!evt) {
                cout <<  "Ntuple2HepMCFiller: Got no event :-(" <<endl;
                filter=false;
                delete  evt;  // @@@
        }       
        if (evtid > input_->getNevhep()) filter = false; 
        evtid++;
        return filter;
}

//-------------------------------------------------------------------------
bool Ntuple2HepMCFiller::setEvent(unsigned int event){
  evtid= event;
  return true;
}

//-------------------------------------------------------------------------

bool Ntuple2HepMCFiller::printHepMcEvent() const{
        if (evt!=0) evt->print();       
        return true;
}

//-------------------------------------------------------------------------
HepMC::GenEvent * Ntuple2HepMCFiller::fillCurrentEventData(){
  if (readCurrentEvent())return evt;
  else return NULL;
  
}


//-------------------------------------------------------------------
bool Ntuple2HepMCFiller::toGenEvent( int evtnum, HepMC::GenEvent* evt ){
        // Written according to HepMC /fio/ IO_HEPEVT.cc( hepmc-01-26 tag at Savannah)
        // 1. set event number
        // evt->set_event_number( input_->getNevhep());
        evt->set_event_number( evtnum ) ; 
        // these do not have the correct defaults if hepev4 is not filled
        //evt->set_event_scale( hepev4()->scalelh[1] );
        //evt->set_alphaQCD( hepev4()->alphaqcdlh );
        //evt->set_alphaQED( hepev4()->alphaqedlh ); 
        std::vector<HepMC::GenParticle*> hepevt_particle(input_->getNhep()+1 );
        //2. create a particle instance for each HEPEVT entry and fill a map
        //    create a vector which maps from the HEPEVT particle index to the 
        //    GenParticle address
        //    (+1 in size accounts for hepevt_particle[0] which is unfilled)
        hepevt_particle[0] = 0;
        for ( int i1 = 1; i1 <= input_->getNhep(); ++i1 ) {
                hepevt_particle[i1] = createParticle(i1);
        }       
        std::set<HepMC::GenVertex*> new_vertices; 
        
        // 3.+4. loop over HEPEVT particles AGAIN, this time creating vertices
        for ( int i = 1; i <= input_->getNhep(); ++i ) {
                // We go through and build EITHER the production or decay 
                // vertex for each entry in hepevt
                // Note: since the HEPEVT pointers are bi-directional, it is
                buildProductionVertex( i, hepevt_particle, evt, 1 );
        }
        //      hepevt_particle.clear();
        return true;
}

//-------------------------------------------------------
HepMC::GenParticle* Ntuple2HepMCFiller::createParticle( int index ) {
        
        // Builds a particle object corresponding to index in HEPEVT
        HepMC::GenParticle* p 
        = new HepMC::GenParticle(FourVector( input_->getPhep(index,0), 
        input_->getPhep(index,1), 
        input_->getPhep(index,2), 
        input_->getPhep(index,3)),
        input_->getIdhep(index), 
        input_->getIsthep(index));
        p->setGeneratedMass( input_->getPhep(index, 4) );
        p->suggest_barcode( index );
        return p;
}

//-------------------------------------------------------
void Ntuple2HepMCFiller::buildProductionVertex( int i, 
std::vector<HepMC::GenParticle*>& hepevt_particle, 
HepMC::GenEvent* evt, bool printInconsistencyErrors )
{
        // 
        // for particle in HEPEVT with index i, build a production vertex
        // if appropriate, and add that vertex to the event
        HepMC::GenParticle* p = hepevt_particle[i];
        // a. search to see if a production vertex already exists
        int mother = input_->getJmohep(i,0);
        // this is dangerous.  Copying null pointer and attempting to fill 
        //  information in prod_vtx......
        HepMC::GenVertex* prod_vtx = p->production_vertex();
        //  no production vertex and mother exist -> produce vertex
        while ( !prod_vtx && mother > 0 ) {
                prod_vtx = hepevt_particle[mother]->end_vertex();
                if ( prod_vtx ) prod_vtx->add_particle_out( p );
                // increment mother for next iteration
                if ( ++mother > input_->getJmohep(i,1) ) mother = 0;
        }
        // b. if no suitable production vertex exists - and the particle
        // has atleast one mother or position information to store - 
        // make one@@@
        FourVector prod_pos( input_->getVhep(i,0), input_->getVhep(i,1), 
        input_->getVhep(i,2), input_->getVhep(i,3));
        // orginal:
        if ( !prod_vtx && (number_parents(i)>0  || prod_pos!=FourVector(0,0,0,0) )){
                prod_vtx = new HepMC::GenVertex();      
                prod_vtx->add_particle_out( p );
                evt->add_vertex( prod_vtx );
        }
        // c. if prod_vtx doesn't already have position specified, fill it
        if ( prod_vtx && prod_vtx->position()==FourVector(0,0,0,0) ) {
                prod_vtx->set_position( prod_pos );
        }
        // d. loop over mothers to make sure their end_vertices are
        //     consistent
        mother = input_->getJmohep(i,0);
        while ( prod_vtx && mother > 0 ) {
                if ( !hepevt_particle[mother]->end_vertex() ) {
                        // if end vertex of the mother isn't specified, do it now
                        prod_vtx->add_particle_in( hepevt_particle[mother] );
                } 
                else if (hepevt_particle[mother]->end_vertex() != prod_vtx ) {
                        // problem scenario --- the mother already has a decay
                        // vertex which differs from the daughter's produciton 
                        // vertex. This means there is internal
                        // inconsistency in the HEPEVT event record. Print an
                        // error
                        // Note: we could provide a fix by joining the two 
                        //       vertices with a dummy particle if the problem
                        //       arrises often with any particular generator.
                        if ( printInconsistencyErrors ) std::cerr
                                << "Ntuple2HepMCFiller:: inconsistent mother/daughter "
                        << "information in HEPEVT event " 
                        << input_->getNevhep()
                        << std::endl;
                }
                if ( ++mother > input_->getJmohep(i,1) ) mother = 0;
        }
}


//-------------------------------------------------------
int Ntuple2HepMCFiller::number_children( int index )
{
        int firstchild = input_->getJdahep(index,0);
        return ( firstchild>0 ) ? 
        ( 1+input_->getJdahep(index,1)-firstchild ) : 0;
}

//-------------------------------------------------------       
int Ntuple2HepMCFiller::number_parents( int index ) {

        // Fix for cmkin single particle ntpls 
        if (input_->getNhep()==1) return 1;
        // Fix for cmkindouble particle ntpls
        if (input_->getNhep()==2) return 1;     
        int firstparent = input_->getJmohep(index,0);
        if( firstparent <= 0 ) return 0;
        int secondparent = input_->getJmohep(index,1);
        if( secondparent <= 0 ) return 1;
        return secondparent - firstparent + 1;
        
} 

/*************************************************************************************
NEVHEP = event number
NHEP   = number of entries (particles, partons)
ISTHEP = status code
IDHEP  = PDG identifier
JMOHEP = position of 1st and 2nd  mother
JDAHEP = position of 1st and last daughter
PHEP   = 4-momentum and mass            (single precision in ntuple file)
VHEP   = vertex xyz and production time (single precision in ntuple file)


IRNMCP      = run number
IEVMCP         = event number
WGTMCP         = event weight
XSECN          = cross section equivalent
IFILTER        = filter pattern
NVRMCP         = number of additional variables
VARMCP(NMXMCP) = list of additional variables
Note: VARMCP(1) = PARI(17) (=s_hat) is reserved.

There are the following default values:

IRNMCP = 1
IEVMCP = NEVHEP
WGTMCP = 1.0
XSECN  = IFILTER = 0
NVRMCP = 1
*****************************************************************************************/