Creation of a ROOT Tree

//*CMZ :  1.00/07 04/04/97  10.13.50  by  Rene Brun
//*-- Author :    Valery Fine(fine@vxcern.cern.ch)   19/01/97

//______________________________________________________________________________
//             A simple example with a ROOT tree
//             =================================
//
//  This program creates :
//    - a ROOT file
//    - a tree
//  Additional arguments can be passed to the program to control the flow
//  of execution. (see comments describing the arguments in the code).
//      Event  nevent comp split fill
//  All arguments are optional: Default is
//      Event  400      1    1     1
//
//  In this example, the tree consists of one single "super branch"
//  The statement ***tree->Branch("event", event, 64000,split);*** below
//  will parse the structure described in Event.h and will make
//  a new branch for each data member of the class if split is set to 1.
//    - 5 branches corresponding to the basic types fNtrack,fNseg,fNvertex
//           ,fFlag and fTemperature.
//    - 3 branches corresponding to the members of the subobject EventHeader.
//    - one branch for each data member of the class Track of TClonesArray.
//    - one branch for the object fH (histogram of class TH1F).
//
//  if split = 0 only one single branch is created and the complete event
//  is serialized in one single buffer.
//  if comp = 0 no compression at all.
//  if comp = 1 event is compressed.
//  if comp = 2 same as 1. In addition branches with floats in the TClonesArray
//                         are also compressed.
//  The 4th argument fill can be set to 0 if one wants to time
//     the percentage of time spent in creating the event structure and
//     not write the event in the file.
//  In this example, one loops over nevent events.
//  The branch "event" is created at the first event.
//  The branch address is set for all other events.
//  For each event, the event header is filled and ntrack tracks
//  are generated and added to the TClonesArray list.
//  For each event the event histogram is saved as well as the list
//  of all tracks.
//
//  The number of events can be given as the first argument to the program.
//  By default 400 events are generated.
//  The compression option can be activated/deactivated via the second argument.
//
//   ---Running/Linking instructions----
//  This program consists of the following files and procedures.
//    - Event.h event class description
//    - Event.C event class implementation
//    - MainEvent.C the main program to demo this class might be used (this file)
//    - EventCint.C  the CINT dictionary for the event and Track classes
//        this file is automatically generated by rootcint (see Makefile),
//        when the class definition in Event.h is modified.
//
//   ---Analyzing the Event.root file with the interactive root
//        example of a simple session
//   Root > TFile f("Event.root")
//   Root > T.Draw("fNtrack")   //histogram the number of tracks per event
//   Root > T.Draw("fPx")       //histogram fPx for all tracks in all events
//   Root > T.Draw("fXfirst:fYfirst","fNtrack>600")
//                              //scatter-plot for x versus y of first point of each track
//   Root > T.Draw("fH.GetRMS()")  //histogram of the RMS of the event histogram
//
//   Look also in the same directory at the following macros:
//     - eventa.C  an example how to read the tree
//     - eventb.C  how to read events conditionally
//
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

#include <stdlib.h>

#include "TROOT.h"
#include "TFile.h"
#include "TRandom.h"
#include "TTree.h"
#include "TBranch.h"
#include "TClonesArray.h"
#include "TStopwatch.h"

#include "Event.h"


