#include <string>
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <TROOT.h>
#include <TTree.h>
#include <TFile.h>  
#include <TLeaf.h> 
#include <TH1.h> 
#include <TGlobal.h>
#include <TCanvas.h>

void rates()
{

/*
Rates marco

Purpose: To calculate the rates of the L1 and L2 triggers on a partcular ROOT file outputed from d0trigsim

By: Josh Kalk

Version: 00-01-00 25-JUN-01

Instructions:  First run D0trigsim and create a ROOTtuple to analyze.  Then open ROOT and open the file that you want to analyze and compile and run the macro:

.x rates.C++  

Now type in the name of the ROOTtuple and the luminosity in 1/(cm^2 s)

Some common luminosities are:

5E30  for 0.14 minbias
25E30 for 0.7 minbias
50E30 for 1.4 minbias

or lum = 250/7 min bias

If you are running d0trigsim p09.00.00 or above you can turn on Level 2 by uncommenting the section marked Level 2 below

Note: The final trigger is the combined rate for a Level 1 accept


*/

  // Read in the filename or the ROOTtuple
  string fname;
  cout<<"Please type in the name of the file :";
  cin>>fname;

  TFile *myfile = new TFile(fname.c_str(),"READ");
  TTree *Global=(TTree *)myfile->Get("Global");

  // Read in the Luminosity
  float lum;
  cout<<"Please input the luminosity in 1/(cm^2 s) (e.g. 5E30): ";
  cin>>lum;

  // Number of specific trigger to be analyzed
  const int num_st = 14;

  // Create some histograms to hold the rates
  TH1F *L1rate = new TH1F("L1Rates","L1Rates",num_st+2,0,num_st+1);
  TH1F *L2rate = new TH1F("L2Rates","L2Rates",num_st+2,0,num_st+1);
  
  // Set axis labels
  L1rate->GetXaxis()->SetTitle("L1 Specific Triggers");
  L1rate->GetYaxis()->SetTitle("Rate (Hz)");
  L2rate->GetXaxis()->SetTitle("L2 Specific Triggers");
  L2rate->GetYaxis()->SetTitle("Rate (Hz)");

  //Turn off the stat box
  L1rate->SetStats(kFALSE);
  L2rate->SetStats(kFALSE);

  // Grab the number of events from the ROOTtuple
  float NumEntries = Global->GetEntries();
  // const Int_t constNumEnt = NumEntries;
  cout << "Number of Events: "<< NumEntries << endl;

  // Create array to store the number of events passed for each specific trigger
  int L1st[num_st];
  int L2st[num_st];

  // initialize to zero
  for (Int_t i=0; i < num_st; i++) 
    {
      L1st[i] = 0;
      L2st[i] = 0;
    }

  // Create total effiency for L1 and L2
  int L1tot = 0;
  int L2tot = 0;

  // Create variable to store the cross section
  float weight = 0;

  //Loop over every event and grab all the info
  for (Int_t i=0; i < NumEntries; i++) 
    {
      Global->GetEntry(i);

      // get the weight for the event
      weight += ((Global->GetBranch("EVENT"))->GetLeaf("weight"))->GetValue();

      // set temp to zero.  to be filled with 1 if event passed L1 or L2
      int l1temp = 0;
      int l2temp = 0;

      for (int j=0; j < num_st; j++) 
	{
	  // Define the name for the specific trigger
	  char st[4] = "ST";
	  st[2] = j%10+48;
	  if(j/10 != 0)
	    {
	      st[2] = j/10+48;
	      st[3] = j%10+48;	      
	    }
	  // If L1 ST fired then set temp to 1 and increment L1st
	  if(((Global->GetBranch("L1frm_Specific_Trigger"))
	      ->GetLeaf(st))->GetValue()==1)
	    {
	      l1temp = 1;
	      L1st[j]++;
	    }
	  //L2 not yet available in p08.12.00 if you want to run uncomment
	  /*

	  // Define the name for the specific trigger
	  char bit[5] = "bit";
	  bit[3] = j%10+48;
	  if(j/10 != 0)
	    {
	      bit[3] = j/10+48;
	      bit[4] = j%10+48;	      
	    }
	  // If L2 ST fired then set temp to 1 and increment L1st
	  if(((Global->GetBranch("L2gblHead"))
	      ->GetLeaf(bit))->GetValue()==1)
	    {
	      l1temp = 1;
	      L1st[j]++;
	    }
	  */
	}

      // If l1 fired increment total 
      if(l1temp == 1)
	L1tot++;
      // If l2 fired increment total 
      if(l2temp == 1)
	L2tot++;

    }//end loop over entries

  // Calculate the Rate

  float rate = lum * weight / NumEntries * 1E-27;

  //Fill histograms with all info

  cout<<"ST               L1Rate                 L2Rate"<<endl;

  for (Int_t i=0; i < num_st; i++) 
    {
      cout<<i+1<<"               "<<rate * L1st[i]/NumEntries<<
	"                 "<<L2st[i]/NumEntries<<endl;
      L1rate->Fill(i,rate * L1st[i]/NumEntries);
      L2rate->Fill(i,rate * L2st[i]/NumEntries);
    }
  cout<<"Tot"<<"              "<<rate * L1tot/NumEntries<<
	"                 "<<L2tot/NumEntries<<endl;
  L1rate->Fill(num_st,rate * L1tot/NumEntries);
  L2rate->Fill(num_st,rate * L2tot/NumEntries);

  // Draw all the histograms
  TCanvas *MyC = new TCanvas("MyC","Combined Rates",1);
  MyC->Divide(2);
  MyC->cd(1);
  L1rate->Draw();
  MyC->cd(2);
  L2rate->Draw();

}//end l2gbl()