Project

General

Profile

Reading ghep.root files

A "genie" executable

When GENIE as a UPS product is setup, an genie executable is made available to $PATH that is really a GENIE-enhanced version of ROOT. It does this by adding most of

Interactive

Browse the resulting file using ROOT/genie

  $ genie gntp.42.ghep.root
  TBrowser tb; // double click on 'genv' and 'gconfig' object
  .q

Start a new interactive 'genie' session to explore how to drill down into
the stored GHEP events:

  $ genie
  genie::PDGLibrary::Instance();  // this initialize the TDatabasePDG
                                  // with GENIE extensions; doing it early
                                  // means that output won't clutter results
  TFile* myfile = TFle::Open("gntp.42.ghep.root");
  TTree* mytree = 0;
  myfile->GetObject("gtree",mytree);

  Long64_t nEntries = mytree->GetEntries();
  cout << nEntries << endl;

  genie::NtpMCEventRecord* myentry = new genie::NtpMCEventRecord();
  mytree->SetBranchAddress("gmcrec",&myentry); // gmcrec is branch name

  mytree->GetEntry(0);  // retrieve an event

  myentry->PrintToStream(cout);  // formatted printout of the event

  // print one particle in the event only
  myentry->event->Particle(0)->Print(); cout << endl;

  myentry->event->XSec()  // leave off ";" to have ROOT print result value

  // vertex info is boring (x,y,z,t)=(0,0,0,0) because gevgen used
  // the PointGeometryAnalyzer and not a real geometry
  // nor did it use a flux that varied in space
  myentry->event->Vertex()->Print();

A script

// a sample script for looping over a GENIE GHEP ntuple such as
// generated by standard GENIE 'gevgen' application
// use:   $ genie plot_gnpt.C
//
// R. Hatcher <rhatcher@fnal.gov> 2014-05-07

void plot_gnpt(string fname="gntp.42.ghep.root", int nprint=10) {
//
// count the number of final state pi0 in each event (print nprint of them)
// make a histogram of neutrino energies of the events

  // trigger this early so that it's screen output isn't disruptive
  genie::PDGLibrary::Instance();

  // open the file, get the TTree
  TFile* myfile = TFile::Open(fname.c_str());
  TTree* mytree = 0;
  myfile->GetObject("gtree",mytree);
  Long64_t nEntries = mytree->GetEntries();
  cout << nEntries << endl;

  // book the histogram
  int nbins = (nEntries>1000) ? 200 : 40;
  TH1D* histe  = new TH1D("enu","E_{#nu};E_{#nu} (GeV)",nbins,0,10.);
  TH2D* histxy = new TH2D("vtx","vertices ; m; m",40,-20.,20.,40,-20.,20.);

  // attach the branch object so we can retrieve each event
  genie::NtpMCEventRecord* myentry = new genie::NtpMCEventRecord();
  mytree->SetBranchAddress("gmcrec",&myentry); // gmcrec is branch name

  for (int i=0; i<nEntries; ++i) {
    mytree->GetEntry(i);

    // info about the probe (i.e. incoming neutrino) & vertex
    double enu = myentry->event->Probe()->P4()->E();
    histe->Fill(enu);
    histxy->Fill(myentry->event->Vertex()->X(),
                 myentry->event->Vertex()->Y());

    int nfinalpi0 = 0;
    // loop over particles in the record for final state pi0
    int nparticles = myentry->event->GetEntries();
    for (int ipart=0; ipart<nparticles; ++ipart ) {
      genie::GHepParticle* apart = myentry->event->Particle(ipart);
      if ( apart->Status() == 1 && apart->Pdg() == 111) ++nfinalpi0;
    } // finished w/ loop over particle in event

    if ( nprint > 0 ) {
      nprint--;
      cout << "i=" << setw(3) << i << " nparticles=" << setw(3)
           << nparticles << " nfinal pi0=" << setw(3) << nfinalpi0 << endl;
    }
  }

  // done with all events, draw the histogram
  TCanvas* c1 = new TCanvas("E#nu","E#nu",600,600);
  c1->cd();
  histe->Draw();
  TCanvas* c2 = new TCanvas("vertices","vertices",600,600);
  c2->cd();
  histxy->Draw("COLZ");
}

Examine the 'plot_gntp.C' script to see how it extends the concepts above
to count the number of pi0's in each event and make a plot of the E_nu
distribution. Run it using the command:

  $ genie 'plot_gnpt.C("gntp.42.ghep.root")'

Use the single quotes (') to hide the ()'s from the shell; use the double
quote to quote the string passed as an argument to the script. This allows
one to use the same script for GHEP ntuple files with different names.

Alternatively one can use "\" to hide special characters from the shell

  $ genie plot_gnpt.C\(\"gntp.42.ghep.root\"\)