//______________________________________________________________________________
main(int argc, char **argv)
{
   TROOT simple("simple","Example of creation of a tree");

   Int_t nevent = 400;     // by default create 400 events
   Int_t comp   = 1;       // by default file is compressed
   Int_t split  = 1;       // by default, split Event in sub branches
   Int_t write  = 1;       // by default the tree is filled
   Int_t hfill  = 0;       // by default histograms are not filled
   Int_t read   = 0;
   Int_t arg4   = 1;
   Int_t arg5   = 600;     //default number of tracks per event

   if (argc > 1)  nevent = atoi(argv[1]);
   if (argc > 2)  comp   = atoi(argv[2]);
   if (argc > 3)  split  = atoi(argv[3]);
   if (argc > 4)  arg4   = atoi(argv[4]);
   if (argc > 5)  arg5   = atoi(argv[5]);
   if (arg4 ==  0) { write = 0; hfill = 0;}
   if (arg4 ==  1) { write = 1; hfill = 0;}
   if (arg4 == 10) { write = 0; hfill = 1;}
   if (arg4 == 11) { write = 1; hfill = 1;}
   if (arg4 == 20) { write = 0; read  = 1;}  //read sequential
   if (arg4 == 25) { write = 0; read  = 2;}  //read random


   TFile *hfile;
   Event *event;
   TTree *tree;
   TBranch *b = 0;

   // Fill event, header and tracks with some random numbers
   //   Create a timer object to benchmark this loop
   TStopwatch timer;
   timer.Start();
   Int_t nb = 0;
   Int_t ev;
   Int_t bufsize;
   Double_t told = 0;
   Double_t tnew = 0;
   Int_t printev = 100;
   if (arg5 < 100) printev = 1000;
   if (arg5 < 10)  printev = 10000;

//         Read case
   if (read) {
      hfile = new TFile("Event.root");
      tree = (TTree*)hfile->Get("T");
      TBranch *branch = tree->GetBranch("event");
      event = new Event();
      branch->SetAddress(&event);
      Int_t nentries = (Int_t)tree->GetEntries();
      nevent = TMath::Max(nevent,nentries);
      if (read == 1) {  //read sequential
         for (ev = 0; ev < nevent; ev++) {
            if (ev%printev == 0) cout<<"event="<<ev<<endl;
            nb += tree->GetEvent(ev);        //read complete event in memory
         }
      } else {    //read random
         Int_t evrandom;
         for (ev = 0; ev < nevent; ev++) {
            if (ev%printev == 0) cout<<"event="<<ev<<endl;
            evrandom = nevent*gRandom->Rndm(1);
            nb += tree->GetEvent(evrandom);  //read complete event in memory
         }
      }
   } else {
//         Write case
      // Create a new ROOT binary machine independent file.
      // Note that this file may contain any kind of ROOT objects, histograms,
      // pictures, graphics objects, detector geometries, tracks, events, etc..
      // This file is now becoming the current directory.
      hfile = new TFile("Event.root","RECREATE","TTree benchmark ROOT file");
      hfile->SetCompressionLevel(comp);

     // Create histogram to show write_time in function of time
     Float_t curtime = -0.5;
     Int_t ntime = nevent/printev;
     TH1F *htime = new TH1F("htime","Real-Time to write versus time",ntime,0,ntime);
     // Create one event
      event = new Event();
      if (hfill) event->Hcreate();

     // Create a ROOT Tree and one superbranch
      TTree *tree = new TTree("T","An example of a ROOT tree");
      tree->SetAutoSave(1000000000);  // autosave when 1 Gbyte written
      bufsize = 256000;
      if (split)  bufsize /= 4;
      TBranch *b = tree->Branch("event", "Event", &event, bufsize,split);

      for (ev = 0; ev < nevent; ev++) {
         if (ev%printev == 0) {
            tnew = timer.RealTime();
            printf("event:%d, rtime=%f s\n",ev,tnew-told);
            htime->Fill(curtime,tnew-told);
            curtime += 1;
            told=tnew;
            timer.Continue();
         }

         Float_t sigmat, sigmas;
         gRandom->Rannor(sigmat,sigmas);
         Int_t ntrack   = Int_t(arg5 +arg5*sigmat/120.);
         Float_t random = gRandom->Rndm(1);

         event->SetHeader(ev, 200, 960312, random);
         event->SetNseg(Int_t(10*ntrack+20*sigmas));
         event->SetNvertex(1);
         event->SetFlag(UInt_t(random+0.5));
         event->SetTemperature(random+20.);

         //  Create and Fill the Track objects
         for (Int_t t = 0; t < ntrack; t++) event->AddTrack(random);

         if (write) nb += tree->Fill();  //fill the tree

         if (hfill) event->Hfill();      //fill histograms

         event->GetTracks()->Clear();
      }
      if (write) { hfile->Write();   tree->Print();}
   }

   //  Stop timer and print results
   timer.Stop();
   Float_t mbytes = 0.000001*nb;
   Double_t rtime = timer.RealTime();
   Double_t ctime = timer.CpuTime();


   printf("\n%d events and %d bytes processed.\n",nevent,nb);
   printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);
   if (read) {
      printf("You read %f Mbytes/Realtime seconds\n",mbytes/rtime);
      printf("You read %f Mbytes/Cputime seconds\n",mbytes/ctime);
   } else {
      printf("compression level=%d, split=%d, arg4=%d\n",comp,split,arg4);
      printf("You write %f Mbytes/Realtime seconds\n",mbytes/rtime);
      printf("You write %f Mbytes/Cputime seconds\n",mbytes/ctime);
      //printf("file compression factor = %f\n",hfile.GetCompressionFactor());
   }
   hfile->Close();
   return 0;
}
//*CMZ :  1.00/00 14/02/97  11.04.07  by  Rene Brun
//*-- Author :    Rene Brun   19/08/96
//______________________________________________________________________________
//                       Event and Track classes
//                       =======================
//
//  The Event class is a naive/simple example of an event structure.
//     public:
//        Int_t          fNtrack;
//        Int_t          fNseg;
//        Int_t          fNvertex;
//        UInt_t         fFlag;
//        Float_t        fTemperature;
//        EventHeader    fEvtHdr;
//        TClonesArray  *fTracks;
//        TH1F          *fH;
//
//   The EventHeader class has 3 data members (integers):
//     public:
//        Int_t          fEvtNum;
//        Int_t          fRun;
//        Int_t          fDate;
//
//
//   The Event data member fTracks is a pointer to a TClonesArray.
//   It is an array of a variable number of tracks per event.
//   Each element of the array is an object of class Track with the members:
//     private:
//        Float_t      fPx;           //X component of the momentum
//        Float_t      fPy;           //Y component of the momentum
//        Float_t      fPz;           //Z component of the momentum
//        Float_t      fRandom;       //A random track quantity
//        Float_t      fMass2;        //The mass square of this particle
//        Float_t      fBx;           //X intercept at the vertex
//        Float_t      fBy;           //Y intercept at the vertex
//        Float_t      fMeanCharge;   //Mean charge deposition of all hits of this track
//        Float_t      fXfirst;       //X coordinate of the first point
//        Float_t      fXlast;        //X coordinate of the last point
//        Float_t      fYfirst;       //Y coordinate of the first point
//        Float_t      fYlast;        //Y coordinate of the last point
//        Float_t      fZfirst;       //Z coordinate of the first point
//        Float_t      fZlast;        //Z coordinate of the last point
//        Float_t      fCharge;       //Charge of this track
//        Int_t        fNpoint;       //Number of points for this track
//        Short_t      fValid;        //Validity criterion
//
//   An example of a batch program to use the Event/Track classes is given
//   in this directory: MainEvent.
//   Look also in the same directory at the following macros:
//     - eventa.C  an example how to read the tree
//     - eventb.C  how to read events conditionally
//
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

