#define TMBTree_bu_cxx
#include "TMBTree_bu.h"
#include "tmb_tree/TMBBCJet.hpp"
#include "tmb_tree/TMBBCMuTagger.hpp"
#include "TLorentzVector.h"
#include "TH2.h"
#include "TCanvas.h"
#include "TRefArray.h"
#include "TString.h"
#include <iostream>

#ifdef __MAKECINT__
#pragma link C++ class TMBTree_bu;
#pragma link C++ class std::vector<TH1*>;
#endif

using namespace std;

void TMBTree_bu::Walk( TMBMCpart * mcpart, Bool_t debug )
{
   Int_t pdgid = mcpart->pdgid();
   if (debug) cout << " Next particle's PDG ID " << pdgid << endl; 
   // Track down the particle
   TMBMCvtx  * mcvtx = mcpart->getDMCvtx();
   if (!mcvtx) {
      if (debug) cout << "   particle is stable " << endl; 
      return;
   }
   Int_t ndaughters = mcvtx->ndaughters();
   if (debug) cout << "    particle has " << ndaughters << " daughters" << endl;
   if (ndaughters == 0) return;
   for (Int_t ii=0; ii<ndaughters; ii++) {
      TMBMCpart * mcpartd = mcvtx->getDaughter(ii);
      pdgid = mcpartd->pdgid();
      if (debug) 
	 cout << "      PDG ID of " << ii << "th daughter " << pdgid << endl;
      Walk (mcpartd, debug);
   }
}


