From Nuclear Physics Group Documentation Pages
Revision as of 19:34, 13 August 2008 by Dan (talk | contribs) (→‎Overview)
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to navigationJump to search


Hall B at JLAB is upgrading their detector to 12GeV. One part of the new detector array is the SVT, or Silicon Vertex Tracker. UNH is helping to design the SVT, and develop the software which will interpret the output from the electronics. For software, GEMC is the main program we will be working on. It is built on GEANT4, which is basically a repository of the last 100 years of particle physics knowledge in easily-wrapped libraries.

Software Locations

We have set-up a suite of software which will be necessary/useful to running GEMC in a shared directory. The top directory is


Under this we have directories with binaries for a handful of platforms. The most interesting one currently is RHEL5_x64, which contains binaries compiled for Red Hat Enterprise Linux 5, AMD64. It is the most interesting because it is the most complete, and it is suitable for use whilst SSHed to Pumpkin, the NPG's hardcore Opteron monster. Note that einstein and mariecurie are also RHEL5_64.

Using GEMC

Speaking of which, do the following:

ssh to Pumpkin, forwarding the X protocol:

ssh -X

Now you will want to set-up your environment. We have an example skeleton script at /net/data/pumpkin2/CLAS12/RHEL5_x64/ As Maurik pointed out, it does not do a lot of things it should. For example, it does not scrub your environment of things which could cause trouble, as the stock GEANT4 environment scrupulously does. In fact, you probably should copy it to your home directory and customize it. One thing you will almost certainly want to do is export a value for the GEMC environment variable before using it. Let's imagine that you don't need to customize it for the moment:

export GEMC=/net/data/pumpkin2/CLAS12/RHEL5_x64/GEMC
source /net/data/pumpkin2/CLAS12/RHEL5_x64/

Great! Things should basically just work now, at least somewhat. Let's try running GEMC:


Now click in the File menu of the resulting G4UI Session and choose Quitter. This should bring up the main GUI of GEMC. Using what settings? Who knows! In a minute we will try specifying a few. First explore the interface. Note that you have some geometry in a viewer, plus a window with controls. Try zooming and panning around, blah blah blah.

Now the useful stuff. Try this:

$GEMC/gemc -USE_QT=0 -GT=BST -N=1000 -PRINT_EVENT=100 -OUTPUT="evio, 1000protons.ev" \
-BEAM_P="proton, 2.0*GeV, 0*deg, 0*deg" -SPREAD_P="1.8*GeV, 180*deg, 360*deg" \
-SPREAD_V="(0.1, 0.1, 2.5)cm" -LUMI_P="e-, 11*GeV, 0, 0" -LUMI_V="(0, 0, -5)cm" \

This will run GEMC with no GUI (-USE_QT=0), using the geometry table BST (-GT=BST), simulating 1000 events (-N=1000), printing the progress every 100 events (-PRINT_EVENT=100), storing the output in EVIO format in the file 1000protons.ev (-OUTPUT="evio, 1000protons.ev"), with a beam of protons of energy 2.0 GeV (-BEAM_P="proton, 2.0*GeV, 0*deg, 0*deg"), then there's something about the vertex I don't understand. Also, some other stuff.

Analysis with ROOT

Right. So we've got this file with data now. What do we do with it? Let's analyze it with ROOT! But first we need to convert it from the EVIO format to that used by ROOT:


Assuming you're still on pumpkin, and have your environment set-up reasonably, you should be able to just run ROOT by typing:


This should start ROOT's interactive interpreter. Start by running

TBrowser B;

Now navigate the directory structure to find your ROOT file, and open it. Try double-clicking on a few trees to display the data in them. To quit ROOT, type


Skipping to ROOT examples: Looking at the relationship between different sets of data is far more useful than just the raw data. Try running this script:

// Daniel Mannarino                         July 14, 2008
// This script reads in rootfiles and makes a plot of phi
// vs. z for protons
// --------------------------------------------------------

// All my standard ROOT includes...
#include <TString.h>
#include <string> 
#include <fstream>
#include <iostream> 
#include <TCanvas.h> 
#include <TH1.h>
#include <TH2.h>
#include <TF1.h>
#include <TH3F.h>
#include <TFile.h>
#include <TTree.h> 
#include <TROOT.h>
#include <TStyle.h>
#include <TApplication.h>
#include <TPaveText.h>
#include <TPaveStats.h>
#include <TPaveLabel.h>
#include <vector>
#include <TProfile.h>
#include <TF1.h>
#include <TGraphErrors.h>
#include <TMultiGraph.h>
#include <TGraph.h>
#include <TDatime.h>
#include <TLegend.h>
#include <TAxis.h>
#include <TGaxis.h>
#include <TPad.h>
#include <TLatex.h>
#include <TRandom.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <TMath.h>
#include <TMinuit.h>
#include <TStopwatch.h>
#include <TSpectrum.h>
#include "TMath.h"

// Interactive includes
#include <TRint.h>
#include <TApplication.h>
#include <TGX11.h>
#include <TVirtualX.h>
#include <TFriendElement.h>

using namespace std ;

// Make tree a global variable so that we can access it outside this function.
// You can then also access it from the command line. 
TTree *tree;

// ************  The main program starts here. ****************//

int main(int argc, char** argv)

  //   Initialize root, get our files and  //
  // make them friends, create our canvas. //
  TApplication theApp("App", &argc, argv);
  // gStyle->SetOptFit(0111);

  //   Proton files

  TFile *pro1 = new TFile("$DATA/proton_500000.root");
  //Get bst tree
  TTree *bst = (TTree*)pro1->Get("bstT");
  bst->AddFriend("genT", "$DATA/proton_500000.root");

  //   Pion files

  TFile *pi1 = new TFile("$DATA/pion_500000.root");
  TTree *bst_pi = (TTree*)pi1->Get("bstT");
  bst_pi->AddFriend("genT", "$DATA/pion_500000.root");
  // Make your canvas

  TCanvas* MyCanvas1 = new TCanvas("MyCanvas1", "Some Graphs", 1);
  TPaveLabel *title1 = NULL;
  title1 = new TPaveLabel(0.1,0.96,0.9,0.99,"A Plot!  Yay!");
  TDatime *now1;
  now1 = new TDatime();
  TPaveLabel *date1 = new TPaveLabel(0.7,0.01,0.9,0.03,now1->AsString());

   TPad* graphPad1 = new TPad("Graphs","Graphs",0.0,0.0,1.0,1.0);
  //   Create and fill our histos          //

  TH2F *Edep_strip = new TH2F("Edep_strip","Edep vs. strip (Protons)",100,0,275,100,0,0.54);

  //   Draw our graphs                     //
 // Draw a legend
  TLegend *leg = new TLegend(0.4,0.6,0.89,0.89);

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

  // Main application eventloop. Calls system dependent eventloop
  // -> runs an interactive interface 


} /**********  End of main *************/

Now look through the code and modify it to display and compare other values from the tree.

For more examples, check out the graphs and code at SVT_Results_at_UNH