[Feedback][Tutorial Contents][Welcome Page]

Full Reconstruction Example

The examples up to now have used the Monte Carlo event information as a basis for analysis. However, to have a full understanding of the capabilities of a proposed detector, one must model and fully simulate the response of a detector design in the presence of all the physics processes. The LCD simulations accomplish this using a program called GISMO, which is an object-oriented program written in C++ offering similar functionality to GEANT. Hits in the sensitive tracking volumes and calorimeters are written out to files with the .sio extension. This tutorial presents an example of a simple reconstruction and analysis of these fully simulated events.

First we show the entire routine, and then follow it with an explanation of different sections.

//
// This is an example reconstruction and simple analysis
// program for LCD events.
//
import java.util.*;

import hep.lcd.util.driver.Driver;
import hep.lcd.util.driver.AbstractProcessor;

import hep.lcd.event.LCDEvent;

import hep.lcd.mc.smear.SmearDriver;
import hep.lcd.recon.tracking.TrackReco;
import hep.lcd.event.Track;
import hep.lcd.event.TrackList;

import hep.lcd.recon.cluster.simple.SimpleClusterBuilder;
import hep.lcd.event.Cluster;
import hep.lcd.event.ClusterList;

import hep.physics.BasicHepLorentzVector;
import hep.physics.jets.AbstractJetFinder;
import hep.physics.jets.JadeEJetFinder;
import hep.physics.jets.DurhamJetFinder;

import hep.lcd.geometry.Detector;

final public class ExampleReconstruction extends Driver
{
   public ExampleReconstruction()
   {
 
      // Smear Tracker hits with default resolutions
      add(new SmearDriver());
      // Find tracks in combined inner tracker
      add(new TrackReco()); 
      // Find Clusters using simple clustering algorithm
      add(new SimpleClusterBuilder());
      // Find jets using tracks
      add(new TrackJetFinder());
      // Find jets using clusters
      add(new ClusterJetFinder());
 
   }
}

//Find Jets using reconstructed Tracks
class TrackJetFinder extends AbstractProcessor
{
   // The jet finder
   AbstractJetFinder _jetFinder;

   // Default Constructor
   public TrackJetFinder()
   {
      _jetFinder = new JadeEJetFinder(0.02);
   }

   // Fully qualified Constructor
   public TrackJetFinder( AbstractJetFinder jf)
   {
      _jetFinder = jf;
   }   

   public void process(final LCDEvent event)
   {
      // We will use the charged tracks in this example to find jets
      // Get the list of tracks from the event
      TrackList tracks = event.getTrackList();
      //feed them to the jet finder
      _jetFinder.setEvent(tracks.getTracks());
 
      // Fill some diagnostic histograms 
      histogram("Tracks Njets").fill(_jetFinder.njets());
 
      // Loop over the found jets
      double totalJetEnergy = 0;
      double totalParticleEnergy = 0;
      for (int ijet =0; ijet<_jetFinder.njets(); ijet++)
      {
         totalJetEnergy += _jetFinder.jet(ijet).t();
         histogram("Tracks Jet Energy").fill(_jetFinder.jet(ijet).t());
         histogram("Number of Tracks in Jet").fill(_jetFinder.nParticlesPerJet(ijet));
         // loop over all the tracks in the jet, and add up their energy
         // (should be the same as the jet energy)
         double etot = 0;
         Enumeration e = _jetFinder.particlesInJet(ijet);
         while (e.hasMoreElements()) 
            etot += ((Track) e.nextElement()).mag();
         histogram("Tracks Ptot").fill(etot);
         totalParticleEnergy += etot;
      }
      histogram("Total Track Jet Energy").fill(totalJetEnergy);
      histogram("Total Track Particle Momentum").fill(totalParticleEnergy);
   }
}

// Find Jets using Calorimeter Clusters
class ClusterJetFinder extends AbstractProcessor
{
   // the jet finder
   AbstractJetFinder _jetFinder;

   // Default Constructor
   public ClusterJetFinder()
   {
      _jetFinder = new DurhamJetFinder(0.02);
   }

