The HepMC library provides a container for the Monte Carlo truth information coming from a Monte Carlo tool. It is the C++ implementation of the fixed-size COMMON BLOCK structures used in FORTRAN tools. This information is the input to a detector simulation (e.g. GEANT), and it enables experimentalists to connect theoretical predictions to "hits" in a detector.
One way to understand the HepMC record is to follow an example of building an event from scratch, taken from the "examples" directory of the HepMC distribution.
// // In this example we will place the following event into HepMC "by hand" // // 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
Note that this is not an entirely realistic event, because it does not account for all of the energy and momentum coming from the beam particles. However, it is good enough for our purposes. An important property of particles is their status. The rule is that:
- status==1 particles are those that should be processed directly by a detector simulation.
- status==2 particles are unstable, such as a pi0 that will later decay to gamma gamma
- status==3 particles label the beam and the hard process
This processes represents d u~ -> W- gamma, with W- -> d u~. The outgoing d and u~ are given status==1. In a realistic event generator, these quarks will not have status==1, but will be an input to a hadronization model that turns partons into particles.
Now we build the graph, which will look like:
// // p7 # // p1 / # // \v1__p3 p5---v4 # // \_v3_/ \ # // / \ p8 # // v2__p4 \ # // / p6 # // p2 # // #
This is accomplished by the following C++ code:
// First create the event container, with Signal Process 20, event number 1
//
// Note that the HepLorentzVectors will be automatically converted to
// HepMC::FourVector within GenParticle and GenVertex
GenEvent* evt = new GenEvent( 20, 1 );
// define the units
evt->use_units(HepMC::Units::GEV, HepMC::Units::MM);
//
// create vertex 1 and vertex 2, together with their in-particles
GenVertex* v1 = new GenVertex();
evt->add_vertex( v1 );
v1->add_particle_in( new GenParticle( HepLorentzVector(0,0,7000,7000),
2212, 3 ) );
GenVertex* v2 = new GenVertex();
evt->add_vertex( v2 );
v2->add_particle_in( new GenParticle( HepLorentzVector(0,0,-7000,7000),
2212, 3 ) );
//
// create the outgoing particles of v1 and v2
GenParticle* p3 =
new GenParticle( HepLorentzVector(.750,-1.569,32.191,32.238), 1, 3 );
v1->add_particle_out( p3 );
GenParticle* p4 =
new GenParticle( HepLorentzVector(-3.047,-19.,-54.629,57.920), -2, 3 );
v2->add_particle_out( p4 );
//
// create v3
GenVertex* v3 = new GenVertex();
evt->add_vertex( v3 );
v3->add_particle_in( p3 );
v3->add_particle_in( p4 );
v3->add_particle_out(
new GenParticle( HepLorentzVector(-3.813,0.113,-1.833,4.233 ), 22, 1 )
);
GenParticle* p5 =
new GenParticle( HepLorentzVector(1.517,-20.68,-20.605,85.925), -24,3);
v3->add_particle_out( p5 );
//
// create v4
GenVertex* v4 = new GenVertex(HepLorentzVector(0.12,-0.3,0.05,0.004));
evt->add_vertex( v4 );
v4->add_particle_in( p5 );
v4->add_particle_out(
new GenParticle( HepLorentzVector(-2.445,28.816,6.082,29.552), 1,1 )
);
v4->add_particle_out(
new GenParticle( HepLorentzVector(3.962,-49.498,-26.687,56.373), -2,1 )
);
//
// tell the event which vertex is the signal process vertex
evt->set_signal_process_vertex( v3 );
HepMC has an internal method to print out a representation of the event record. For this event, it looks like this:
GenEvent: #1 ID=20 SignalProcessGenVertex Barcode: -3
Entries this event: 4 vertices, 8 particles.
Beam Particles are not defined.
RndmState(0)=
Wgts(0)=
EventScale -1 [energy] alphaQCD=-1 alphaQED=-1
GenParticle Legend
Barcode PDG ID ( Px, Py, Pz, E ) Stat DecayVtx
________________________________________________________________________________
GenVertex: -1 ID: 0 (X,cT):0
I: 1 10001 2212 +0.00e+00,+0.00e+00,+7.00e+03,+7.00e+03 3 -1
O: 1 10003 1 +7.50e-01,-1.57e+00,+3.22e+01,+3.22e+01 3 -3
GenVertex: -2 ID: 0 (X,cT):0
I: 1 10002 2212 +0.00e+00,+0.00e+00,-7.00e+03,+7.00e+03 3 -2
O: 1 10004 -2 -3.05e+00,-1.90e+01,-5.46e+01,+5.79e+01 3 -3
GenVertex: -3 ID: 0 (X,cT):0
I: 2 10003 1 +7.50e-01,-1.57e+00,+3.22e+01,+3.22e+01 3 -3
10004 -2 -3.05e+00,-1.90e+01,-5.46e+01,+5.79e+01 3 -3
O: 2 10005 22 -3.81e+00,+1.13e-01,-1.83e+00,+4.23e+00 1
10006 -24 +1.52e+00,-2.07e+01,-2.06e+01,+8.59e+01 3 -4
Vertex : -4 ID: 0 (X,cT)=+1.20e-01,-3.00e-01,+5.00e-02,+4.00e-03
I: 1 10006 -24 +1.52e+00,-2.07e+01,-2.06e+01,+8.59e+01 3 -4
O: 2 10007 1 -2.44e+00,+2.88e+01,+6.08e+00,+2.96e+01 1
10008 -2 +3.96e+00,-4.95e+01,-2.67e+01,+5.64e+01 1
________________________________________________________________________________
The record can be traced either through the vertices, or the particles. The most naive way to access the event information is to loop (or iterate) over the particles. Here is an example of how to sum the four momentum and print out the information about stable particles:
HepLorentzVector sum;
for ( GenEvent::particle_const_iterator p = evt->particles_begin();
p != evt->particles_end(); ++p ){
if( (*p)->status() == 1 ) {
sum += convertTo( (*p)->momentum() );
(*p)->print();
}
}
std::cout << "Vector Sum: " << sum << std::endl;
This is also the way to identify all of the stable particles which may be used as the input to a jet algorithm. Of course, you would want to exclude particles such as neutrinos (or neutralinos for a SUSY signal) from such a list.