void TMBTree_bu::Loop()
{
  //   In a ROOT session, you can do:
  //      Root > .L T.C+
  //      Root > T 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
  //       To read only selected branches, Insert statements like:
  // disable all branches:
  //    fChain->SetBranchStatus("*",0);  
  // activate branchname (DON'T FORGET THE "*" AT THE END!)
  //    fChain->SetBranchStatus("branchname*",1);  
  
  if (fChain == 0) return;
  
  
  // make structures for the information you 
  // want to put in your tuples
  // I dont know how to feed the tree Int_t's so
  // call everything a float (or double)

  // inclusive muon informatio
  typedef struct {
    Float_t pt, phi, eta, dca1, dca2, sctimeA, chi2;
    Float_t ncft, nsmt, nseg, nwhitA, nshitA, nwhitBC, nshitBC, centralr, charge;
  } my_muon;
  
  // dimuon system
  typedef struct {
    Float_t pt, phi, eta, y, mass, mkk;
    Float_t charge;
    Float_t pt_new, phi_new, eta_new, mass_new;
    Float_t lxy, svov, chisq, colin, znew, xzee;
    Float_t dtheta;
    Float_t ec1, ec2, ec4, ec6, ec15, dr, dz1, dz2;
  } my_dimuon;

  // generic event information
  typedef struct {
    Float_t evtno, runno, used_runno, nmuon, ncand, nDcand;
    Int_t  Fastz, LM_north, LM_south, Ahalo, Phalo;
  } my_event;
  
  // generic event information
  typedef struct {
    Float_t energy, ieta, iphi, ilyr;
  } my_emcell;

  // track information 
  typedef struct {
    Float_t pt, phi, eta, dca1, dca2;
    Float_t ncft, nsmt, charge;
  } my_track;

  // intermidiate resonance information 
  typedef struct {
    Float_t pt, phi, eta, mkpi, mpipi, mkk, mppi, dphi, deta, dr;
  } my_intermidiate;

  // b information 
  typedef struct {
    Float_t pt, phi, eta, mkpi, mpipi, mkk, mppi;
  } my_b;

  typedef struct {
    Float_t esum_n0, esum_n10, esum_n50, esum_n100;
    Float_t esum_s0, esum_s10, esum_s50, esum_s100;
    Float_t esum_0, esum_10, esum_50, esum_100;
  } my_esum;
  

  // declare structures you will use
  static my_muon All_muons;
  static my_muon muon2;
  static my_dimuon dimu;
  static my_event Event;
  static my_emcell EmCell;
  static my_esum  Esum;
  static my_intermidiate hh;
  static my_intermidiate bcand;


  // define your root tree
  TTree *t100 = new TTree("t100","muon data tree");
   
  // define a branch, feed it the struct, and give it variable names

  // fill once per muon
  t100->Branch("All_muons",&All_muons,"pt:phi:eta:dca1:dca2:sctimeA:chi2:ncft:nsmt:nseg:nwhitA:nshitA:nwhitBC:nshitBC:centralr:charge:sig1:sig2");
  
  
  // new tree
  // I define a new tree depending on when I want to fill the tree
  // fill once for every dimuon candidate
  TTree *t101 = new TTree("t101","dimuon data tree");

  t101->Branch("muon1",&All_muons,"pt1:phi1:eta1:dca11:dca21:sctimeA1:chi21:ncft1:nsmt1:nseg1:nwhitA1:nshitA1:nwhitBC1:nshitBC1:centralr1:charge1");
  t101->Branch("muon2",&muon2,"pt2:phi2:eta2:dca12:dca22:sctimeA2:chi22:ncft2:nsmt2:nseg2:nwhitA2:nshitA2:nwhitBC2:nshitBC2:centralr2:charge2");
  t101->Branch("dimu",&dimu,"pt:phi:eta:y:mass:mkk:charge:pt_new:phi_new:eta_new:mass_new:lxy:svov:chisq:colin:znew:xzee:dtheta:ec1:ec2:ec4:ec6:ec15:dr:dz1:dz2");
  t101->Branch("event",&Event,"evtno:runno:used_runno:nmuon:ncand:ndcand:Fastz/I:LM_north/I:LM_south/I:Ahalo/I:Phalo/I");
  t101->Branch("esum",&Esum,"esum_n0:esum_n10:esum_n50:esum_n100:esum_s0:esum_s10:esum_s50:esum_s100:esum_0:esum_10:esum_50:esum_100");


  //new tree
  // fill once for cell, for every dimuon candidate
  TTree *t103 = new TTree("t103","cell tree");

  t103->Branch("emcell",&EmCell,"energy:ieta:iphi:ilyr");
  t103->Branch("event",&Event,"evtno:runno:used_runno:nmuon:ncand:ndcand:Fastz/I:LM_north/I:LM_south/I:Ahalo/I:Phalo/I");
  

  
  //new tree
  //TTree *t102 = new TTree("t102","event tree");

  //t102->Branch("event",&Event,"evtno:runno:nmuon:ncand:ndcand:Fastz:LM_north:LM_south:Ahalo:Phalo");






  // histograms
  TH1F* h1 = new TH1F("h1","m(Kpi)1",200,0.,2.0);
  TH1F* h2 = new TH1F("h2","m(Kpi)2",200,0.,2.0);
  TH1F* h3 = new TH1F("h3","m(KK)",200,0.,2.0);
  TH1F* h4 = new TH1F("h4","m(pipi)",200,0.,2.0);
  TH1F* h5 = new TH1F("h5","m(p pi)1",200,0.,2.0);
  TH1F* h6 = new TH1F("h6","m(p pi)2",200,0.,2.0);

  TH1F* h7 = new TH1F("h7","track pt",1000,0.,20.0);
  TH1F* h8 = new TH1F("h8","n track",200,0.,200.0);

  
  /// Adding my own histograms ///
  
  //EM cells studies...
  TH1F* h80 = new TH1F("h80","n emcells",200,0.,300.0);
  TH1F* h81 = new TH1F("h81","energy",200,0.,20.0);
  TH2F* h82 = new TH2F("h82","etaXE",200,-40,40,200,-5,5);
  TH1F* h83 = new TH1F("h83","layer",200,0,20);
  
  TH2F* h85_lego = new TH2F("h85_lego","eta X sum_energy",100,0,20,100,-35,-28);
  
  TH1F* h85 = new TH1F("h85","sum of energy on north side",100,0,200);
  TH1F* h85_log = new TH1F("h85_log","sum of energy on north side (log10 scale)",100,-4,10);
  
  TH1F* h85Noff = new TH1F("h85Noff","sum of energy on north side with N off",100,-4,10);
  TH1F* h85Noff_lgy = new TH1F("h85Noff_lgy","sum of energy on north side with N off",100,0,200);
  TH1F* h85Noff_gapselection = new TH1F("h85Noff_gapselection","sum of energy (log10(sumn_E)<1.0) on north side with N off",100,-2,4);
  
  TH1F* h85Soff = new TH1F("h85Soff","energy on north side S off",100,-4,10);
  TH1F* h85Soff_lgy = new TH1F("h85Soff_lgy","energy on north side S off",100,0,200);
  TH1F* h85NSoff = new TH1F("h85NSoff","energy on north side both detectors off",100,-4,10);
  TH1F* h85NSoff_lgy = new TH1F("h85NSoff_lgy","energy on north side both detectors off",100,0,200);
  TH1F* h85fastZ = new TH1F("h85fastZ","energy on north side both detectors on",100,-4,10);
  TH1F* h85fastZ_lgy = new TH1F("h85fastZ_lgy","energy on north side both detectors on",100,0,200);
  
  TH2F* h87_lego = new TH2F("h87_lego","eta X sum_energy",100,0,20,100,28,35);
  
  TH1F* h87 = new TH1F("h87","sum of energiy on south side",100,0,200);
  TH1F* h87_log = new TH1F("h87_log","sum of energy on south side (log10 scale)",100,-4,10);
  
  TH1F* h87Soff = new TH1F("h87Soff","energy on south side S off",100,-4,10);
  TH1F* h87Soff_lgy = new TH1F("h87Soff_lgy","energy on south side S off",100,0,200);
  TH1F* h87Soff_gapselection = new TH1F("h87Soff_gapselection","sum of energy (log10(sums_E)<1.0) on south side with S off",100,-4,10);
  
  TH1F* h87Noff = new TH1F("h87Noff","energy on south side N off",100,-4,10);
  TH1F* h87Noff_lgy = new TH1F("h87Noff_lgy","energy on south side N off",100,0,200);
  TH1F* h87NSoff = new TH1F("h87NSoff","energy on south side both detectors off",100,-4,10);
  TH1F* h87NSoff_lgy = new TH1F("h87NSoff_lgy","energy on south side both detectors off",100,0,200);
  TH1F* h87fastZ = new TH1F("h87fastZ","energy on south side both detectors on",100,-4,10);
  TH1F* h87fastZ_lgy = new TH1F("h87fastZ_lgy","energy on south side both detectors on",100,0,200);
  
  TH1F* h88 = new TH1F("h88","total energy sum",100,0,200);  
  TH1F* h88_log = new TH1F("h88_log","sum of energy on north and south side (log10 scale)",100,-4,10);
  //end of EM cells histograms 
  
  TH2F* etabyphi = new TH2F("etabyphi","eta by phi",100,-6,6,100,0,7);
  TH1F* pt_distribution = new TH1F("pt_distribution","pt distribution for dimuons candidates",100,0,5);
  
  //gap information
  TH1F* mass_gapN = new TH1F("mass_gapN","m(dimu) with gapN",100,1,5);
  TH1F* mass_NoffSon = new TH1F("mass_NoffSon","m(dimu) with LM_north off and LM_south on",100,1,5);
  TH1F* eta_gapN = new TH1F("eta_gapN","eta(dimu) with gapN",100,-6,6);
  TH1F* eta_NoffSon = new TH1F("eta_NoffSon","eta(dimu) with LM_north off and LM_south on",100,-6,6);
  
  TH1F* mass_gapS = new TH1F("mass_gapS","m(dimu) with gapS",100,1,5);
  TH1F* mass_NonSoff = new TH1F("mass_NonSoff","m(dimu) with LM_north on and LM_south off",100,1,5);
  TH1F* eta_gapS = new TH1F("eta_gapS","eta(dimu) with gapS",100,-6,6); 
  TH1F* eta_NonSoff = new TH1F("eta_NonSoff","eta(dimu) with LM_north on and LM_south off",100,-6,6);
  
  TH1F* mass_gapNS = new TH1F("mass_gapNS","m(dimu) with gapNS",100,1,5);
  TH1F* eta_gapNS = new TH1F("eta_gapNS","m(dimu) with gapNS",100,-6,6);
  
  TH1F* mass_fastZ = new TH1F("mass_fastZ","m(dimu) in FastZ coincidence",100,1,5);
  TH1F* mass_Phalo = new TH1F("mass_Phalo","m(dimu) in Phalo coincidence",100,1,5);
  TH1F* mass_Ahalo = new TH1F("mass_Ahalo","m(dimu) in Ahalo coincidence",100,1,5);
  
  TH1F* eta_forneg_mass_NoffSon = new TH1F("eta_forneg_mass_NoffSon","mass in -2<eta<-1 (Noff Son)",100,1,5);
  TH1F* eta_centneg_mass_NoffSon = new TH1F("eta_centneg_mass_NoffSon","mass in -1<eta<0 (Noff Son)",100,1,5);
  TH1F* eta_centpos_mass_NoffSon = new TH1F("eta_centpos_mass_NoffSon","mass in 0<eta<1 (Noff Son)",100,1,5);
  TH1F* eta_forpos_mass_NoffSon = new TH1F("eta_forpos_mass_NoffSon","mass in 1<eta<2 (Noff Son)",100,1,5);
  
  TH1F* eta_forneg_mass_NonSoff = new TH1F("eta_forneg_mass_NonSoff","mass in -2<eta<-1 (Non Soff)",100,1,5);
  TH1F* eta_centneg_mass_NonSoff = new TH1F("eta_centneg_mass_NonSoff","mass in -1<eta<0 (Non Soff)",100,1,5);
  TH1F* eta_centpos_mass_NonSoff = new TH1F("eta_centpos_mass_NonSoff","mass in 0<eta<1 (Non Soff)",100,1,5);
  TH1F* eta_forpos_mass_NonSoff = new TH1F("eta_forpos_mass_NonSoff","mass in 1<eta<2 (Non Soff)",100,1,5);
  
  TH1F* eta_forneg_mass_gapNS = new TH1F("eta_forneg_mass_gapNS","mass in -2<eta<-1 (double gap)",100,1,5);
  TH1F* eta_centneg_mass_gapNS = new TH1F("eta_centneg_mass_gapNS","mass in -1<eta<0 (double gap)",100,1,5);
  TH1F* eta_centpos_mass_gapNS = new TH1F("eta_centpos_mass_gapNS","mass in 0<eta<1 (double gap)",100,1,5);
  TH1F* eta_forpos_mass_gapNS = new TH1F("eta_forpos_mass_gapNS","mass in 1<eta<2 (double gap)",100,1,5);
  
  TH1F* eta_forneg_mass_fastZ = new TH1F("eta_forneg_mass_fastZ","mass in -2<eta<-1 (fastZ)",100,1,5);
  TH1F* eta_centneg_mass_fastZ = new TH1F("eta_centneg_mass_fastZ","mass in -1<eta<0 (fastZ)",100,1,5);
  TH1F* eta_centpos_mass_fastZ = new TH1F("eta_centpos_mass_fastZ","mass in 0<eta<1 (fastZ)",100,1,5);
  TH1F* eta_forpos_mass_fastZ = new TH1F("eta_forpos_mass_fastZ","mass in 1<eta<2 (fastZ)",100,1,5);
  
  TH1F* y_forneg_mass_NoffSon = new TH1F("y_forneg_mass_NoffSon","mass in -2<y<-1 (Noff Son)",100,1,5);
  TH1F* y_centneg_mass_NoffSon = new TH1F("y_centneg_mass_NoffSon","mass in -1<y<0 (Noff Son)",100,1,5);
  TH1F* y_centpos_mass_NoffSon = new TH1F("y_centpos_mass_NoffSon","mass in 0<y<1 (Noff Son)",100,1,5);
  TH1F* y_forpos_mass_NoffSon = new TH1F("y_forpos_mass_NoffSon","mass in 1<y<2 (Noff Son)",100,1,5);
  
  TH1F* y_forneg_mass_NonSoff = new TH1F("y_forneg_mass_NonSoff","mass in -2<y<-1 (Non Soff)",100,1,5);
  TH1F* y_centneg_mass_NonSoff = new TH1F("y_centneg_mass_NonSoff","mass in -1<y<0 (Non Soff)",100,1,5);
  TH1F* y_centpos_mass_NonSoff = new TH1F("y_centpos_mass_NonSoff","mass in 0<y<1 (Non Soff)",100,1,5);
  TH1F* y_forpos_mass_NonSoff = new TH1F("y_forpos_mass_NonSoff","mass in 1<y<2 (Non Soff)",100,1,5);
  
  TH1F* y_forneg_mass_gapNS = new TH1F("y_forneg_mass_gapNS","mass in -2<y<-1 (double gap)",100,1,5);
  TH1F* y_centneg_mass_gapNS = new TH1F("y_centneg_mass_gapNS","mass in -1<y<0 (double gap)",100,1,5);
  TH1F* y_centpos_mass_gapNS = new TH1F("y_centpos_mass_gapNS","mass in 0<y<1 (double gap)",100,1,5);
  TH1F* y_forpos_mass_gapNS = new TH1F("y_forpos_mass_gapNS","mass in 1<ey<2 (double gap)",100,1,5);
  
  TH1F* y_forneg_mass_fastZ = new TH1F("y_forneg_mass_fastZ","mass in -2<y<-1 (fastZ)",100,1,5);
  TH1F* y_centneg_mass_fastZ = new TH1F("y_centneg_mass_fastZ","mass in -1<y<0 (fastZ)",100,1,5);
  TH1F* y_centpos_mass_fastZ = new TH1F("y_centpos_mass_fastZ","mass in 0<y<1 (fastZ)",100,1,5);
  TH1F* y_forpos_mass_fastZ = new TH1F("y_forpos_mass_fastZ","mass in 1<y<2 (fastZ)",100,1,5);
  
  
  TH1F* y_firstneg_mass_NoffSon = new TH1F("y_firstneg_mass_NoffSon","mass in -0.6<y<0 (NoffSon)",100,1,5);
  TH1F* y_secondneg_mass_NoffSon = new TH1F("y_secondneg_mass_NoffSon","mass in -1.2<y<-0.6 (NoffSon)",100,1,5);
  TH1F* y_thirdneg_mass_NoffSon = new TH1F("y_thirdneg_mass_NoffSon","mass in -1.8<y<-1.2 (NoffSon)",100,1,5);
  
  TH1F* y_firstpos_mass_NoffSon = new TH1F("y_firstpos_mass_NoffSon","mass in 0<y<0.6 (NoffSon)",100,1,5);
  TH1F* y_secondpos_mass_NoffSon = new TH1F("y_secondpos_mass_NoffSon","mass in 0.6<y<1.2 (NoffSon)",100,1,5);
  TH1F* y_thirdpos_mass_NoffSon = new TH1F("y_thirdpos_mass_NoffSon","mass in 1.2<y<1.8 (NoffSon)",100,1,5);
  
  TH1F* y_firstneg_mass_NonSoff = new TH1F("y_firstneg_mass_NonSoff","mass in -0.6<y<0 (NonSoff)",100,1,5);
  TH1F* y_secondneg_mass_NonSoff = new TH1F("y_secondneg_mass_NonSoff","mass in -1.2<y<-0.6 (NonSoff)",100,1,5);
  TH1F* y_thirdneg_mass_NonSoff = new TH1F("y_thirdneg_mass_NonSoff","mass in -1.8<y<-1.2 (NonSoff)",100,1,5);
  
  TH1F* y_firstpos_mass_NonSoff = new TH1F("y_firstpos_mass_NonSoff","mass in 0<y<0.6 (NonSoff)",100,1,5);
  TH1F* y_secondpos_mass_NonSoff = new TH1F("y_secondpos_mass_NonSoff","mass in 0.6<y<1.2 (NonSoff)",100,1,5);
  TH1F* y_thirdpos_mass_NonSoff = new TH1F("y_thirdpos_mass_NonSoff","mass in 1.2<y<1.8 (NonSoff)",100,1,5);
  
  TH1F* y_firstneg_mass_gapNS = new TH1F("y_firstneg_mass_gapNS","mass in -0.6<y<0 (gapNS)",100,1,5);
  TH1F* y_secondneg_mass_gapNS = new TH1F("y_secondneg_mass_gapNS","mass in -1.2<y<-0.6 (gapNS)",100,1,5);
  TH1F* y_thirdneg_mass_gapNS = new TH1F("y_thirdneg_mass_gapNS","mass in -1.8<y<-1.2 (gapNS)",100,1,5);
  
  TH1F* y_firstpos_mass_gapNS = new TH1F("y_firstpos_mass_gapNS","mass in 0<y<0.6 (gapNS)",100,1,5);
  TH1F* y_secondpos_mass_gapNS = new TH1F("y_secondpos_mass_gapNS","mass in 0.6<y<1.2 (gapNS)",100,1,5);
  TH1F* y_thirdpos_mass_gapNS = new TH1F("y_thirdpos_mass_gapNS","mass in 1.2<y<1.8 (gapNS)",100,1,5);
  
  TH1F* y_firstneg_mass_fastZ = new TH1F("y_firstneg_mass_fastZ","mass in -0.6<y<0 (fastZ)",100,1,5);
  TH1F* y_secondneg_mass_fastZ = new TH1F("y_secondneg_mass_fastZ","mass in -1.2<y<-0.6 (fastZ)",100,1,5);
  TH1F* y_thirdneg_mass_fastZ = new TH1F("y_thirdneg_mass_fastZ","mass in -1.8<y<-1.2 (fastZ)",100,1,5);
  
  TH1F* y_firstpos_mass_fastZ = new TH1F("y_firstpos_mass_fastZ","mass in 0<y<0.6 (fastZ)",100,1,5);
  TH1F* y_secondpos_mass_fastZ = new TH1F("y_secondpos_mass_fastZ","mass in 0.6<y<1.2 (fastZ)",100,1,5);
  TH1F* y_thirdpos_mass_fastZ = new TH1F("y_thirdpos_mass_fastZ","mass in 1.2<y<1.8 (fastZ)",100,1,5);
  
  TH1F* y_first_mass_NoffSon = new TH1F("y_first_mass_NoffSon","mass in -2<y<-1 (NoffSon)",100,1,5);
  TH1F* y_second_mass_NoffSon = new TH1F("y_second_mass_NoffSon","mass in -2<y<-1 (NoffSon)",100,1,5);
  TH1F* y_third_mass_NoffSon = new TH1F("y_third_mass_NoffSon","mass in -2<y<-1 (NoffSon)",100,1,5);
  
  TH1F* y_first_mass_NonSoff = new TH1F("y_first_mass_NonSoff","mass in 0<|y|<0.6 (NonSoff)",100,1,5);
  TH1F* y_second_mass_NonSoff = new TH1F("y_second_mass_NonSoff","mass in 0.6<|y|<1.2 (NonSoff)",100,1,5);
  TH1F* y_third_mass_NonSoff = new TH1F("y_third_mass_NonSoff","mass in 1.2<|y|<1.8 (NonSoff)",100,1,5);
  
  TH1F* y_first_mass_gapNS = new TH1F("y_first_mass_gapNS","mass in 0<|y|<0.6 (gapNS)",100,1,5);
  TH1F* y_second_mass_gapNS = new TH1F("y_second_mass_gapNS","mass in 0.6<|y|<1.2 (gapNS)",100,1,5);
  TH1F* y_third_mass_gapNS = new TH1F("y_third_mass_gapNS","mass in 1.2<|y|<1.8 (gapNS)",100,1,5);
  
  TH1F* y_first_mass_fastZ = new TH1F("y_first_mass_fastZ","mass in 0<|y|<0.6 (fastZ)",100,1,5);
  TH1F* y_second_mass_fastZ = new TH1F("y_second_mass_fastZ","mass in 0.6<|y|<1.2 (fastZ)",100,1,5);
  TH1F* y_third_mass_fastZ = new TH1F("y_third_mass_fastZ","mass in 1.2<|y|<1.8 (fastZ)",100,1,5);
  //end of gap information
  
  
  
  // fChain is the chain of files you got (I think)
  Int_t nentries = Int_t(fChain->GetEntries());
  
  Int_t nbytes = 0, nb = 0;
  

  //Defining the AndOrTerms
  Int_t aot_30[256];
  for(Int_t jjj = 0; jjj < 256 ; jjj++){aot_30[jjj] = 0;}
  
  //Looking for PV
  TString used_evt;
  
  
  // loop over all events
  for (Int_t jentry=0; jentry<nentries; jentry++) {
    
    
    //in case of a TChain, ientry is the entry number in the current file
    //    Int_t ientry = LoadTree(jentry); 
    // ???
    nb = fChain->GetEntry(jentry);   
    nbytes += nb;
    
    
    if(jentry%500 ==0) cout<<"events analyzed "<<jentry<<endl;


    // Vertex information
    Int_t nfVrts = fVrts->GetLast();
    Int_t vert_type = 0;
    Double_t ntr[10000],vertZ[10000], Dz = 0;
    
    for(Int_t i=0; i<10000; i++){
      ntr[i] = 0.0;
      vertZ[i] = 0.0;
    } 
       

    
      
    cout << jentry << " with " << nfVrts << " vertices " << endl;
	
    used_evt = "not_useful";
    for(Int_t tq=0; tq<nfVrts; ++tq){
	
	    TMBVrts* Vert = ( TMBVrts* ) fVrts->At(tq);
	  
	    ntr[tq] = Vert->Ntr();
	    vertZ[tq] = Vert->vz();
	    
	    if(nfVrts==1 && ntr[0]>2){
	    
	      used_evt = "1PV";
	      //use_evt[0] = "1PV"
	      //cout << nfVrts << " 1PV " << endl;
	      
	    }
	    if(nfVrts>1 && tq>0){
	    
	       if(ntr[tq]<=2){
	       
	       used_evt = "1PV";
	       //use_evt[tq] = "1PV";
	       //cout << nfVrts << "  " << tq << " trk<=2 => 1PV " << endl;
	       
	       }
	       else{
	       
	         Dz=abs(vertZ[tq] - vertZ[0]);
		 
		 if(Dz<0.375){
		 
		 used_evt = "1PV";
		 //use_evt[tq] = "1PV";
		 //cout << nfVrts <<  "  " << tq << " Dz " << Ddz << " 0.375 => 1PV " << endl;
		 
		 }
		 else{
		 
		   used_evt = "not_useful";
		   //cout << nfVrts << " just leave it " << endl;
		   tq=nfVrts;
		 
		 }
	       
	       }
	    
	    }
	    //if(tq==nfVrts-1){cout << used_evt << endl;}    
    }
    cout << used_evt << endl;
    if(used_evt=="1PV"){
  


    
            //Looking at the Glob Class
            Int_t run_number = fGlob->runno();
            Int_t event_number = fGlob->evtno();
    
            Event.runno = run_number;
            Event.evtno = event_number;
	    
	    //cout << event_number << endl;
	    

            // fMuons and fTracks are copies of the
            // root tree info for muons and tracks
            // see TMBTree_bu.h
            Int_t nfMuons = fMuons->GetLast();
            Int_t nfTrk = fTrks->GetLast();
            Int_t nfJets = fJets->GetLast();
    
    
            // Trigger Information
            Int_t nfTrig = fTrig->GetLast();
            TString trig_name;
    

            Int_t nmuons = 0;
    
            // This is how I check against double counting muons
            // also, as far as I can tell, TArrayD is the equivalent of a 
            // stl::vector in root
            TArrayD used_muons(nfMuons);

            Int_t n_used_muons = 0;
        
            Int_t ncand = 0;
            Int_t nDcand = 0;
    
    
            //Looking at the AndOrTerms for each event 
            for (Int_t ttt = 0; ttt<256; ttt++){

              aot_30[ttt] = fL1Muon->L1MuGetAndOrTerms(ttt);

              if(ttt== 217){Event.Fastz = aot_30[ttt];}
              if(ttt== 220){Event.LM_south = aot_30[ttt];}
              if(ttt== 221){Event.LM_north = aot_30[ttt];}
              if(ttt== 222){Event.Ahalo = aot_30[ttt];}
              if(ttt== 223){Event.Phalo = aot_30[ttt];}

            }	
    

            // copying also fEmCells from the root tree
            Int_t nfEmCell = fEmCells->GetLast();
            Double_t energia;
            Double_t i_eta;
            Double_t i_phi;
            Int_t i_layer;


            Double_t sumn_E = 0.0;
            Double_t sums_E = 0.0;
            Double_t sumns_E = 0.0;
    

            Double_t sumn_E0 = 0.0;
            Double_t sumn_E10 = 0.0;
            Double_t sumn_E50 = 0.0;
            Double_t sumn_E100 = 0.0;
	    
	    Double_t sums_E0 = 0.0;
            Double_t sums_E10 = 0.0;
            Double_t sums_E50 = 0.0;
            Double_t sums_E100 = 0.0;
	    
	    Double_t sum_E0 = 0.0;
            Double_t sum_E10 = 0.0;
            Double_t sum_E50 = 0.0;
            Double_t sum_E100 = 0.0;



	    //loop over triggers
            for(Int_t b=0; b<nfTrig; ++b){
      
	      TMBTrig* Trigger = ( TMBTrig* ) fTrig->At(b);
	
	      trig_name = Trigger->getTrgName();

	      //if(trig_name!="2MT1_2TRK_GAPN"&&trig_name!="2MT1_2TRK_GAPS"&&trig_name!="2MT1_2TRK_GAPSN"){
	      //if(trig_name=="2MT1_2TRK_GAPN" || trig_name=="2MT1_2TRK_GAPS" || trig_name=="2MT1_2TRK_GAPSN"){
	      //if(trig_name=="2MU_A_L2M0" || trig_name=="2MU_A_L2ETAPHI" || trig_name=="2MT1_2TRK_GAPN" || trig_name=="2MT1_2TRK_GAPS" ||trig_name=="2MT1_2TRK_GAPSN"){
	      
		
	      if(trig_name=="2MU_A_L2M0" || trig_name=="2MU_A_L2ETAPHI"){  
		
    
		//}
	    
	    
		Event.used_runno = run_number;
    
    
		//I only want to look for runs after  
		// loop over all muons
		for(Int_t ii = 0; ii < nfMuons; ++ii){
	      


                TMBMuon* Muon = ( TMBMuon* ) fMuons->At(ii);

                // Muon is now an object of type TMBMuon 
                // it contains all the information about the muon stored in the 
                // tree at location ii
                // For a list of variables, see
      
                // http://www-d0.fnal.gov/d0dist/dist/packages/tmb_tree/v00-07-10/tmb_tree/TMBMuon.hpp
      
                if(Muon){  // you got the pointer but you have to check that its
	          // pointing to something or you get a segmentation violation
	
	          TMBTrks* MatchedTrack = ( TMBTrks* ) Muon->GetChargedTrack();
	          // MatchedTrack is now an object of type TMBTrks
	          // this is the track object that is best associated 
	          // with the muon candidate
	          // see the above web site with TMBMuon.hpp changed to TMBTrks.hpp
	
	          if(MatchedTrack){

	            if(Muon->centralrank()==1){  // muon variable soon to be redundant

	              // here I'm checking that I havent used this track already
	              // i.e. two muons were associated with this track
	              Int_t OK_to_use = 1;
	    
	              for( Int_t iii = 0; iii < n_used_muons; iii++){

	                // I use px instead of the index just incase the indecies 
	                // have been screwed up
	                if(MatchedTrack->px() == used_muons.At(iii)) OK_to_use = 0;

	              }
	 
	              if(OK_to_use == 1){
	      
	                used_muons.AddAt(MatchedTrack->px(),n_used_muons);

	                n_used_muons++;

	                // this is for checking later
	                Double_t px1 = MatchedTrack->px();

	                // fill the struct with the variables
	                All_muons.pt = MatchedTrack->pT();
	                All_muons.eta = MatchedTrack->eta();	      
	                All_muons.phi = MatchedTrack->phi();
	                All_muons.dca1 = MatchedTrack->dca1();
	                All_muons.dca2 = MatchedTrack->dca2();
	                All_muons.nsmt = MatchedTrack->nsmt();
	                All_muons.ncft = MatchedTrack->ncft();
	                All_muons.charge = MatchedTrack->charge();
	                All_muons.nseg = Muon->nseg();
	                All_muons.sctimeA = Muon->sctimeA();
	                All_muons.chi2 = Muon->chisqloc();
	                All_muons.centralr = Muon->centralrank();
	      
	                // unpack the muon nhit variable into the scintillator
	                // and wire hit information for each layer
	                Int_t idnhit = Muon->nhit();
	       
	                Int_t nshitBC = idnhit/10000;
	                Int_t nshitA = (idnhit - 10000*nshitBC)/1000;
	                Int_t nwhitBC = (idnhit - 10000*nshitBC - 1000*nshitA)/10;
	                Int_t nwhitA = (idnhit - 10000*nshitBC - 1000*nshitA-10*nwhitBC);
	      
	                All_muons.nshitA = nshitA;
	                All_muons.nwhitA = nwhitA;
	                All_muons.nshitBC = nshitBC;
	                All_muons.nwhitBC = nwhitBC;
	      
	                nmuons++;
	      
	                // make a 4 vector of the muon
	                TLorentzVector p4mu1(0.0, 0.0, 0.0, 0.0);     
	      
	                p4mu1.SetPx(MatchedTrack->px());
	                p4mu1.SetPy(MatchedTrack->py());
	                p4mu1.SetPz(MatchedTrack->pz());
			      
	                // assign muon mass hypothesis
	                Float_t E_sig_mu = sqrt((MatchedTrack->pT()*MatchedTrack->pT()) +
			       	               (MatchedTrack->pz()*MatchedTrack->pz()) +
				                0.0112);

	                p4mu1.SetE(E_sig_mu);
	      
	                // assign kaon mass hypothesis
	                Float_t E_sig_k = sqrt((MatchedTrack->pT()*MatchedTrack->pT()) +
				              (MatchedTrack->pz()*MatchedTrack->pz()) +
				               0.244);
			
	                TLorentzVector p4k1(p4mu1.X(),p4mu1.Y(),p4mu1.Z(),E_sig_k);          

	      
	                // this is filled once for every muon
	                // to make it smaller, I only sample some events
	      
	                // if(jentry%10000 == 0) t100->Fill();
	      
	      
	                // make some basic cuts to get rid of junk
	                //	      if(All_muons.pt > 1.0 &&
	                //	 All_muons.nseg > 0 &&
	                //	 All_muons.ncft > 0 &&
	                //	 All_muons.nsmt > 2 &&
	                //	 TMath::Abs(All_muons.sctimeA) < 20){
	      
	                //double counting checking for the second muon
	                TArrayD used_muons2(nfMuons);

	                Int_t n_used_muons2 = 0;


	                // second muon loop, second verse, same as the first
	                for(Int_t iii = ii+1; iii < nfMuons; ++iii){
			  
		        TMBMuon* Muon2 = ( TMBMuon* ) fMuons->At(iii);
		
		        if(Muon2){
		  
		          TMBTrks* MatchedTrack2 = ( TMBTrks* ) Muon2->GetChargedTrack();
		  
		          if(MatchedTrack2){

		            if(Muon2->centralrank()==1){

		              Int_t OK_to_use2 = 1;
		      
		              for( Int_t iiii = 0; iiii < n_used_muons; iiii++){

			        if(MatchedTrack2->px() == used_muons.At(iiii)) OK_to_use2 = 0;
			
		              }
	 
		              if(OK_to_use2 == 1){

			        for( Int_t iiii = 0; iiii < n_used_muons2; iiii++){

			          if(MatchedTrack2->px() == used_muons2.At(iiii)) OK_to_use2 = 0;

			        }

			        if(OK_to_use2 == 1){

			          used_muons2.AddAt(MatchedTrack2->px(),n_used_muons2);
			    
			          n_used_muons2++;
			    
			  
			          Double_t deta = MatchedTrack->eta() - MatchedTrack2->eta(); 
			  
			          Double_t dphi = 
			            TMath::Abs((MatchedTrack->phi() - MatchedTrack2->phi())); 
			  
			          if(dphi>3.14159) dphi = 2*(3.14159) - dphi;
			    
			          Double_t mydr = sqrt(dphi*dphi + deta*deta);
			    
			          //if(mydr<2){

			          Double_t px2 = MatchedTrack2->px();
			      
			          muon2.pt = MatchedTrack2->pT();
			          muon2.eta = MatchedTrack2->eta();	      
			          muon2.phi = MatchedTrack2->phi();
			          muon2.dca1 = MatchedTrack2->dca1();
			          muon2.dca2 = MatchedTrack2->dca2();
			          muon2.nsmt = MatchedTrack2->nsmt();
			          muon2.ncft = MatchedTrack2->ncft();
			          muon2.charge = MatchedTrack2->charge();
			          muon2.nseg = Muon2->nseg();
			          muon2.sctimeA = Muon2->sctimeA();
			          muon2.chi2 = Muon2->chisqloc();
			          muon2.centralr = Muon2->centralrank();
			  
			          Int_t idnhit2 = Muon2->nhit();
			  
			          Int_t nshitBC2 = idnhit2/10000;
			          Int_t nshitA2 = (idnhit2 - 10000*nshitBC2)/1000;
			          Int_t nwhitBC2 = (idnhit2 - 10000*nshitBC2 - 1000*nshitA2)/10;
			          Int_t nwhitA2 = (idnhit2 - 10000*nshitBC2 - 1000*nshitA2-10*nwhitBC2);
			  
			          muon2.nshitA = nshitA2;
			          muon2.nwhitA = nwhitA2;
			          muon2.nshitBC = nshitBC2;
			          muon2.nwhitBC = nwhitBC2;
			  
			          //			  if(//muon2.pt > 1.0 &&
			          //muon2.nseg > 0 &&
			          // muon2.ncft > 0 &&
			          //   muon2.nsmt > 2 &&
			          //   TMath::Abs(muon2.sctimeA) < 20){
			  
			          TLorentzVector p4mu2(0.0, 0.0, 0.0, 0.0);     
			  
			          p4mu2.SetPx(MatchedTrack2->px());
			          p4mu2.SetPy(MatchedTrack2->py());
			          p4mu2.SetPz(MatchedTrack2->pz());
			  
			          Float_t E_sig_mu2 = 
			            sqrt((MatchedTrack2->pT()*MatchedTrack2->pT()) +
				         (MatchedTrack2->pz()*MatchedTrack2->pz()) +
				         0.0112);
			  
			          p4mu2.SetE(E_sig_mu2);

			          Double_t dimu_angle = p4mu2.Angle(p4mu1.Vect());

			          Double_t dimu_dr = p4mu2.DeltaR(p4mu1);
			      			      
			          Float_t E_sig_k2 
			            = sqrt((MatchedTrack2->pT()*MatchedTrack2->pT()) +
				           (MatchedTrack2->pz()*MatchedTrack2->pz()) +
				           0.244);
			  
			          TLorentzVector p4k2(p4mu2.X(),p4mu2.Y(),p4mu2.Z(),E_sig_k2); 
			  
				        // add the two 4 vectors
			          TLorentzVector P4mumu(p4mu1.X() + p4mu2.X(),
						        p4mu1.Y() + p4mu2.Y(),
						        p4mu1.Z() + p4mu2.Z(),
						        p4mu1.T() + p4mu2.T());
			  
			          TLorentzVector P4kk(p4k1.X() + p4k2.X(),
					              p4k1.Y() + p4k2.Y(),
					              p4k1.Z() + p4k2.Z(),
					              p4k1.T() + p4k2.T());
			  
				        // fill the dimuon struct
			          dimu.mass = P4mumu.M();
			          dimu.mkk = P4kk.M();
			          dimu.eta = P4mumu.PseudoRapidity();
				  dimu.y = P4mumu.Rapidity();
			          dimu.phi = P4mumu.Phi();
			          dimu.pt = P4mumu.Pt();
			          dimu.charge = All_muons.charge + muon2.charge;
			          dimu.dtheta = dimu_angle;
			          dimu.dr = dimu_dr;

			  
			          if(MatchedTrack->pT()>MatchedTrack2->pT()){
			    
			            dimu.ec1 = Muon->EInCone1();
			            dimu.ec2 = Muon->EInCone2();
			            dimu.ec4 = Muon->EInCone4();
			            dimu.ec6 = Muon->EInCone6();
			            dimu.ec15 = Muon->EInCone15();
			          }
			          else{
			    
			            dimu.ec1 = Muon2->EInCone1();
			            dimu.ec2 = Muon2->EInCone2();
			            dimu.ec4 = Muon2->EInCone4();
			            dimu.ec6 = Muon2->EInCone6();
			            dimu.ec15 = Muon2->EInCone15();
			    
			          }


			          if(dimu.mass > 1.9 && dimu.mass < 4.1 && dimu.charge ==0){
			    
			            ncand++;
				    
				    
				    //loop over the calorimeter cells
				    for(Int_t a=0; a<nfEmCell; ++a){
      
				      TMBEmCells* EmCells = ( TMBEmCells* ) fEmCells->At(a);
	
				      energia = EmCells->E();
				      i_eta = EmCells->ieta();
				      i_phi = EmCells->iphi();
				      i_layer = EmCells->ilyr();				       
				      
				      Int_t fill_me = 0;
				      
				      if(i_layer>=1 && i_layer<=7){
					
					//north
					if(i_eta>=-35 && i_eta<=-28){
					  
					  sumn_E = sumn_E + energia;

					  h85_lego->Fill(energia,i_eta);
					  h85_lego->SetOption("lego");
					  
					  EmCell.energy = energia;
					  EmCell.ieta = i_eta;
					  EmCell.iphi = i_phi;
					  EmCell.ilyr = i_layer;
					  
					  fill_me = 1;

					}
					
					//south
					if(i_eta>=28 && i_eta<=35){
					  
					  sums_E = sums_E + energia;

					  h87_lego->Fill(energia,i_eta);
					  h87_lego->SetOption("lego");
					  
					  EmCell.energy = energia;
					  EmCell.ieta = i_eta;
					  EmCell.iphi = i_phi;
					  EmCell.ilyr = i_layer;
					  
					  fill_me = 1;
					  
					}
					
					if(fill_me == 1) t103->Fill(); 

						 
					//No threshold
					//north
					if(i_eta>=-35 && i_eta<=-28){sumn_E0 = sumn_E0 + energia;}
					  
					//south
					if(i_eta>=28 && i_eta<=35){sums_E0 = sums_E0 + energia;}

					
					//energy above 10MeV
					if(energia>0.01){
					  
					  //north
					  if(i_eta>=-35 && i_eta<=-28){sumn_E10 = sumn_E10 + energia;}
					  
					  //south
					  if(i_eta>=28 && i_eta<=35){sums_E10 = sums_E10 + energia;}
					  
					}
					
					//energy above 50MeV
					if(energia>0.05){
					  
					  //north
					  if(i_eta>=-35 && i_eta<=-28){sumn_E50 = sumn_E50 + energia;}
					  
					  //south
					  if(i_eta>=28 && i_eta<=35){sums_E50 = sums_E50 + energia;}
					  
					}
					
					//energy above 100MeV
					if(energia>0.1){
					  
					  //north
					  if(i_eta>=-35 && i_eta<=-28){sumn_E100 = sumn_E100 + energia;}
					  
					  //south
					  if(i_eta>=28 && i_eta<=35){sums_E100 = sums_E100 + energia;}
					  
					}
					
				      }
				    }
				    
				    sum_E0 = sumn_E0 + sums_E0;
				    sum_E10 = sumn_E10 + sums_E10;
				    sum_E50 = sumn_E50 + sums_E50;
				    sum_E100 = sumn_E100 + sums_E100;
				    
				    Esum.esum_n0 = sumn_E0;
				    Esum.esum_s0 = sums_E0;
				    Esum.esum_0 = sum_E0;
				    
				    Esum.esum_n10 = sumn_E10;
				    Esum.esum_s10 = sums_E10;
				    Esum.esum_10 = sum_E10;
				    
				    Esum.esum_n50 = sumn_E50;
				    Esum.esum_s50 = sums_E50;
				    Esum.esum_50 = sum_E50;
				    
				    Esum.esum_n100 = sumn_E100;
				    Esum.esum_s100 = sums_E100;
				    Esum.esum_100 = sum_E100;
				    
				    
				    h85->Fill(sumn_E);
				    h85_log->Fill(log10(sumn_E));
				   	    
	                            h87->Fill(sums_E);
				    h87_log->Fill(log10(sums_E));
	    
	                            sumns_E = sumn_E + sums_E;
                                    h88->Fill(sumns_E);
				    h88_log->Fill(log10(sumns_E));
	    
	                            //North Information
                                    if(Event.LM_north==0 && Event.LM_south==1 && Event.Fastz==0 && Event.Ahalo==0 && Event.Phalo==0){
                                     
				      h85Noff->Fill(log10(sumn_E));
				      h85Noff_lgy->Fill(sumn_E);
				     
				      if(log10(sumn_E)<1.0){h85Noff_gapselection->Fill(log10(sumn_E));}
				   
                                    }
                                    if(Event.LM_north==1 && Event.LM_south==0 && Event.Fastz==0 && Event.Ahalo==0 && Event.Phalo==0){
                                      
				      h85Soff->Fill(log10(sumn_E));
				      h85Soff_lgy->Fill(sumn_E);
				      
                                    }
				    
				    if(Event.LM_north==0 && Event.LM_south==0 && Event.Fastz==0 && Event.Ahalo==0 && Event.Phalo==0){
                                      
				      h85NSoff->Fill(log10(sumn_E));
				      h85NSoff_lgy->Fill(sumn_E);
				      
                                    }
				    
                                    if(Event.LM_north==1 && Event.LM_south==1 && Event.Fastz==1){
                                      
				      h85fastZ->Fill(log10(sumn_E));
				      h85fastZ_lgy->Fill(sumn_E);
				      
                                    }      
	    
	                            //South Information
                                    if(Event.LM_north==0 && Event.LM_south==1 && Event.Fastz==0 && Event.Ahalo==0 && Event.Phalo==0){
                                      
				      h87Noff->Fill(log10(sums_E));
				      h87Noff_lgy->Fill(sums_E);
                                    
				    }
                                    
				    if(Event.LM_north==1 && Event.LM_south==0 && Event.Fastz==0 && Event.Ahalo==0 && Event.Phalo==0){
                                     
				      h87Soff->Fill(log10(sums_E));
				      h87Soff_lgy->Fill(sums_E);
				     
				      if(log10(sums_E)<1.0){h87Soff_gapselection->Fill(log10(sums_E));}
				   
                                    }
				    
				    if(Event.LM_north==0 && Event.LM_south==0 && Event.Fastz==0 && Event.Ahalo==0 && Event.Phalo==0){
                                      
				      h87NSoff->Fill(log10(sums_E));
				      h87NSoff_lgy->Fill(sums_E);
                                    
				    }
                                    
				    if(Event.LM_north==1 && Event.LM_south==1 && Event.Fastz==1){
                                      
				      h87fastZ->Fill(log10(sums_E));
				      h87fastZ_lgy->Fill(sums_E);
                                    
				    }



			            //
				        // vertex stuff for the two muons
				        //
			    
			    
				        // fill the dimuon struct
			            dimu.mass_new = 0.0;
			            dimu.eta_new = 0.0;
			            dimu.phi_new = 0.0;
			            dimu.pt_new = 0.0;
			            dimu.lxy = 0.0;
			            dimu.svov = 0.0;
			            dimu.chisq = 0.0;
				
			            dimu.znew = 0.0;
			            dimu.colin = 0.0;
			            dimu.dz1 = 0.0;
			            dimu.dz2 = 0.0;
			            dimu.xzee = 0.0;

			    
			    
			            // fill the dimuon struct
			    
			    

				        // filling the tree here gives me one tree entry
				        // for every dimuon candidate
			            t101->Fill();
			    
			   
			            h8->Fill(nfTrk);
			            h80->Fill(nfEmCell);

			    
			            etabyphi->Fill(dimu.eta,dimu.phi);
			            etabyphi->SetOption("lego");
				    
				    pt_distribution->Fill(dimu.pt);
			    
			    
			            if(Event.LM_north==0 && Event.Fastz==0 && Event.Phalo==0 && Event.Ahalo==0){
			      
			              //North off
			      
			              mass_gapN->Fill(dimu.mass);
			              eta_gapN->Fill(dimu.eta);
			      
			            }
			    
			            if(Event.LM_north==0 && Event.LM_south==1 && Event.Fastz==0 && Event.Phalo==0 && Event.Ahalo==0){
			      
			              //North off AND South on
			      
			              mass_NoffSon->Fill(dimu.mass);
			              eta_NoffSon->Fill(dimu.eta);
				      
			       
			              if(dimu.eta>-2 && dimu.eta<-1){eta_forneg_mass_NoffSon->Fill(dimu.mass);} 
			      
			              if(dimu.eta>-1 && dimu.eta<0){eta_centneg_mass_NoffSon->Fill(dimu.mass);}
			       
			              if(dimu.eta>0 && dimu.eta<1){eta_centpos_mass_NoffSon->Fill(dimu.mass);}
			       
			              if(dimu.eta>1 && dimu.eta<2){eta_forpos_mass_NoffSon->Fill(dimu.mass);} 
				      
				      
				      
				      if(dimu.y>-2 && dimu.y<-1){y_forneg_mass_NoffSon->Fill(dimu.mass);} 
			      
			              if(dimu.y>-1 && dimu.y<0){y_centneg_mass_NoffSon->Fill(dimu.mass);}
			       
			              if(dimu.y>0 && dimu.y<1){y_centpos_mass_NoffSon->Fill(dimu.mass);}
			       
			              if(dimu.y>1 && dimu.y<2){y_forpos_mass_NoffSon->Fill(dimu.mass);} 
				      
				      
				      
				      if(dimu.y>-0.6 && dimu.y<0){y_firstneg_mass_NoffSon->Fill(dimu.mass);} 
			      
			              if(dimu.y>-1.2 && dimu.y<-0.6){y_secondneg_mass_NoffSon->Fill(dimu.mass);}
			       
			              if(dimu.y>-1.8 && dimu.y<-1.2){y_thirdneg_mass_NoffSon->Fill(dimu.mass);}
				      
				      
				      if(dimu.y>0 && dimu.y<0.6){y_firstpos_mass_NoffSon->Fill(dimu.mass);} 
			      
			              if(dimu.y>0.6 && dimu.y<1.2){y_secondpos_mass_NoffSon->Fill(dimu.mass);}
			       
			              if(dimu.y>1.2 && dimu.y<1.8){y_thirdpos_mass_NoffSon->Fill(dimu.mass);}
				      
			       
			              if((dimu.y>-0.6 && dimu.y<0) || (dimu.y>0 && dimu.y<0.6)){y_first_mass_NoffSon->Fill(dimu.mass);}
				      
				      if((dimu.y>-1.2 && dimu.y<-0.6) || (dimu.y>0.6 && dimu.y<1.2)){y_second_mass_NoffSon->Fill(dimu.mass);}
				      
				      if((dimu.y>-1.8 && dimu.y<-1.2) || (dimu.y>1.2 && dimu.y<1.8)){y_third_mass_NoffSon->Fill(dimu.mass);} 
			      		
			      
			            }
			    
			            if(Event.LM_south==0 && Event.Fastz==0 && Event.Phalo==0 && Event.Ahalo==0){
			      
			              //South off
			      
			              mass_gapS->Fill(dimu.mass);
			              eta_gapS->Fill(dimu.eta);
			      
			            }
			    
			            if(Event.LM_south==0 && Event.LM_north==1 && Event.Fastz==0 && Event.Phalo==0 && Event.Ahalo==0){
			      
			              //South off AND North on
			      
			              mass_NonSoff->Fill(dimu.mass);
			              eta_NonSoff->Fill(dimu.eta);
			      
			              if(dimu.eta>-2 && dimu.eta<-1){eta_forneg_mass_NonSoff->Fill(dimu.mass);} 
			      
			              if(dimu.eta>-1 && dimu.eta<0){eta_centneg_mass_NonSoff->Fill(dimu.mass);}
			       
			              if(dimu.eta>0 && dimu.eta<1){eta_centpos_mass_NonSoff->Fill(dimu.mass);}
			       
			              if(dimu.eta>1 && dimu.eta<2){eta_forpos_mass_NonSoff->Fill(dimu.mass);} 
				      
				      
				      
				      if(dimu.y>-2 && dimu.y<-1){y_forneg_mass_NonSoff->Fill(dimu.mass);} 
			      
			              if(dimu.y>-1 && dimu.y<0){y_centneg_mass_NonSoff->Fill(dimu.mass);}
			       
			              if(dimu.y>0 && dimu.y<1){y_centpos_mass_NonSoff->Fill(dimu.mass);}
			       
			              if(dimu.y>1 && dimu.y<2){y_forpos_mass_NonSoff->Fill(dimu.mass);} 
				      
				      
				      
				      if(dimu.y>-0.6 && dimu.y<0){y_firstneg_mass_NonSoff->Fill(dimu.mass);} 
			      
			              if(dimu.y>-1.2 && dimu.y<-0.6){y_secondneg_mass_NonSoff->Fill(dimu.mass);}
			       
			              if(dimu.y>-1.8 && dimu.y<-1.2){y_thirdneg_mass_NonSoff->Fill(dimu.mass);}
				      
				      
				      if(dimu.y>0 && dimu.y<0.6){y_firstpos_mass_NonSoff->Fill(dimu.mass);} 
			      
			              if(dimu.y>0.6 && dimu.y<1.2){y_secondpos_mass_NonSoff->Fill(dimu.mass);}
			       
			              if(dimu.y>1.2 && dimu.y<1.8){y_thirdpos_mass_NonSoff->Fill(dimu.mass);}
				      
			       
			              if((dimu.y>-0.6 && dimu.y<0) || (dimu.y>0 && dimu.y<0.6)){y_first_mass_NonSoff->Fill(dimu.mass);}
				      
				      if((dimu.y>-1.2 && dimu.y<-0.6) || (dimu.y>0.6 && dimu.y<1.2)){y_second_mass_NonSoff->Fill(dimu.mass);}
				      
				      if((dimu.y>-1.8 && dimu.y<-1.2) || (dimu.y>1.2 && dimu.y<1.8)){y_third_mass_NonSoff->Fill(dimu.mass);}
			      
			      
			            }
			    
			            if(Event.LM_north==0 && Event.LM_south==0 && Event.Fastz==0 && Event.Phalo==0 && Event.Ahalo==0){
			      
			              //North and South off
			      
			              mass_gapNS->Fill(dimu.mass);
			              eta_gapNS->Fill(dimu.eta);
			      
			              if(dimu.eta>-2 && dimu.eta<-1){eta_forneg_mass_gapNS->Fill(dimu.mass);} 
			      
			              if(dimu.eta>-1 && dimu.eta<0){eta_centneg_mass_gapNS->Fill(dimu.mass);}
			       
			              if(dimu.eta>0 && dimu.eta<1){eta_centpos_mass_gapNS->Fill(dimu.mass);}
			       
			              if(dimu.eta>1 && dimu.eta<2){eta_forpos_mass_gapNS->Fill(dimu.mass);} 
				      
				      

				      if(dimu.y>-2 && dimu.y<-1){y_forneg_mass_gapNS->Fill(dimu.mass);} 
			      
			              if(dimu.y>-1 && dimu.y<0){y_centneg_mass_gapNS->Fill(dimu.mass);}
			       
			              if(dimu.y>0 && dimu.y<1){y_centpos_mass_gapNS->Fill(dimu.mass);}
			       
			              if(dimu.y>1 && dimu.y<2){y_forpos_mass_gapNS->Fill(dimu.mass);} 
				      
				      
				      
				      if(dimu.y>-0.6 && dimu.y<0){y_firstneg_mass_gapNS->Fill(dimu.mass);} 
			      
			              if(dimu.y>-1.2 && dimu.y<-0.6){y_secondneg_mass_gapNS->Fill(dimu.mass);}
			       
			              if(dimu.y>-1.8 && dimu.y<-1.2){y_thirdneg_mass_gapNS->Fill(dimu.mass);}
				      
				      
				      if(dimu.y>0 && dimu.y<0.6){y_firstpos_mass_gapNS->Fill(dimu.mass);} 
			      
			              if(dimu.y>0.6 && dimu.y<1.2){y_secondpos_mass_gapNS->Fill(dimu.mass);}
			       
			              if(dimu.y>1.2 && dimu.y<1.8){y_thirdpos_mass_gapNS->Fill(dimu.mass);}
				      
			       
			              if((dimu.y>-0.6 && dimu.y<0) || (dimu.y>0 && dimu.y<0.6)){y_first_mass_gapNS->Fill(dimu.mass);}
				      
				      if((dimu.y>-1.2 && dimu.y<-0.6) || (dimu.y>0.6 && dimu.y<1.2)){y_second_mass_gapNS->Fill(dimu.mass);}
				      
				      if((dimu.y>-1.8 && dimu.y<-1.2) || (dimu.y>1.2 && dimu.y<1.8)){y_third_mass_gapNS->Fill(dimu.mass);}
			      
			      
			            }
			    
			            //North and South on...
			    
			            //FastZ (in coincidence)
			            if(Event.LM_north==1 && Event.LM_south==1 && Event.Fastz==1){
			      
			              mass_fastZ->Fill(dimu.mass);
			      
			              if(dimu.eta>-2 && dimu.eta<-1){eta_forneg_mass_fastZ->Fill(dimu.mass);} 
			      
			              if(dimu.eta>-1 && dimu.eta<0){eta_centneg_mass_fastZ->Fill(dimu.mass);}
			       
			              if(dimu.eta>0 && dimu.eta<1){eta_centpos_mass_fastZ->Fill(dimu.mass);}
			       
			              if(dimu.eta>1 && dimu.eta<2){eta_forpos_mass_fastZ->Fill(dimu.mass);}
				      
				      
				      
				      if(dimu.y>-2 && dimu.y<-1){y_forneg_mass_fastZ->Fill(dimu.mass);} 
			      
			              if(dimu.y>-1 && dimu.y<0){y_centneg_mass_fastZ->Fill(dimu.mass);}
			       
			              if(dimu.y>0 && dimu.y<1){y_centpos_mass_fastZ->Fill(dimu.mass);}
			       
			              if(dimu.y>1 && dimu.y<2){y_forpos_mass_fastZ->Fill(dimu.mass);}  
				      
				      
				      
				      if(dimu.y>-0.6 && dimu.y<0){y_firstneg_mass_fastZ->Fill(dimu.mass);} 
			      
			              if(dimu.y>-1.2 && dimu.y<-0.6){y_secondneg_mass_fastZ->Fill(dimu.mass);}
			       
			              if(dimu.y>-1.8 && dimu.y<-1.2){y_thirdneg_mass_fastZ->Fill(dimu.mass);}
				      
				      
				      if(dimu.y>0 && dimu.y<0.6){y_firstpos_mass_fastZ->Fill(dimu.mass);} 
			      
			              if(dimu.y>0.6 && dimu.y<1.2){y_secondpos_mass_fastZ->Fill(dimu.mass);}
			       
			              if(dimu.y>1.2 && dimu.y<1.8){y_thirdpos_mass_fastZ->Fill(dimu.mass);}
				      
			       
			              if((dimu.y>-0.6 && dimu.y<0) || (dimu.y>0 && dimu.y<0.6)){y_first_mass_fastZ->Fill(dimu.mass);}
				      
				      if((dimu.y>-1.2 && dimu.y<-0.6) || (dimu.y>0.6 && dimu.y<1.2)){y_second_mass_fastZ->Fill(dimu.mass);}
				      
				      if((dimu.y>-1.8 && dimu.y<-1.2) || (dimu.y>1.2 && dimu.y<1.8)){y_third_mass_fastZ->Fill(dimu.mass);}
			      
			      
			            }
			   
			    
			            //Phalo (North fired early)
			            if(Event.LM_north==1 && Event.LM_south==1 && Event.Phalo==1){
			      
			              mass_Phalo->Fill(dimu.mass);
			      
			            }
			    
			            //Ahalo (South fired early)
			            if(Event.LM_north==1 && Event.LM_south==1 && Event.Ahalo==1){
			      
			              mass_Ahalo->Fill(dimu.mass);
			      
			            }


			          }
			        }
		              }
		            }
		          }
		        }
	              }
	            }
	          }
	        }
              }
            }
    
          }
        }
	
          
        Event.nmuon = nmuons*1.0;
        Event.ncand = ncand*1.0;
        Event.nDcand = nDcand*1.0;
	
    
        // this is filled once per event
        //t102->Fill();
    
        //    } // entry > 600
    }
  }
  
  // define your output file
  TFile file("test.root","recreate");
  
  
  // put the trees in the file
  t100->Write();
  t101->Write();
  //t102->Write();
  t103->Write();

  h1->Write();
  h2->Write();
  h3->Write();
  h4->Write();
  h5->Write();
  h6->Write();
  h7->Write();
  h8->Write();
  
  h80->Write();
  h81->Write();
  h82->Write();
  h83->Write();
  
  h85_lego->Write();
  h87_lego->Write();
  
  h85->Write();
  h85_log->Write();
  
  h85Noff->Write();
  h85Noff_lgy->Write();
  h85Noff_gapselection->Write();
  h85Soff->Write();
  h85Soff_lgy->Write();
  h85NSoff->Write();
  h85NSoff_lgy->Write();
  h85fastZ->Write();
  h85fastZ_lgy->Write();
  
  h87->Write();
  h87_log->Write();
  
  h87Noff->Write();
  h87Noff_lgy->Write();
  h87Soff->Write();
  h87Soff_lgy->Write();
  h87Soff_gapselection->Write();
  h87NSoff->Write();
  h87NSoff_lgy->Write();
  h87fastZ->Write();
  h87fastZ_lgy->Write();
  
  h88->Write();
  h88_log->Write();
  
  etabyphi->Write();
  pt_distribution->Write();
  
  mass_gapN->Write();
  mass_NoffSon->Write();
  eta_gapN->Write();
  eta_NoffSon->Write();
  
  mass_gapS->Write();
  mass_NonSoff->Write();
  eta_gapS->Write();
  eta_NonSoff->Write();
  
  mass_gapNS->Write();
  eta_gapNS->Write();
  
  mass_fastZ->Write();
  mass_Phalo->Write();
  mass_Ahalo->Write();
  
  eta_forneg_mass_NoffSon->Write();
  eta_centneg_mass_NoffSon->Write();
  eta_centpos_mass_NoffSon->Write();
  eta_forpos_mass_NoffSon->Write();
  
  eta_forneg_mass_NonSoff->Write();
  eta_centneg_mass_NonSoff->Write();
  eta_centpos_mass_NonSoff->Write();
  eta_forpos_mass_NonSoff->Write();
  
  eta_forneg_mass_gapNS->Write();
  eta_centneg_mass_gapNS->Write();
  eta_centpos_mass_gapNS->Write();
  eta_forpos_mass_gapNS->Write();
  
  eta_forneg_mass_fastZ->Write();
  eta_centneg_mass_fastZ->Write();
  eta_centpos_mass_fastZ->Write();
  eta_forpos_mass_fastZ->Write();
  
  y_forneg_mass_NoffSon->Write();
  y_centneg_mass_NoffSon->Write();
  y_centpos_mass_NoffSon->Write();
  y_forpos_mass_NoffSon->Write();
  
  y_forneg_mass_NonSoff->Write();
  y_centneg_mass_NonSoff->Write();
  y_centpos_mass_NonSoff->Write();
  y_forpos_mass_NonSoff->Write();
  
  y_forneg_mass_gapNS->Write();
  y_centneg_mass_gapNS->Write();
  y_centpos_mass_gapNS->Write();
  y_forpos_mass_gapNS->Write();
  
  y_forneg_mass_fastZ->Write();
  y_centneg_mass_fastZ->Write();
  y_centpos_mass_fastZ->Write();
  y_forpos_mass_fastZ->Write();
  
  
  y_firstneg_mass_NoffSon->Write();
  y_secondneg_mass_NoffSon->Write();
  y_thirdneg_mass_NoffSon->Write();

  y_firstpos_mass_NoffSon->Write();
  y_secondpos_mass_NoffSon->Write();
  y_thirdpos_mass_NoffSon->Write();
  
  y_firstneg_mass_NonSoff->Write();
  y_secondneg_mass_NonSoff->Write();
  y_thirdneg_mass_NonSoff->Write();
  
  y_firstpos_mass_NonSoff->Write();
  y_secondpos_mass_NonSoff->Write();
  y_thirdpos_mass_NonSoff->Write();
  
  y_firstneg_mass_gapNS->Write();
  y_secondneg_mass_gapNS->Write();
  y_thirdneg_mass_gapNS->Write();
  
  y_firstpos_mass_gapNS->Write();
  y_secondpos_mass_gapNS->Write();
  y_thirdpos_mass_gapNS->Write();
  
  y_firstneg_mass_fastZ->Write();
  y_secondneg_mass_fastZ->Write();
  y_thirdneg_mass_fastZ->Write();
  
  y_firstpos_mass_fastZ->Write();
  y_secondpos_mass_fastZ->Write();
  y_thirdpos_mass_fastZ->Write();

  y_first_mass_NoffSon->Write();
  y_second_mass_NoffSon->Write();
  y_third_mass_NoffSon->Write();

  y_first_mass_NonSoff->Write();
  y_second_mass_NonSoff->Write();
  y_third_mass_NonSoff->Write();
  
  y_first_mass_gapNS->Write();
  y_second_mass_gapNS->Write();
  y_third_mass_gapNS->Write();
  
  y_first_mass_fastZ->Write();
  y_second_mass_fastZ->Write();
  y_third_mass_fastZ->Write();

  
  
  // close the file
  file.Close();
  
  
} // TMBTree_bu::Loop