   // Fully qualified Constructor
   public ClusterJetFinder( AbstractJetFinder jf)
   {
      _jetFinder = jf;
   }

   public void process(final LCDEvent event)
   {   
      // We will use the calorimeter clusters in this example to find jets      
      Vector finals = new Vector();
      // Get the EM clusters from the event...
      ClusterList em = event.getEMClusters();
      Enumeration eem = em.getClusters();
      // Get the hadronic clusters from the event...
      ClusterList had = event.getHADClusters();
      Enumeration ehad = had.getClusters();

      // which detector are we using? Needed to apply sampling fractions.
      // This will be automatic in the future.
      Detector det = event.getCurrentDetector();
      // defaults...
      double emscale = 1/.0697; //cng
      double hadscale = 1/.0189; //cng      
      if(det.toString().equals("ldmar01"))
      {
         emscale = 1/.0263; 
         hadscale = 1/.028;
      }
      else if(det.toString().equals("sdmar01"))
      {
         emscale = 1/.0222;  
         hadscale = 1/.0555;     
      }
      else if(det.toString().equals("pdmar01"))
      {
         emscale = 1/.0697;  
         hadscale = 1/.0189;     
      }

 
      double cEn = 0;
      double pt = 0;
      double px = 0;
      double py = 0;
      double pz = 0;
 
      // loop over the EM clusters
      while (eem.hasMoreElements())
      {
         Cluster c = (Cluster) eem.nextElement();
         cEn = emscale*c.getEnergy();
         histogram("Total Cluster Energy").fill(cEn);
         histogram("EM Cluster Energy").fill(cEn);
         histogram("EM Theta").fill(c.getEnergyTheta());
         histogram("EM CosTheta").fill(Math.cos(c.getEnergyTheta()));
         pt = cEn;
         px = cEn*Math.sin(c.getEnergyPhi())*Math.sin(c.getEnergyTheta());
         py = cEn*Math.cos(c.getEnergyPhi())*Math.sin(c.getEnergyTheta());
         pz = cEn*Math.cos(c.getEnergyTheta());
         // cut out low energy clusters...
         if ( cEn > 0.5 )
         {
            finals.addElement(new BasicHepLorentzVector(pt,px,py,pz));   
         }
      }
      // Loop over the hadronic clusters
      while (ehad.hasMoreElements())
      {
         Cluster c = (Cluster) ehad.nextElement();
         // scale the energy
         cEn = hadscale*c.getEnergy();
         histogram("HAD Cluster Energy").fill(cEn);
         histogram("Total Cluster Energy").fill(cEn);
         histogram("HAD Theta").fill(c.getEnergyTheta());
         histogram("HAD CosTheta").fill(Math.cos(c.getEnergyTheta()));
         pt = cEn;
         px = cEn*Math.sin(c.getEnergyPhi())*Math.sin(c.getEnergyTheta());
         py = cEn*Math.cos(c.getEnergyPhi())*Math.sin(c.getEnergyTheta());
         pz = cEn*Math.cos(c.getEnergyTheta());
         // cut out low energy clusters...
         if ( cEn > 0.5 )
         {
            finals.addElement(new BasicHepLorentzVector(pt + 0.137,px,py,pz));   
         }
      }  
      // Send the list of four-vectors to the jet finder... 
      _jetFinder.setEvent(finals.elements());
      int njets = _jetFinder.njets();
      histogram("Durham clusters Njets").fill(njets);
      double totalJetEnergy = 0;
      double totalParticleEnergy = 0;
      for (int ijet =0; ijet<_jetFinder.njets(); ijet++)
      {
         totalJetEnergy += _jetFinder.jet(ijet).t();
         histogram("Durham Cluster Jet Energy").fill(_jetFinder.jet(ijet).t());
         histogram("Durham Number of Clusters in Jet").fill(_jetFinder.nParticlesPerJet(ijet));
      }
      histogram("Durham Total Jet Energy").fill(totalJetEnergy);     
   }
}

Track Reconstruction