#include "TRandom.h"

#include "Event.h"


ClassImp(EventHeader)
ClassImp(Event)
ClassImp(Track)

TClonesArray *gTracks;
TH1F *gHist;
TH1F *hNtrack,*hNseg, *hTemperature;
TH1F *hPx, *hPy, *hPz, *hRandom, *hMass2, *hBx, *hBy;
TH1F *hMeanCharge,*hXfirst,*hXlast,*hYfirst,*hYlast;
TH1F *hZfirst,*hZlast,*hCharge,*hNpoint,*hValid;

//______________________________________________________________________________
Event::Event()
{
//   Create one Event object
//   When the constructor is invoked for the first time, the global
//   variable gTracks is NULL. The TClonesArray gTracks is created.
//   The histogram fH is created.
//
   fNtrack = 0;
   if (!gTracks) gTracks = new TClonesArray("Track", 1000);
   fTracks = gTracks;
   fH = 0;
   }

//______________________________________________________________________________
Event::~Event()
{
   Clear();
}

//______________________________________________________________________________
void Event::AddTrack(Float_t random)
{
// Add a new track to the list of tracks for this event.
// To avoid calling the very time consuming operator new for each track,
// the standard but not well know C++ operator "new with placement" is called.
// if tracks[i] is NULL, a new Track object will be created
// otherwise the previous Track[i] will be overwritten.

   TClonesArray &tracks = *fTracks;
   new(tracks[fNtrack++]) Track(random);
}

//______________________________________________________________________________
void Event::Clear()
{
   fTracks->Clear();
   delete fH;
   fH = 0;
}

