HowToReadTree()
{

  // declare  tree
  TChain fChain("tree");

  // chain all input files
  fChain.Add("pythia.H.gamgam.conf.root");


  // set address to your variables

  // Declaration of leaf types
   Int_t           mc_n;
   std::vector<double>  *mc_px;
   std::vector<double>  *mc_py;
   std::vector<double>  *mc_pz;
   std::vector<double>  *mc_E;
   std::vector<float>   *mc_pt;
   std::vector<float>   *mc_m;
   std::vector<float>   *mc_eta;
   std::vector<float>   *mc_phi;
   std::vector<int>     *mc_status;
   std::vector<int>     *mc_barcode;
   std::vector<int>     *mc_pdgId;

  // List of branches
   TBranch        *b_mc_n;   //!
   TBranch        *b_mc_px;   //!
   TBranch        *b_mc_py;   //!
   TBranch        *b_mc_pz;   //!
   TBranch        *b_mc_E;   //!
   TBranch        *b_mc_pt;   //!
   TBranch        *b_mc_m;   //!
   TBranch        *b_mc_eta;   //!
   TBranch        *b_mc_phi;   //!
   TBranch        *b_mc_status;   //!
   TBranch        *b_mc_barcode;   //!
   TBranch        *b_mc_pdgId;   //!


   // Set object pointer
   mc_px = 0;
   mc_py = 0;
   mc_pz = 0;
   mc_E = 0;
   mc_pt = 0;
   mc_m = 0;
   mc_eta = 0;
   mc_phi = 0;
   mc_pdgId = 0;
   mc_status = 0;
   mc_barcode = 0;

   fChain.SetBranchAddress("mc_n", &mc_n, &b_mc_n);
   fChain.SetBranchAddress("mc_px", &mc_px, &b_mc_px);
   fChain.SetBranchAddress("mc_py", &mc_py, &b_mc_py);
   fChain.SetBranchAddress("mc_pz", &mc_pz, &b_mc_pz);
   fChain.SetBranchAddress("mc_E", &mc_E, &b_mc_E);
   fChain.SetBranchAddress("mc_pt", &mc_pt, &b_mc_pt);
   fChain.SetBranchAddress("mc_m", &mc_m, &b_mc_m);
   fChain.SetBranchAddress("mc_eta", &mc_eta, &b_mc_eta);
   fChain.SetBranchAddress("mc_phi", &mc_phi, &b_mc_phi);
   fChain.SetBranchAddress("mc_status", &mc_status, &b_mc_status);
   fChain.SetBranchAddress("mc_barcode", &mc_barcode, &b_mc_barcode);
   fChain.SetBranchAddress("mc_pdgId", &mc_pdgId, &b_mc_pdgId);

  // create histograms
  TH1F *hHiggs_mass   = new TH1F("hHiggs_mass"  ,"Higgs boson mass",50,123., 127.);
  TH1F *hHiggs_pt     = new TH1F("hHiggs_pt"    ,"Higgs boson transverse momenta",50,0.0, 200.);
  TH1F *hHiggs_eta    = new TH1F("hHiggs_eta"   ,"Higgs boson pseudorapidity",50,-10.0, 10.0);

  // loop over all entries in the fChain
  Int_t nevent = fChain.GetEntries(); 
  for (Int_t i=0;i<nevent;i++) {
    fChain.GetEvent(i);  
    // loop over all particles in the event
    for(int ipart = 0; ipart < mc_n; ipart++){
      if( mc_pdgId->at(ipart) == 25 && mc_status->at(ipart) > 60 ){
        //fill histogram for Higgs boson mass and transverse momenta
	hHiggs_mass->Fill(mc_m->at(ipart));
	// here recalculate pt
	float pt= sqrt(mc_px->at(ipart) * mc_px->at(ipart) + mc_py->at(ipart) * mc_py->at(ipart));
	hHiggs_pt->Fill(pt);
	hHiggs_eta->Fill(mc_eta->at(ipart));
      }
    }
  }
  
  // draw your histogram
  //////////////////////////////////////////////////
  TCanvas *c = new TCanvas("c1","canvas",800,601);
  c->SetFillColor(0);
  gPad->SetFillColor(0);

  hHiggs_mass->Draw();
  c->Print("hHiggs_mass.eps"); 

  hHiggs_pt->Draw();
  c->Print("hHiggs_pt.eps"); 

  hHiggs_eta->Draw();
  c->Print("hHiggs_eta.eps"); 


  // open file to store your histogram
  TFile  outputFile("histos.root","RECREATE");
  outputFile.cd();
  hHiggs_mass->Write();
  hHiggs_pt->Write();
  hHiggs_eta->Write();
  outputFile.Close();

}