To increase the flexibility of the simulations, the hits as recorded by the sensitive volumes in the full simulation are not digitized. Therefore the first step is to smear these "ideal" hits by the detector resolutions. We will use the default settings and simply add a SmearDriver.

Having smeared the hits, we will now find and fit the tracks in the central tracking chambers. To do so, we add a TrackReco Driver. This performs pattern recognition in the sub-detectors, links the track segments and performs a global fit to the hits. The output tracks are then added back to the event.

      // Smear Tracker hits with default resolutions
      add(new SmearDriver());
      // Find tracks in combined inner tracker
      add(new TrackReco());

The tracks will now show up in the LCD Event Display as red tracks.

The Wired event display is also set up to display the reconstructed tracks.

The individual tracks can be toggled on or off by clicking on the boxes. An example of some of the more advanced features of the display can be seen here.

Calorimeter Cluster Reconstruction

We will cluster the calorimeter cells using a simple clustering algorithm by adding a SimpleClusterBuilder Driver. The resultant clusters found in the electromagnetic and hadronic calorimeters will be added to the event.

      // Find Clusters using simple clustering algorithm
      add(new SimpleClusterBuilder());

That's it for the reconstruction; we've found and fit tracks and clustered the energy depositions in the calorimeter cells.

Example Analyses

In the previous tutorial the Monte Carlo particles were used as inputs to the jet finding algorithms. Here we will provide examples for using the reconstructed tracks and the calorimeter clusters as input. (The use of Energy Flow objects, where tracks are used for charged particles and unassociated calorimeter clusters are used for neutrals, is still under development.)

Using the Jet Finder with Tracks

We construct an AbstractProcessor TrackJetFinder to use the reconstructed tracks for jet finding. The default constructor initializes the member jet finder to be a JadeEJetFinder with a ycut of 0.02 just as an example. The user can pass a different jet finder to the fully qualified constructor. The list of tracks is obtained from the LCDEvent via its getTrackList() method. The resulting TrackList is then sent to the jet finder.

      // Get the list of tracks from the event
      TrackList tracks = event.getTrackList();
      //feed them to the jet finder
      _jetFinder.setEvent(tracks.getTracks());

Analysis of the found jets proceeds exactly as in the previous tutorial.

Using the Jet Finder with Calorimeter Clusters

We construct an AbstractProcessor ClusterJetFinder to use the reconstructed calorimeter clusters for jet finding. Using the SimpleClusterBuilder results in separate clusters being reconstructed in the EM and hadronic calorimeters. The lists of clusters are obtained from the LCDEvent via its getEMClusters() and getHADClusters() methods. We then extract the Enumeration over which we will loop to access the Clusters.

      ClusterList em = event.getEMClusters();
      Enumeration eem = em.getClusters();
      // Get the hadronic clusters from the event...
      ClusterList had = event.getHADClusters();
      Enumeration ehad = had.getClusters();

As was the case for the tracking hits, the calorimeter deposited energies are written out as measured in the full detector simulation. We need, therefore, to convert the sampled energy into deposited energy. To do so, we need to know which detector is being simulated and then set the EM and hadronic sampling fractions appropriately.

     // Which detector are we using? 
     Detector det = event.getCurrentDetector();

The LCDEvent contains the information about the detector which was used to simulate the response. Here we only need to know which detector was used, so we compare against our known detectors. (Of course, in the future this will be done automatically).

We then loop over the Enumeration of clusters and histogram some example quantities and fill a vector of BasicHepLorentzVector objects.

// Creating the vector of final clusters to be sent to the jet finder
Vector finals = new Vector();

// EM
finals.addElement(new BasicHepLorentzVector(pt,px,py,pz));   
// Hadronic
finals.addElement(new BasicHepLorentzVector(pt + 0.137,px,py,pz));

EM objects are considered to be massless, whereas we add the pion rest mass to the hadronic clusters when creating the four-vectors. The resulting list is then sent to the jet finder.

      // Send the list of four-vectors to the jet finder... 
      _jetFinder.setEvent(finals.elements());

Analysis of the found jets proceeds exactly as in the previous tutorial.