//______________________________________________________________________________
void Event::Hcreate()
{

//       Create histograms
   hNtrack      = new TH1F("hNtrack",    "Ntrack",100,575,625);
   hNseg        = new TH1F("hNseg",      "Nseg",100,5800,6200);
   hTemperature = new TH1F("hTemperature","Temperature",100,19.5,20.5);
   hPx          = new TH1F("hPx",        "Px",100,-4,4);
   hPy          = new TH1F("hPy",        "Py",100,-4,4);
   hPz          = new TH1F("hPz",        "Pz",100,0,5);
   hRandom      = new TH1F("hRandom",    "Random",100,0,1000);
   hMass2       = new TH1F("hMass2",     "Mass2",100,0,12);
   hBx          = new TH1F("hBx",        "Bx",100,-0.5,0.5);
   hBy          = new TH1F("hBy",        "By",100,-0.5,0.5);
   hMeanCharge  = new TH1F("hMeanCharge","MeanCharge",100,0,0.01);
   hXfirst      = new TH1F("hXfirst",    "Xfirst",100,-40,40);
   hXlast       = new TH1F("hXlast",     "Xlast",100,-40,40);
   hYfirst      = new TH1F("hYfirst",    "Yfirst",100,-40,40);
   hYlast       = new TH1F("hYlast",     "Ylast",100,-40,40);
   hZfirst      = new TH1F("hZfirst",    "Zfirst",100,0,80);
   hZlast       = new TH1F("hZlast",     "Zlast",100,0,250);
   hCharge      = new TH1F("hCharge",    "Charge",100,-1.5,1.5);
   hNpoint      = new TH1F("hNpoint",    "Npoint",100,50,80);
   hValid       = new TH1F("hValid",     "Valid",100,0,1.2);
}

//______________________________________________________________________________
void Event::Hfill()
{

//       Fill histograms

   hNtrack->Fill(fNtrack);
   hNseg->Fill(fNseg);
   hTemperature->Fill(fTemperature);

   for (Int_t itrack=0;itrack<fNtrack;itrack++) {
      Track *track = (Track*)fTracks->UncheckedAt(itrack);
      hPx->Fill(track->GetPx());
      hPy->Fill(track->GetPy());
      hPz->Fill(track->GetPz());
      hRandom->Fill(track->GetRandom());
      hMass2->Fill(track->GetMass2());
      hBx->Fill(track->GetBx());
      hBy->Fill(track->GetBy());
      hMeanCharge->Fill(track->GetMeanCharge());
      hXfirst->Fill(track->GetXfirst());
      hXlast->Fill(track->GetXlast());
      hYfirst->Fill(track->GetYfirst());
      hYlast->Fill(track->GetYlast());
      hZfirst->Fill(track->GetZfirst());
      hZlast->Fill(track->GetZlast());
      hCharge->Fill(track->GetCharge());
      hNpoint->Fill(track->GetNpoint());
      hValid->Fill(track->GetValid());
   }
}

//______________________________________________________________________________
void Event::SetHeader(Int_t i, Int_t run, Int_t date, Float_t random)
{
   fNtrack = 0;
   fEvtHdr.Set(i, run, date);
   if (!gHist)   gHist   = new TH1F("hstat","Event Histogram",100,0,1);
   fH = gHist;
   fH->Fill(random);
}

//______________________________________________________________________________
Track::Track(Float_t random) :TObject()
{
//*-*-*-*-*-*-*-*-*-*-*Track normal constructor*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
//*-*                  =======================
//   Note that in this example, data members do not have any physical meaning

  Float_t a,b,sigmat,sigmas,px,py,pz;
  gRandom->Rannor(px,py);
  fPx = px;
  fPy = py;
  fPz = TMath::Sqrt(px*px+py*py);
  fRandom = 1000*random;
  if (fRandom < 10) fMass2 = 0.08;
  else if (fRandom < 100) fMass2 = 0.8;
  else if (fRandom < 500) fMass2 = 4.5;
  else if (fRandom < 900) fMass2 = 8.9;
  else  fMass2 = 9.8;
  gRandom->Rannor(a,b);
  fBx = 0.1*a;
  fBy = 0.1*b;
  fMeanCharge = 0.01*gRandom->Rndm(1);
  gRandom->Rannor(a,b);
  fXfirst = a*10;
  fXlast  = b*10;
  gRandom->Rannor(a,b);
  fYfirst = a*12;
  fYlast  = b*16;
  gRandom->Rannor(a,b);
  fZfirst = 50 + 5*a;
  fZlast  = 200 + 10*b;
  fCharge = Float_t(Int_t(3*gRandom->Rndm(1)) - 1);
  fNpoint = Int_t(60+10*gRandom->Rndm(1));
  fValid  = Int_t(0.6+gRandom->Rndm(1));

}