---

HepMCEventWriter.cxx

Go to the documentation of this file.

// ================================================
// HepMCEventWriter class Implementation
// Puts a HepMC event into the TES. For Atlfast Tests
// ================================================
//

#include "AtlfastAlgs/HepMCEventWriter.h"
#include "AtlfastEvent/MsgStreamDefs.h"
#include "AtlfastEvent/MCparticleCollection.h"
#include "AtlfastUtils/HeaderPrinter.h"

#include <cmath>
// Gaudi includes
#include "GaudiKernel/DataSvc.h"
#include "StoreGate/DataHandle.h"

// CLHEP,HepMC
#include "GeneratorObjects/McEvent.h"
#include "GeneratorObjects/McEventCollection.h"
#include "HepMC/GenEvent.h"
#include "HepMC/ParticleDataTable.h"

//

namespace Atlfast {

  HepMCEventWriter::HepMCEventWriter 
  ( const std::string& name, ISvcLocator* pSvcLocator )
    : Algorithm( name, pSvcLocator ){}
 
  HepMCEventWriter::~HepMCEventWriter() {
    MsgStream log( messageService(), name() ) ;
    log << MSG::INFO << "destructor called" << endreq;
  }
  
  StatusCode HepMCEventWriter::initialize()
  {
    MsgStream log( messageService(), name() ) ;
    log<<"Initialise...."<<endreq;
    HeaderPrinter hp("Atlfast HepMCEventWriter:", log);
    hp.print();

    m_tesIO = new TesIO();

    return StatusCode::SUCCESS;
  }

  StatusCode HepMCEventWriter::finalize(){
    
    MsgStream log( messageService(), name() ) ;
    
    log << MSG::INFO << "Finalizing" << endreq;  
    
    return StatusCode::SUCCESS ;
  }
  
  StatusCode HepMCEventWriter::execute( ){
    MsgStream log( messageService(), name() ) ;
    log << MSG::DEBUG << "Execute() " << endreq;

    McEventCollection* mcCollection = new McEventCollection;
    McEvent* mcEvent = new McEvent("Test",20,1);
    HepMC::GenEvent* evt = mcEvent->pGenEvt();
    this->makeGenEvent(evt);

    mcCollection->push_back(mcEvent);
    m_tesIO->store(mcCollection);
    return StatusCode::SUCCESS;
  }
    
