Analysis Code

The code here failed to work initially, I'm assuming it's faily old, and probably a bit out of date. Below is my updated version, it works for me, and produces four graphs.

#define h101_cxx

#include "h101.h"

#include 

#include 

#include 

                                                                            

void h101::Loop()

{

//   In a ROOT session, you can do:

//      Root > .L h101.C

//      Root > h101 t

//      Root > t.GetEntry(12); // Fill t data members with entry number 12

//      Root > t.Show();       // Show values of entry 12

//      Root > t.Show(16);     // Read and show values of entry 16

//      Root > t.Loop();       // Loop on all entries

//

                                                                            

//     This is the loop skeleton where:

//    jentry is the global entry number in the chain

//    ientry is the entry number in the current Tree

//  Note that the argument to GetEntry must be:

//    jentry for TChain::GetEntry

//    ientry for TTree::GetEntry and TBranch::GetEntry

//

//       To read only selected branches, Insert statements like:

// METHOD1:

//    fChain->SetBranchStatus("*",0);  // disable all branches

//    fChain->SetBranchStatus("branchname",1);  // activate branchname

// METHOD2: replace line

//    fChain->GetEntry(jentry);       //read all branches

//by  b_branchname->GetEntry(ientry); //read only this branch

                                                                            

Int_t num_events = h101->GetEntries();

Float_t eta, theta, p, pt;

Int_t num_Z = 0;

Int_t num_g = 0;

// Make an empty histogram named "eta"

// (declared on stack memory so it can be used after the fuction is over)

                                                                            

cout << "\t Number of Events = " << num_events << endl;

cout << "\t Nhep = " << this->Nhep << endl;

TH1F *myZHisto1 = new TH1F("Z eta","Z in CMS eta (H0 150)",25, -5, 5);

TH1F *myZHisto2 = new TH1F("Z theta","Z in CMS theta (H0 150)", 10, -3, 3);

TH1F *mygHisto1 = new TH1F("g eta","gluons in CMS eta (H0 150)",25, -5, 5);

TH1F *mygHisto2 = new TH1F("g theta","glouns in CMS theta (H0 150)", 10, -3, 3);

                                                                                                      

   for (Int_t i=0; i< num_events;i++)

   {

      GetEntry(i);  // This sets pointer to a particular Event #

                                                                                                      

      for (Int_t j=0; j < this->Nhep; j++)

      {

                                                                                                      

          if(abs(this->Idhep[j]) == 23) //Sets Particle type (Idhep)

         {

                //cout << "\t Event (" << i << ", " << j << ") has ID = " << this->Idhep[j] << endl;

                //cout << "Px=" << Phep[0][j]<< " Py=" << Phep[1][j]<< " Pz=" << Phep[2][j]<                p = sqrt(this->Phep[0][j]*this->Phep[0][j] + this->Phep[1][j]*this->Phep[1][j] + this->Phep[2][j]*this->Phep[2][j]);

                pt = sqrt(this->Phep[0][j]*this->Phep[0][j] + this->Phep[1][j]*this->Phep[1][j]);

                //cout << "p=" << p << endl << "pt=" << pt << endl;

                theta = asin(pt/p);

                eta = -log(tan(theta/2));

                if(this->Phep[2][j] < 0) eta = -1*eta;  // put sign on theta

                if(this->Phep[2][j] < 0) theta = -1*theta;

                //cout << "====>  eta for " << this->Idhep[j] << "\t is " << eta << endl;

                myZHisto1->Fill(eta);

                myZHisto2->Fill(theta);

                num_Z++;

          } // end if particle type

          if(abs(this->Idhep[j]) == 21) //Sets Particle type (Idhep)

         {

                if (Phep[0][j] != 0){

                        //cout << "\t Event (" << i << ", " << j << ") has ID = " << this->Idhep[j] << endl;

                        //cout << "Px=" << Phep[0][j]<< " Py=" << Phep[1][j]<< " Pz=" << Phep[2][j]<                        p = sqrt(this->Phep[0][j]*this->Phep[0][j] + this->Phep[1][j]*this->Phep[1][j] + this->Phep[2][j]*this->Phep[2][j]);

                        pt = sqrt(this->Phep[0][j]*this->Phep[0][j] + this->Phep[1][j]*this->Phep[1][j]);

                        //cout << "p=" << p << endl << "pt=" << pt << endl;

                        theta = asin(pt/p);

                        eta = -log(tan(theta/2));

                        if(this->Phep[2][j] < 0) eta = -1*eta;  // put sign on theta

                        if(this->Phep[2][j] < 0) theta = -1*theta;

                        //cout << "====>  eta for " << this->Idhep[j] << "\t is " << eta << endl;

                        mygHisto1->Fill(eta);

                        mygHisto2->Fill(theta);

                        num_g++;

                }

          } // end if particle type

                                                                                                      

      } // end loop on particles(j)

   } // end loop on events (i)

                                                                                                      

  TCanvas *can1 = new TCanvas("can1", "Z Eta", 1);

  TCanvas *can2 = new TCanvas("can2", "Z Theta", 1);

  TCanvas *can3 = new TCanvas("can3", "g Eta", 1);

  TCanvas *can4 = new TCanvas("can4", "g Theta", 1);

                                                                                                      

//Plot Z's

  can1->cd();

  myZHisto1->Draw();

  can2->cd();

  myZHisto2->Draw();

  can3->cd();

  mygHisto1->Draw();

  can4->cd();

  mygHisto2->Draw();

                                                                                                      

  cout << "number of Z's = " << num_Z << endl << "number of g's = " << num_g << endl;

                                                                                                      

} // END myclass::Loop