---

RMCMacro.cpp

{
/*** This macro gives an example of how to navigate the ReactorTree data
 * structure and make root plots.  What it does is create charge spectra 
 * for neutrons and positrons with radii > 200 cm, and superimposes them on
 * the same plot, and then creates an eps file which can be looked at with
 * your favorite postscript viewer.  
 *
 * Remember you need to load the shared libraray  libRMCEvent.so in order to use 
 * the RMCEvent classes.  To do this, make  sure there is a libRMCEvent.so in your 
 * root directory, and do  .L libRMCEvent.so before you execute this macro 
 * (with .x RMCMacro.cpp). The shared library is created by the ReactorFsim Makefile
 * and placed in the ./src directory.  You  only need to load the shared library 
 * once per root session.
 ***/

//
// First, let's re-set some graphical options, to overcome root's 
// disasterous defaults.
//
     gROOT->SetStyle("Plain");

     gStyle->SetOptStat(0);
     gStyle->SetOptFit(0);
     gStyle->GetAttDate()->SetTextColor(1);
     gStyle->SetOptTitle(0); // no title; comment out if you want a title
     gStyle->SetLabelFont(132,"XYZ");
     gStyle->SetTextFont(132);
     gStyle->SetTitleFont(132,"XYZ");

     gROOT->ForceStyle();

// Now we read in the root file using the TFile class, create the root tree,
// and set the first branch to be an RMCEvent object.  RMCEvent holds an 
// array of secondary objects (for example, neutron and positron) and each
// of these secondary objects will hold an array of PMT hits.

      TFile hfile("Events.root");
      TTree *tree = (TTree*)hfile.Get("RTree");
      RMCEvent *event = new RMCEvent();
      tree->SetBranchAddress("RMCEvent",&event);

// Now create histograms for the neutrons and positrons

      TH1F *pmtqn = new TH1F("pmtqn"," ",50,-0.1,4.0);
      pmtqn->SetLineColor(kBlue);
      pmtqn->SetLineWidth(3);
      pmtqn->SetLineStyle(1);
      pmtqn->SetXTitle("Charge (photoelectrons)");
      pmtqn->SetYTitle("Hits/0.08 pe");
      pmtqn->SetTitle("Neutron and Positron Charge Spectra for R > 200 cm");

      TH1F *pmtqpos = new TH1F("pmtqpos"," ",50,-0.1,4.0);
      pmtqpos->SetLineColor(kRed);
      pmtqpos->SetLineWidth(3);
      pmtqpos->SetLineStyle(2);
      pmtqpos->SetXTitle("Charge (photoelectrons)");
      pmtqpos->SetYTitle("Hits/0.08 pe");
      pmtqpos->SetTitle("Neutron and Positron Charge Spectra for R > 200 cm");

// And a legend to tell us which is which

      TLegend *leg = new TLegend(0.60,0.70,0.88,0.88);
      leg->AddEntry(pmtqpos,"Positrons","L");
      leg->AddEntry(pmtqn,"Neutrons","L");

      Int_t nentries = (Int_t)tree->GetEntries(); // this gets the number of events
                                                 // in the tree

// Now let's define some useful quantities we want to get out of the tree
// structure, like the event number, the number of secondaries, the event
// type, the number of hits, energy, and generated position.
      Int_t num;
      Int_t numsecs;
      Int_t type,hits;
      Double_t energy,rgen;
      Double_t* vertex;

// Now loop over all the events held in the tree
      for (Int_t i=0; i < nentries; i++){
        tree->GetEntry(i);          //This gets the next event
        num = event->GetEventNum();
        numsecs = event->GetNumSecondaries();
        vertex = event->GetInteractionVertex();
        rgen = sqrt(vertex[0]*vertex[0]+vertex[1]*vertex[1]+
                        vertex[0]*vertex[0]);

// Now let's create a NULL pointer of the RMCSecondary class, to use to
// access information about the secondaries
        RMCSecondary *secondary = NULL;
        for (Int_t isec=0; isec<numsecs; isec++){

// Now we actually get the secondary information, which is the isec element
// in the array held in the event 
         secondary = (RMCSecondary *) event.fSecondaries->At(isec);
         type = secondary->GetSecondaryType();
         energy = secondary->GetSecondaryKE();
         hits = secondary->GetNumPMTs();
         int npe;
         double Time;
         Float_t Q;

// Now that we know how many PMT hits there are, let's loop through them and
// put the measured charge in the relevant histograms
         for (Int_t j=0; j<hits; j++){
                RMCPMT *pmt = NULL;
                pmt = (RMCPMT *) secondary.fPMTs->At(j);
                Q = (Float_t) pmt->GetPMTQ();
                npe = pmt->GetPMTnpe();
                if (npe > 0){      //Let's require the tube to be hit by > 0 photons
                 if (rgen > 200){  //And the event to have been created at r>200 cm
                   if (isec == 0){
                    pmtqpos->Fill(Q);  //fill positron histogram
                   }
                   else 
                    pmtqn->Fill(Q);   //fill neutron histogram
                   }
                 }
                 Time = pmt->GetPMTtime();  //in case we want to look at time
          delete pmt;
          }
          delete secondary;
        }
//      tree->Show();   //uncomment this if you'd like to see summary info on tree
      }
          pmtqn->DrawNormalized(); 
          pmtqpos->DrawNormalized("same");   //we'll draw normalized histos, for overlaying


          leg->SetTextFont(132);  // use Times-Roman font for legend
          leg->Draw("same");     // and draw it

          c1->Print("Qnp.eps");   //make a postscript file
          
        delete event;
//      hfile.Close();
//
   return 0;
}

---