  StatusCode HepMCEventWriter::makeGenEvent(HepMC::GenEvent* evt){
    
    MsgStream log( messageService(), name() ) ;
    log << MSG::DEBUG << "MakeGenEvent() " << endreq;

    //
    //     name status pdg_id  parent Px       Py    Pz       Energy      Mass
    //  1  !p+!    3   2212    0,0    0.000    0.000 7000.000 7000.000    0.938
    //  2  !p+!    3   2212    0,0    0.000    0.000-7000.000 7000.000    0.938
    //=========================================================================
    //  3  !d!     3      1    1,1    0.750   -1.569   32.191   32.238    0.000
    //  4  !u~!    3     -2    2,2   -3.047  -19.000  -54.629   57.920    0.000
    //  5  !W-!    3    -24    1,2    1.517   -20.68  -20.605   85.925   80.799
    //  6  !gamma! 1     22    1,2   -3.813    0.113   -1.833    4.233    0.000
    //  7  !d!     1      1    5,5   -2.445   28.816    6.082   29.552    0.010
    //  8  !u~!    1     -2    5,5    3.962  -49.498  -26.687   56.373    0.006

    // first we construct a ParticleDataTable with all the particles we need
    HepMC::ParticleDataTable pdt("my particle data table");
    // create a particle data entry for the proton and add it to pdt at the
    // same time
    pdt.insert( new HepMC::ParticleData( "p+", 2212,   +1, 0.938,  -1, .5 ) );
    pdt.insert( new HepMC::ParticleData( "d",  1,  -2./3., 0,      -1, .5 ) );
    pdt.insert( new HepMC::ParticleData( "u~", -2, -1./3., 0,      -1, .5 ) );
    pdt.insert( new HepMC::ParticleData( "W-", -24,    -1, 80.396,
                                  HepMC::clifetime_from_width(2.06), 1 )    );
    pdt.insert( new HepMC::ParticleData( "gamma", 22,   0, 0,      -1, 1  ) );

    // print out the GenParticle Data to the screen
    pdt.print();

    // now we build the graph, which will look like
    //                       p7
    // p1                   /
    //   \v1__p3      p5---v4
    //         \_v3_/       \ 
    //         /    \        p8
    //    v2__p4     \ 
    //   /            p6
    // p2
    //

    // First create the event container, with Signal Process 20, event number 1
    // create vertex 1 and vertex 2, together with their inparticles
    HepMC::GenVertex* v1 = new HepMC::GenVertex();
    evt->add_vertex( v1 );
    v1->
      add_particle_in( new HepMC::GenParticle( HepLorentzVector(0,0,7000,7000),
                                               2212, 
                                               3 ) 
                       );
    HepMC::GenVertex* v2 = new HepMC::GenVertex();
    evt->add_vertex( v2 );
    v2->
      add_particle_in( 
                      new HepMC::GenParticle( HepLorentzVector(0,0,-7000,7000),
                                              2212, 
                                              3 ) 
                      );
    //
    // create the outgoing particles of v1 and v2
    HepMC::GenParticle* p3 = 
        new HepMC::GenParticle( 
                               HepLorentzVector(.750,-1.569,32.191,32.238), 
                               1, 
                               3 );
    v1->add_particle_out( p3 );
    HepMC::GenParticle* p4 = 
        new HepMC::GenParticle( HepLorentzVector(-3.047,-19.,-54.629,57.920),
                                -2, 
                                3 );
    v2->add_particle_out( p4 );
    //
    // create v3
    HepMC::GenVertex* v3 = new HepMC::GenVertex();
    evt->add_vertex( v3 );
    v3->add_particle_in( p3 );
    v3->add_particle_in( p4 );
    v3->add_particle_out( 
        new HepMC::GenParticle( HepLorentzVector(-3.813,0.113,-1.833,4.233 ), 
                                22, 
                                1 )
        );
    HepMC::GenParticle* p5 = 
        new HepMC::GenParticle( HepLorentzVector(1.517,-20.68,-20.605,85.925), 
                                -24,
                                3);
    v3->add_particle_out( p5 );
    //
    // create v4
    HepMC::GenVertex* v4 = new HepMC::GenVertex();
    evt->add_vertex( v4 );
    v4->add_particle_in( p5 );
    v4->add_particle_out( 
        new HepMC::GenParticle( HepLorentzVector(-2.445,28.816,6.082,29.552), 
                                1,
                                1 )
        );
    v4->
      add_particle_out( new HepMC::GenParticle( 
                                               HepLorentzVector(3.962,
                                                                -49.498,
                                                                -26.687,
                                                                56.373),
                                               -2,
                                               1 )
                        );
    int i, j, k;
    double delta = 2.5;
    double px, py, pz, e;
    int maxstep = 21;
    for(i=1;i<maxstep;++i){
      px = i*delta;
      for(j=1;j<maxstep;++j){
        py = j*delta;
        for(k=1;k<maxstep;++k){
          pz = k*delta;
          e=sqrt(px*px+py*py+pz*pz);
          log << MSG::INFO << "px,pypz,e" 
              <<px<<" "<<py<<" "<<pz<<" "<<e<<" "<<endreq;
          HepMC::GenParticle* temp = 
            new HepMC::GenParticle( HepLorentzVector(px,py,pz,e), 13, 1);
          v4->add_particle_out(temp);
          }
        }
      }
    //    
    // tell the event which vertex is the signal process vertex
    evt->set_signal_process_vertex( v3 );
    // the event is complete, we now print it out to the screen
    evt->print();
    
    // now clean-up by deleteing all objects from memory
    //
    // deleting the event deletes all contained vertices, and all particles
    // contained in those vertices
    //    delete evt;

    // delete all particle data objects in the particle data table pdt
    pdt.delete_all();
    return StatusCode::SUCCESS;

  }
} // Atlfast

---