#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, 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, nmuon, ncand, nDcand;
    Int_t  Fastz, LM_north, LM_south, Ahalo, Phalo;
  } my_event;

  // 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;

  

  // declare structures you will use
  static my_muon All_muons;
  static my_muon muon2;
  static my_dimuon dimu;
  static my_event Event;
  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
  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
  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: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");
  
  //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);
  
  TH1F* h80 = new TH1F("h80","n emcells",200,0.,300.0);
  TH1F* h81 = new TH1F("h81","energia",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);
  TH1F* h84 = new TH1F("h84","energia",200,0.,20.0);
  TH1F* h85 = new TH1F("h85","energia",200,0.,!00);
  TH1F* h86 = new TH1F("h86","energia",200,0.,20.0);
  TH1F* h87 = new TH1F("h87","energia",200,0.,100);
  TH1F* h88 = new TH1F("h88","energia",200,0.,100);

  
  /// Adding my own histograms ///
  TH2F* h9 = new TH2F("h9","eta by phi",100,-6,6,100,0,7);
  TH2F* h10 = new TH2F("h10","pt by phi",100,0,80,100,0,7);
  TH2F* h11 = new TH2F("h11","pt by eta",100,0,80,100,-6,6);
  TH1F* h12 = new TH1F("h12","muon1_eta",100,-6,6);
  TH1F* h13 = new TH1F("h13","muon2_eta",100,-6,6);
  TH1F* h14 = new TH1F("h14","dimu_eta",100,-6,6);
  TH2F* h15 = new TH2F("h15","mass by eta",100,-6,6,100,1,5);
  
  TH1F* h16 = new TH1F("h16","mass with GapN",200,1,5);
  TH1F* h161 = new TH1F("h161","mass with LM_North off AND LM_South on",200,1,5);
  TH2F* h17 = new TH2F("h17","eta by phi with GapN",100,-6,6,100,0,7);
  TH2F* h171 = new TH2F("h171","eta by phi with LM_North off AND LM_South on",100,-6,6,100,0,7);
  TH1F* h18 = new TH1F("h18","mass with GapS",200,1,5);
  TH1F* h181 = new TH1F("h181","mass with LM_South off AND LM_North on",200,1,5);
  TH2F* h19 = new TH2F("h19","eta by phi with GapS",100,-6,6,100,0,7);
  TH2F* h191 = new TH2F("h191","eta by phi with LM_South off AND LM_North on",100,-6,6,100,0,7);
  TH1F* h20 = new TH1F("h20","mass with GapNS",200,1,5);
  TH2F* h21 = new TH2F("h21","eta by phi with GapNS",100,-6,6,100,0,7);
  TH1F* h22 = new TH1F("h22","mass and the FastZ",200,1,5);
  TH1F* h23 = new TH1F("h23","mass and the Phalo",200,1,5);
  TH1F* h24 = new TH1F("h24","mass and the Ahalo",200,1,5);
  
  TH1F* h25 = new TH1F("h25","pt mu1",100,0,10);
  TH1F* h26 = new TH1F("h26","pt mu2",100,0,10);
  TH1F* h27 = new TH1F("h27","pt dimu",100,0,10);
  
  
  // 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;}
  
  
  // 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;


    // 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();

    
    // 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;
    
    //Sum EM cells?
    Double_t sumn_E = 0.0;
    Double_t sums_E = 0.0;
    Double_t sum;
    
    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];}

    }			  
    
    //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;
    
    
    //loop over the calorimeter cells
    for(Int_t a=0; a<nfEmCell; ++a){
      
        TMBEmCells* EmCell = ( TMBEmCells* ) fEmCells->At(a);
	
        energia = EmCell->E();
	i_eta = EmCell->ieta();
	i_phi = EmCell->iphi();
	i_layer = EmCell->ilyr();
	
        //cout << energia << " " << i_eta << endl;
	
        h81->Fill(energia);
	h82->Fill(i_eta,energia);
	h82->SetOption("lego");
	//h83->Fill(i_layer);
	if(energia>0.1){
	  h83->Fill(i_layer);
	  
	  //north
	  if(i_eta>=-35 && i_eta<=-28){
	    h84->Fill(energia);
	    //cout << energia << " " << sumn_E << endl;
	    sumn_E = sumn_E + energia;
	  }
	  
	  //south
	  if(i_eta>=28 && i_eta<=35){
	    h86->Fill(energia);
	    sums_E = sums_E + energia;
	  }
	  
	}  

    }
    h85->Fill(sumn_E);
    h87->Fill(sums_E);
    sum = sumn_E + sums_E;
    h88->Fill(sum);

    //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();
	      h25->Fill(All_muons.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();
			  h26->Fill(muon2.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.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++;


			    //
				// 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);
			    h9->Fill(dimu.eta,dimu.phi);
			    h9->SetOption("lego");
			    h10->Fill(dimu.pt,dimu.phi);
			    h11->Fill(dimu.pt,dimu.eta);
			    
			    h12->Fill(All_muons.eta);
			    h12->SetLineColor(kBlue);
			    h13->Fill(muon2.eta);
			    h13->SetLineColor(kGreen);
			    h14->Fill(dimu.eta);
			    h14->SetLineColor(kMagenta);
			    
			    h15->Fill(dimu.eta,dimu.mass);
			    
			    h27->Fill(dimu.pt);
			    
			    
			    if(Event.LM_north==0 && Event.Fastz==0 && Event.Phalo==0 && Event.Ahalo==0){
			      
			      //North off
			      h16->Fill(dimu.mass);
			      
			      h17->Fill(dimu.eta,dimu.phi);
			      
			    }
			    
			    if(Event.LM_north==0 && Event.LM_south==1 && Event.Fastz==0 && Event.Phalo==0 && Event.Ahalo==0){
			      
			      //North off AND South on
			      h161->Fill(dimu.mass);
			      
			      h171->Fill(dimu.eta,dimu.phi);
			      
			    }
			    
			    if(Event.LM_south==0 && Event.Fastz==0 && Event.Phalo==0 && Event.Ahalo==0){
			      
			      //South off
			      h18->Fill(dimu.mass);
			      
			      h19->Fill(dimu.eta,dimu.phi);
			      
			    }
			    
			    if(Event.LM_south==0 && Event.LM_north==1 && Event.Fastz==0 && Event.Phalo==0 && Event.Ahalo==0){
			      
			      //South off AND North on
			      h181->Fill(dimu.mass);
			      
			      h191->Fill(dimu.eta,dimu.phi);
			      
			    }
			    
			    if(Event.LM_north==0 && Event.LM_south==0 && Event.Fastz==0 && Event.Phalo==0 && Event.Ahalo==0){
			      
			      //North and South off
			      h20->Fill(dimu.mass);
			      
			      h21->Fill(dimu.eta,dimu.phi);
			      
			    }
			    
			    //North and South on...
			    
			    //FastZ (in coincidence)
			    if(Event.LM_north==1 && Event.LM_south==1 && Event.Fastz==1){
			      
			      h22->Fill(dimu.mass);
			      
			    }
			    
			    //Phalo (North fired early)
			    if(Event.LM_north==1 && Event.LM_south==1 && Event.Phalo==1){
			      
			      h23->Fill(dimu.mass);
			      
			    }
			    
			    //Ahalo (South fired early)
			    if(Event.LM_north==1 && Event.LM_south==1 && Event.Ahalo==1){
			      
			      h24->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();
  

  h1->Write();
  h2->Write();
  h3->Write();
  h4->Write();
  h5->Write();
  h6->Write();
  h7->Write();
  h8->Write();
  
  h80->Write();
  h81->Write();
  h82->Write();
  h83->Write();
  h84->Write();
  h85->Write();
  h86->Write();
  h87->Write();
  h88->Write();
  
  h9->Write();
  h10->Write();
  h11->Write();
  h12->Write();
  h13->Write();
  h14->Write();
  h15->Write();
  h16->Write();
  h161->Write();
  h17->Write();
  h171->Write();
  h18->Write();
  h181->Write();
  h19->Write();
  h191->Write();
  h20->Write();
  h21->Write();
  h22->Write();
  h23->Write();
  h24->Write();
  h25->Write();
  h26->Write();
  h27->Write();
  
  //TF1 *f1=new TF1("f1","(1/[0]*sqrt(2*3.1415))*exp(-(pow(x-[1],2))/(2*pow([2],2)))",2.5,3.5);
  //f1->SetParameters(1,2,3);
  //f1->SetLineColor(kBlue);
  //h30->Fit("f1","R");
  
  // close the file
  file.Close();

} // TMBTree_bu::Loop