FitNtuple/src/FitNtuple.cc

// -*- C++ -*-
//
// Package:    FitNtuple
// Class:      FitNtuple
// 
/**\class FitNtuple FitNtuple.cc AnalysisTools/FitNtuple/src/FitNtuple.cc

 Description: [one line class summary]

 Implementation:
     [Notes on implementation]
*/
//
// Original Author:  James Pivarski
//         Created:  Fri Nov 19 17:04:33 CST 2010
// $Id: FitNtuple.cc,v 1.13 2010/12/15 20:55:09 pivarski Exp $
//
//


// system include files
#include <memory>

// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EDAnalyzer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/InputTag.h"

#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
#include "DataFormats/Candidate/interface/Candidate.h"
#include "AnalysisDataFormats/MuJetAnalysis/interface/MultiMuon.h"
#include "DataFormats/PatCandidates/interface/Muon.h"
#include "DataFormats/PatCandidates/interface/TriggerEvent.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/VertexReco/interface/VertexFwd.h"
#include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
#include "DataFormats/MuonReco/interface/Muon.h"
#include "DataFormats/MuonReco/interface/MuonChamberMatch.h"
#include "DataFormats/MuonReco/interface/MuonSegmentMatch.h"

#include "CommonTools/UtilAlgos/interface/TFileService.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "TTree.h"
#include "TRandom3.h"

//
// class declaration
//

class FitNtuple : public edm::EDAnalyzer {
   public:
      explicit FitNtuple(const edm::ParameterSet&);
      ~FitNtuple();


   private:
      virtual void beginJob() ;
      virtual void analyze(const edm::Event&, const edm::EventSetup&);
      virtual void endJob() ;

      // ----------member data ---------------------------

      // input parameters
      bool m_getQscale;
      bool m_getAlternating;
      bool m_fillGenLevel;

      TRandom3 m_trandom3;

      // all ntuples need these variables
      Float_t m_qscale;  // also known as "pt-hat" ($\hat{p}_T$)
      Int_t m_genLevel;  // gen-level information about all mu-jets
                         //                     1 = b-jets
                         //                     2 = resonances
                         //                     4 = any decays-in-flight
                         //                     8 = any non-muons
      Int_t m_run;       // run and event number so that we can look at
      Int_t m_event;     //     interesting cases in the event display
      Int_t m_trigger;   // satisfied trigger?  0 = failed all
                         //                     1 = passed HLT_Mu5
                         //                     2 = passed HLT_Mu9
                         //                     4 = passed HLT_Mu11
                         //                     8 = passed HLT_Mu15_v1
                         //                    16 = passed HLT_Mu20_NoVertex
      Float_t m_muon1pt;
      Float_t m_muon1eta;
      Float_t m_muon2pt;
      Float_t m_muon2eta;
      Float_t m_muon3pt;
      Float_t m_muon3eta;
      Float_t m_muon4pt;
      Float_t m_muon4eta;
      Float_t m_muontrigpt;
      Float_t m_muontrigeta;

      // for scaling to the Z peak
      TTree *m_alldimuons;
      Float_t m_alldimuons_mass;
      Float_t m_alldimuons_iso;

      // control sample "lowdimuon" (single mu-jet with two muons and low pT)
      TTree *m_lowdimuon;
      Int_t m_lowdimuon_containstrig;
      Float_t m_lowdimuon_genmass;
      Float_t m_lowdimuon_mass;
      Float_t m_lowdimuon_pt;
      Float_t m_lowdimuon_eta;
      Float_t m_lowdimuon_phi;
      Float_t m_lowdimuon_dr;
      Float_t m_lowdimuon_pluspx;
      Float_t m_lowdimuon_pluspy;
      Float_t m_lowdimuon_pluspz;
      Float_t m_lowdimuon_minuspx;
      Float_t m_lowdimuon_minuspy;
      Float_t m_lowdimuon_minuspz;
      Float_t m_lowdimuon_vprob;
      Float_t m_lowdimuon_vx;
      Float_t m_lowdimuon_vy;
      Float_t m_lowdimuon_vz;
      Float_t m_lowdimuon_iso;
      Float_t m_lowdimuon_compluspx;
      Float_t m_lowdimuon_compluspy;
      Float_t m_lowdimuon_compluspz;
      Float_t m_lowdimuon_comrotpluspx;
      Float_t m_lowdimuon_comrotpluspy;
      Float_t m_lowdimuon_comrotpluspz;
      Int_t m_lowdimuon_plusmatches;
      Float_t m_lowdimuon_plusst1x;
      Float_t m_lowdimuon_plusst1xsig;
      Float_t m_lowdimuon_plusst1dxdz;
      Float_t m_lowdimuon_plusst1dxdzsig;
      Int_t m_lowdimuon_plushits;
      Float_t m_lowdimuon_plusnormchi2;
      Int_t m_lowdimuon_minusmatches;
      Float_t m_lowdimuon_minusst1x;
      Float_t m_lowdimuon_minusst1xsig;
      Float_t m_lowdimuon_minusst1dxdz;
      Float_t m_lowdimuon_minusst1dxdzsig;
      Int_t m_lowdimuon_minushits;
      Float_t m_lowdimuon_minusnormchi2;

      // signal region "highdimuon" (single mu-jet with two muons and high pT)
      TTree *m_highdimuon;

      // control sample "mujetplustracks" (single mu-jet with 2-3 muons and 1-2 random, nearby tracks)
      TTree *m_mujetplustracks;
      Int_t m_mujetplustracks_containstrig;
      Int_t m_mujetplustracks_muons;
      Int_t m_mujetplustracks_fakes;
      Int_t m_mujetplustracks_charge;
      Float_t m_mujetplustracks_mass;
      Float_t m_mujetplustracks_massa;
      Float_t m_mujetplustracks_massb;
      Float_t m_mujetplustracks_pt;
      Float_t m_mujetplustracks_eta;
      Float_t m_mujetplustracks_phi;
      Float_t m_mujetplustracks_dr;
      Float_t m_mujetplustracks_vprob;
      Float_t m_mujetplustracks_lxy;
      Float_t m_mujetplustracks_lxyz;
      Float_t m_mujetplustracks_iso;
      Float_t m_mujetplustracks_mupt1;
      Float_t m_mujetplustracks_mupt2;
      Float_t m_mujetplustracks_mupt3;
      Float_t m_mujetplustracks_fakept1;
      Float_t m_mujetplustracks_fakept2;

      // signal region "quadmuon" (single mu-jet with four or more muons)
      TTree *m_quadmuon;
      Int_t m_quadmuon_containstrig;
      Int_t m_quadmuon_muons;
      Int_t m_quadmuon_orphans;
      Int_t m_quadmuon_charge;
      Float_t m_quadmuon_mass;
      Float_t m_quadmuon_massa;
      Float_t m_quadmuon_massb;
      Float_t m_quadmuon_massc;
      Float_t m_quadmuon_massd;
      Float_t m_quadmuon_masse;
      Float_t m_quadmuon_massf;
      Float_t m_quadmuon_pt;
      Float_t m_quadmuon_eta;
      Float_t m_quadmuon_phi;
      Float_t m_quadmuon_dr;
      Float_t m_quadmuon_vchi2;
      Int_t m_quadmuon_vndof;
      Float_t m_quadmuon_vprob;
      Float_t m_quadmuon_lxy;
      Float_t m_quadmuon_lxyz;
      Float_t m_quadmuon_iso;

      // control sample "dimuorphan" (single, two-muon mu-jet and a single muon)
      TTree *m_dimuorphan;
      Float_t m_dimuorphan_deltaphi;
      Float_t m_dimuorphan_orphanpt;
      Float_t m_dimuorphan_orphaneta;
      Float_t m_dimuorphan_orphanphi;
      Int_t m_dimuorphan_orphanisglobal;
      Int_t m_dimuorphan_orphanmatches;
      Int_t m_dimuorphan_orphanhits;
      Float_t m_dimuorphan_orphanchi2;
      Float_t m_dimuorphan_orphanst1x;
      Float_t m_dimuorphan_orphanst1xsig;
      Float_t m_dimuorphan_orphanst1dxdz;
      Float_t m_dimuorphan_orphanst1dxdzsig;
      Int_t m_dimuorphan_containstrig;
      Float_t m_dimuorphan_mass;
      Float_t m_dimuorphan_pt;
      Float_t m_dimuorphan_eta;
      Float_t m_dimuorphan_phi;
      Float_t m_dimuorphan_dr;
      Float_t m_dimuorphan_vprob;
      Float_t m_dimuorphan_lxy;
      Float_t m_dimuorphan_lxyz;
      Float_t m_dimuorphan_iso;

      // signal region "dimudimu" (two mu-jets with two muons each)
      TTree *m_dimudimu;
      Int_t m_dimudimu_orphans;
      Float_t m_dimudimu_wholemass;
      Float_t m_dimudimu_wholept;
      Float_t m_dimudimu_deltaphi;
      Int_t m_dimudimu_containstrig;
      Float_t m_dimudimu_massC;
      Float_t m_dimudimu_ptC;
      Float_t m_dimudimu_etaC;
      Float_t m_dimudimu_phiC;
      Float_t m_dimudimu_drC;
      Float_t m_dimudimu_vprobC;
      Float_t m_dimudimu_lxyC;
      Float_t m_dimudimu_lxyzC;
      Float_t m_dimudimu_isoC;
      Float_t m_dimudimu_massF;
      Float_t m_dimudimu_ptF;
      Float_t m_dimudimu_etaF;
      Float_t m_dimudimu_phiF;
      Float_t m_dimudimu_drF;
      Float_t m_dimudimu_vprobF;
      Float_t m_dimudimu_lxyF;
      Float_t m_dimudimu_lxyzF;
      Float_t m_dimudimu_isoF;

      // signal region "twomujets"
      TTree *m_twomujets;
      Int_t m_twomujets_muonsC;
      Int_t m_twomujets_muonsF;
      Float_t m_twomujets_massCa;
      Float_t m_twomujets_massCb;
      Float_t m_twomujets_massCc;
      Float_t m_twomujets_massFa;
      Float_t m_twomujets_massFb;
      Float_t m_twomujets_massFc;

      // signal region "moremujets"
      TTree *m_moremujets;
      Int_t m_moremujets_mujets;
      Int_t m_moremujets_orphans;
      Float_t m_moremujets_wholemass;
      Float_t m_moremujets_wholept;
      Float_t m_moremujets_massa;
      Float_t m_moremujets_massb;
      Float_t m_moremujets_massc;
      Float_t m_moremujets_massd;
      Float_t m_moremujets_masse;
      Float_t m_moremujets_massf;
      Int_t m_moremujets_muons1;
      Int_t m_moremujets_containstrig1;
      Float_t m_moremujets_mass1;
      Float_t m_moremujets_pt1;
      Float_t m_moremujets_iso1;
      Int_t m_moremujets_muons2;
      Int_t m_moremujets_containstrig2;
      Float_t m_moremujets_mass2;
      Float_t m_moremujets_pt2;
      Float_t m_moremujets_iso2;
      Int_t m_moremujets_muons3;
      Int_t m_moremujets_containstrig3;
      Float_t m_moremujets_mass3;
      Float_t m_moremujets_pt3;
      Float_t m_moremujets_iso3;
      Int_t m_moremujets_muons4;
      Int_t m_moremujets_containstrig4;
      Float_t m_moremujets_mass4;
      Float_t m_moremujets_pt4;
      Float_t m_moremujets_iso4;
      Int_t m_moremujets_muons5;
      Int_t m_moremujets_containstrig5;
      Float_t m_moremujets_mass5;
      Float_t m_moremujets_pt5;
      Float_t m_moremujets_iso5;
      Int_t m_moremujets_muons6;
      Int_t m_moremujets_containstrig6;
      Float_t m_moremujets_mass6;
      Float_t m_moremujets_pt6;
      Float_t m_moremujets_iso6;

      // background-enriched and signal for displaced muons
      TTree *m_displacedFake;
      TTree *m_displaced;
      Int_t m_displaced_mujets;
      Int_t m_displaced_orphans;
      Int_t m_displaced_muons;
      Float_t m_displaced_mass;
      Float_t m_displaced_pt;
      Float_t m_displaced_eta;
      Float_t m_displaced_phi;
      Float_t m_displaced_vx;
      Float_t m_displaced_vy;
      Float_t m_displaced_vz;
      Float_t m_displaced_zclosest;
      Float_t m_displaced_zcloseness;
      Float_t m_displaced_lxy;
      Float_t m_displaced_lxyz;
      Float_t m_displaced_iso;
      Float_t m_displaced_vprob;
      Float_t m_displaced_mupt1;
      Float_t m_displaced_mupt2;
      Int_t m_displaced_plusmatches;
      Float_t m_displaced_plusst1x;
      Float_t m_displaced_plusst1xsig;
      Float_t m_displaced_plusst1dxdz;
      Float_t m_displaced_plusst1dxdzsig;
      Int_t m_displaced_plushits;
      Float_t m_displaced_plusnormchi2;
      Int_t m_displaced_minusmatches;
      Float_t m_displaced_minusst1x;
      Float_t m_displaced_minusst1xsig;
      Float_t m_displaced_minusst1dxdz;
      Float_t m_displaced_minusst1dxdzsig;
      Int_t m_displaced_minushits;
      Float_t m_displaced_minusnormchi2;
      Int_t m_displacedFake_fakes;
};

//
// constants, enums and typedefs
//

//
// static data member definitions
//

//
// constructors and destructor
//
FitNtuple::FitNtuple(const edm::ParameterSet& iConfig)
   : m_getQscale(iConfig.getParameter<bool>("getQscale"))
   , m_getAlternating(iConfig.getParameter<bool>("getAlternating"))
   , m_fillGenLevel(iConfig.getParameter<bool>("fillGenLevel"))
{
   //now do what ever initialization is needed
   edm::Service<TFileService> tFile;

   m_alldimuons = tFile->make<TTree>("alldimuons", "alldimuons");
   m_alldimuons->Branch("trigger", &m_trigger, "trigger/I");
   m_alldimuons->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
   m_alldimuons->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
   m_alldimuons->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
   m_alldimuons->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
   m_alldimuons->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
   m_alldimuons->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
   m_alldimuons->Branch("mass", &m_alldimuons_mass, "mass/F");
   m_alldimuons->Branch("iso", &m_alldimuons_iso, "iso/F");

   m_lowdimuon = tFile->make<TTree>("lowdimuon", "lowdimuon");
   m_lowdimuon->Branch("qscale", &m_qscale, "qscale/F");
   m_lowdimuon->Branch("genLevel", &m_genLevel, "genLevel/I");
   m_lowdimuon->Branch("run", &m_run, "run/I");
   m_lowdimuon->Branch("event", &m_event, "event/I");
   m_lowdimuon->Branch("trigger", &m_trigger, "trigger/I");
   m_lowdimuon->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
   m_lowdimuon->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
   m_lowdimuon->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
   m_lowdimuon->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
   m_lowdimuon->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
   m_lowdimuon->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
   m_lowdimuon->Branch("containstrig", &m_lowdimuon_containstrig, "containstrig/I");
   m_lowdimuon->Branch("genmass", &m_lowdimuon_genmass, "genmass/F");
   m_lowdimuon->Branch("mass", &m_lowdimuon_mass, "mass/F");
   m_lowdimuon->Branch("pt", &m_lowdimuon_pt, "pt/F");
   m_lowdimuon->Branch("eta", &m_lowdimuon_eta, "eta/F");
   m_lowdimuon->Branch("phi", &m_lowdimuon_phi, "phi/F");
   m_lowdimuon->Branch("dr", &m_lowdimuon_dr, "dr/F");
   m_lowdimuon->Branch("pluspx", &m_lowdimuon_pluspx, "pluspx/F");
   m_lowdimuon->Branch("pluspy", &m_lowdimuon_pluspy, "pluspy/F");
   m_lowdimuon->Branch("pluspz", &m_lowdimuon_pluspz, "pluspz/F");
   m_lowdimuon->Branch("minuspx", &m_lowdimuon_minuspx, "minuspx/F");
   m_lowdimuon->Branch("minuspy", &m_lowdimuon_minuspy, "minuspy/F");
   m_lowdimuon->Branch("minuspz", &m_lowdimuon_minuspz, "minuspz/F");
   m_lowdimuon->Branch("vprob", &m_lowdimuon_vprob, "vprob/F");
   m_lowdimuon->Branch("vx", &m_lowdimuon_vx, "vx/F");
   m_lowdimuon->Branch("vy", &m_lowdimuon_vy, "vy/F");
   m_lowdimuon->Branch("vz", &m_lowdimuon_vz, "vz/F");
   m_lowdimuon->Branch("iso", &m_lowdimuon_iso, "iso/F");
   m_lowdimuon->Branch("compluspx", &m_lowdimuon_compluspx, "compluspx/F");
   m_lowdimuon->Branch("compluspy", &m_lowdimuon_compluspy, "compluspy/F");
   m_lowdimuon->Branch("compluspz", &m_lowdimuon_compluspz, "compluspz/F");
   m_lowdimuon->Branch("comrotpluspx", &m_lowdimuon_comrotpluspx, "comrotpluspx/F");
   m_lowdimuon->Branch("comrotpluspy", &m_lowdimuon_comrotpluspy, "comrotpluspy/F");
   m_lowdimuon->Branch("comrotpluspz", &m_lowdimuon_comrotpluspz, "comrotpluspz/F");
   m_lowdimuon->Branch("plusmatches", &m_lowdimuon_plusmatches, "plusmatches/I");
   m_lowdimuon->Branch("plusst1x", &m_lowdimuon_plusst1x, "plusst1x/F");
   m_lowdimuon->Branch("plusst1xsig", &m_lowdimuon_plusst1xsig, "plusst1xsig/F");
   m_lowdimuon->Branch("plusst1dxdz", &m_lowdimuon_plusst1dxdz, "plusst1dxdz/F");
   m_lowdimuon->Branch("plusst1dxdzsig", &m_lowdimuon_plusst1dxdzsig, "plusst1dxdzsig/F");
   m_lowdimuon->Branch("plushits", &m_lowdimuon_plushits, "plushits/I");
   m_lowdimuon->Branch("plusnormchi2", &m_lowdimuon_plusnormchi2, "plusnormchi2/F");
   m_lowdimuon->Branch("minusmatches", &m_lowdimuon_minusmatches, "minusmatches/I");
   m_lowdimuon->Branch("minusst1x", &m_lowdimuon_minusst1x, "minusst1x/F");
   m_lowdimuon->Branch("minusst1xsig", &m_lowdimuon_minusst1xsig, "minusst1xsig/F");
   m_lowdimuon->Branch("minusst1dxdz", &m_lowdimuon_minusst1dxdz, "minusst1dxdz/F");
   m_lowdimuon->Branch("minusst1dxdzsig", &m_lowdimuon_minusst1dxdzsig, "minusst1dxdzsig/F");
   m_lowdimuon->Branch("minushits", &m_lowdimuon_minushits, "minushits/I");
   m_lowdimuon->Branch("minusnormchi2", &m_lowdimuon_minusnormchi2, "minusnormchi2/F");

   m_highdimuon = tFile->make<TTree>("highdimuon", "highdimuon");
   m_highdimuon->Branch("qscale", &m_qscale, "qscale/F");
   m_highdimuon->Branch("genLevel", &m_genLevel, "genLevel/I");
   m_highdimuon->Branch("run", &m_run, "run/I");
   m_highdimuon->Branch("event", &m_event, "event/I");
   m_highdimuon->Branch("trigger", &m_trigger, "trigger/I");
   m_highdimuon->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
   m_highdimuon->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
   m_highdimuon->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
   m_highdimuon->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
   m_highdimuon->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
   m_highdimuon->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
   m_highdimuon->Branch("containstrig", &m_lowdimuon_containstrig, "containstrig/I");
   m_highdimuon->Branch("genmass", &m_lowdimuon_genmass, "genmass/F");
   m_highdimuon->Branch("mass", &m_lowdimuon_mass, "mass/F");
   m_highdimuon->Branch("pt", &m_lowdimuon_pt, "pt/F");
   m_highdimuon->Branch("eta", &m_lowdimuon_eta, "eta/F");
   m_highdimuon->Branch("phi", &m_lowdimuon_phi, "phi/F");
   m_highdimuon->Branch("dr", &m_lowdimuon_dr, "dr/F");
   m_highdimuon->Branch("pluspx", &m_lowdimuon_pluspx, "pluspx/F");
   m_highdimuon->Branch("pluspy", &m_lowdimuon_pluspy, "pluspy/F");
   m_highdimuon->Branch("pluspz", &m_lowdimuon_pluspz, "pluspz/F");
   m_highdimuon->Branch("minuspx", &m_lowdimuon_minuspx, "minuspx/F");
   m_highdimuon->Branch("minuspy", &m_lowdimuon_minuspy, "minuspy/F");
   m_highdimuon->Branch("minuspz", &m_lowdimuon_minuspz, "minuspz/F");
   m_highdimuon->Branch("vprob", &m_lowdimuon_vprob, "vprob/F");
   m_highdimuon->Branch("vx", &m_lowdimuon_vx, "vx/F");
   m_highdimuon->Branch("vy", &m_lowdimuon_vy, "vy/F");
   m_highdimuon->Branch("vz", &m_lowdimuon_vz, "vz/F");
   m_highdimuon->Branch("iso", &m_lowdimuon_iso, "iso/F");
   m_highdimuon->Branch("compluspx", &m_lowdimuon_compluspx, "compluspx/F");
   m_highdimuon->Branch("compluspy", &m_lowdimuon_compluspy, "compluspy/F");
   m_highdimuon->Branch("compluspz", &m_lowdimuon_compluspz, "compluspz/F");
   m_highdimuon->Branch("comrotpluspx", &m_lowdimuon_comrotpluspx, "comrotpluspx/F");
   m_highdimuon->Branch("comrotpluspy", &m_lowdimuon_comrotpluspy, "comrotpluspy/F");
   m_highdimuon->Branch("comrotpluspz", &m_lowdimuon_comrotpluspz, "comrotpluspz/F");
   m_highdimuon->Branch("plusmatches", &m_lowdimuon_plusmatches, "plusmatches/I");
   m_highdimuon->Branch("plusst1x", &m_lowdimuon_plusst1x, "plusst1x/F");
   m_highdimuon->Branch("plusst1xsig", &m_lowdimuon_plusst1xsig, "plusst1xsig/F");
   m_highdimuon->Branch("plusst1dxdz", &m_lowdimuon_plusst1dxdz, "plusst1dxdz/F");
   m_highdimuon->Branch("plusst1dxdzsig", &m_lowdimuon_plusst1dxdzsig, "plusst1dxdzsig/F");
   m_highdimuon->Branch("plushits", &m_lowdimuon_plushits, "plushits/I");
   m_highdimuon->Branch("plusnormchi2", &m_lowdimuon_plusnormchi2, "plusnormchi2/F");
   m_highdimuon->Branch("minusmatches", &m_lowdimuon_minusmatches, "minusmatches/I");
   m_highdimuon->Branch("minusst1x", &m_lowdimuon_minusst1x, "minusst1x/F");
   m_highdimuon->Branch("minusst1xsig", &m_lowdimuon_minusst1xsig, "minusst1xsig/F");
   m_highdimuon->Branch("minusst1dxdz", &m_lowdimuon_minusst1dxdz, "minusst1dxdz/F");
   m_highdimuon->Branch("minusst1dxdzsig", &m_lowdimuon_minusst1dxdzsig, "minusst1dxdzsig/F");
   m_highdimuon->Branch("minushits", &m_lowdimuon_minushits, "minushits/I");
   m_highdimuon->Branch("minusnormchi2", &m_lowdimuon_minusnormchi2, "minusnormchi2/F");

   m_mujetplustracks = tFile->make<TTree>("mujetplustracks", "mujetplustracks");
   m_mujetplustracks->Branch("qscale", &m_qscale, "qscale/F");
   m_mujetplustracks->Branch("run", &m_run, "run/I");
   m_mujetplustracks->Branch("event", &m_event, "event/I");
   m_mujetplustracks->Branch("trigger", &m_trigger, "trigger/I");
   m_mujetplustracks->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
   m_mujetplustracks->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
   m_mujetplustracks->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
   m_mujetplustracks->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
   m_mujetplustracks->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
   m_mujetplustracks->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
   m_mujetplustracks->Branch("containstrig", &m_mujetplustracks_containstrig, "containstrig/I");
   m_mujetplustracks->Branch("muons", &m_mujetplustracks_muons, "muons/I");
   m_mujetplustracks->Branch("fakes", &m_mujetplustracks_fakes, "fakes/I");
   m_mujetplustracks->Branch("charge", &m_mujetplustracks_charge, "charge/I");
   m_mujetplustracks->Branch("mass", &m_mujetplustracks_mass, "mass/F");
   m_mujetplustracks->Branch("massa", &m_mujetplustracks_massa, "massa/F");
   m_mujetplustracks->Branch("massb", &m_mujetplustracks_massb, "massb/F");
   m_mujetplustracks->Branch("pt", &m_mujetplustracks_pt, "pt/F");
   m_mujetplustracks->Branch("eta", &m_mujetplustracks_eta, "eta/F");
   m_mujetplustracks->Branch("phi", &m_mujetplustracks_phi, "phi/F");
   m_mujetplustracks->Branch("dr", &m_mujetplustracks_dr, "dr/F");
   m_mujetplustracks->Branch("vprob", &m_mujetplustracks_vprob, "vprob/F");
   m_mujetplustracks->Branch("lxy", &m_mujetplustracks_lxy, "lxy/F");
   m_mujetplustracks->Branch("lxyz", &m_mujetplustracks_lxyz, "lxyz/F");
   m_mujetplustracks->Branch("iso", &m_mujetplustracks_iso, "iso/F");
   m_mujetplustracks->Branch("mupt1", &m_mujetplustracks_mupt1, "mupt1/F");
   m_mujetplustracks->Branch("mupt2", &m_mujetplustracks_mupt2, "mupt2/F");
   m_mujetplustracks->Branch("mupt3", &m_mujetplustracks_mupt3, "mupt3/F");
   m_mujetplustracks->Branch("fakept1", &m_mujetplustracks_fakept1, "fakept1/F");
   m_mujetplustracks->Branch("fakept2", &m_mujetplustracks_fakept2, "fakept2/F");

   m_quadmuon = tFile->make<TTree>("quadmuon", "quadmuon");
   m_quadmuon->Branch("qscale", &m_qscale, "qscale/F");
   m_quadmuon->Branch("genLevel", &m_genLevel, "genLevel/I");
   m_quadmuon->Branch("run", &m_run, "run/I");
   m_quadmuon->Branch("event", &m_event, "event/I");
   m_quadmuon->Branch("trigger", &m_trigger, "trigger/I");
   m_quadmuon->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
   m_quadmuon->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
   m_quadmuon->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
   m_quadmuon->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
   m_quadmuon->Branch("muon3pt", &m_muon3pt, "muon3pt/F");
   m_quadmuon->Branch("muon3eta", &m_muon3eta, "muon3eta/F");
   m_quadmuon->Branch("muon4pt", &m_muon4pt, "muon4pt/F");
   m_quadmuon->Branch("muon4eta", &m_muon4eta, "muon4eta/F");
   m_quadmuon->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
   m_quadmuon->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
   m_quadmuon->Branch("containstrig", &m_quadmuon_containstrig, "containstrig/I");
   m_quadmuon->Branch("muons", &m_quadmuon_muons, "muons/I");
   m_quadmuon->Branch("orphans", &m_quadmuon_orphans, "orphans/I");
   m_quadmuon->Branch("charge", &m_quadmuon_charge, "charge/I");
   m_quadmuon->Branch("mass", &m_quadmuon_mass, "mass/F");
   m_quadmuon->Branch("massa", &m_quadmuon_massa, "massa/F");
   m_quadmuon->Branch("massb", &m_quadmuon_massb, "massb/F");
   m_quadmuon->Branch("massc", &m_quadmuon_massc, "massc/F");
   m_quadmuon->Branch("massd", &m_quadmuon_massd, "massd/F");
   m_quadmuon->Branch("masse", &m_quadmuon_masse, "masse/F");
   m_quadmuon->Branch("massf", &m_quadmuon_massf, "massf/F");
   m_quadmuon->Branch("pt", &m_quadmuon_pt, "pt/F");
   m_quadmuon->Branch("eta", &m_quadmuon_eta, "eta/F");
   m_quadmuon->Branch("phi", &m_quadmuon_phi, "phi/F");
   m_quadmuon->Branch("dr", &m_quadmuon_dr, "dr/F");
   m_quadmuon->Branch("vchi2", &m_quadmuon_vchi2, "vchi2/F");
   m_quadmuon->Branch("vndof", &m_quadmuon_vndof, "vndof/I");
   m_quadmuon->Branch("vprob", &m_quadmuon_vprob, "vprob/F");
   m_quadmuon->Branch("lxy", &m_quadmuon_lxy, "lxy/F");
   m_quadmuon->Branch("lxyz", &m_quadmuon_lxyz, "lxyz/F");
   m_quadmuon->Branch("iso", &m_quadmuon_iso, "iso/F");

   m_dimuorphan = tFile->make<TTree>("dimuorphan", "dimuorphan");
   m_dimuorphan->Branch("qscale", &m_qscale, "qscale/F");
   m_dimuorphan->Branch("genLevel", &m_genLevel, "genLevel/I");
   m_dimuorphan->Branch("run", &m_run, "run/I");
   m_dimuorphan->Branch("event", &m_event, "event/I");
   m_dimuorphan->Branch("trigger", &m_trigger, "trigger/I");
   m_dimuorphan->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
   m_dimuorphan->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
   m_dimuorphan->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
   m_dimuorphan->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
   m_dimuorphan->Branch("muon3pt", &m_muon3pt, "muon3pt/F");
   m_dimuorphan->Branch("muon3eta", &m_muon3eta, "muon3eta/F");
   m_dimuorphan->Branch("muon4pt", &m_muon4pt, "muon4pt/F");
   m_dimuorphan->Branch("muon4eta", &m_muon4eta, "muon4eta/F");
   m_dimuorphan->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
   m_dimuorphan->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
   m_dimuorphan->Branch("deltaphi", &m_dimuorphan_deltaphi, "deltaphi/F");
   m_dimuorphan->Branch("orphanpt", &m_dimuorphan_orphanpt, "orphanpt/F");
   m_dimuorphan->Branch("orphaneta", &m_dimuorphan_orphaneta, "orphaneta/F");
   m_dimuorphan->Branch("orphanphi", &m_dimuorphan_orphanphi, "orphanphi/F");
   m_dimuorphan->Branch("orphanisglobal", &m_dimuorphan_orphanisglobal, "orphanisglobal/I");
   m_dimuorphan->Branch("orphanmatches", &m_dimuorphan_orphanmatches, "orphan/I");
   m_dimuorphan->Branch("orphanhits", &m_dimuorphan_orphanhits, "orphanhits/I");
   m_dimuorphan->Branch("orphanchi2", &m_dimuorphan_orphanchi2, "orphanchi2/F");
   m_dimuorphan->Branch("orphanst1x", &m_dimuorphan_orphanst1x, "orphanst1x/F");
   m_dimuorphan->Branch("orphanst1xsig", &m_dimuorphan_orphanst1xsig, "orphanst1xsig/F");
   m_dimuorphan->Branch("orphanst1dxdz", &m_dimuorphan_orphanst1dxdz, "orphanst1dxdz/F");
   m_dimuorphan->Branch("orphanst1dxdzsig", &m_dimuorphan_orphanst1dxdzsig, "orphanst1dxdzsig/F");
   m_dimuorphan->Branch("containstrig", &m_dimuorphan_containstrig, "containstrig/I");
   m_dimuorphan->Branch("mass", &m_dimuorphan_mass, "mass/F");
   m_dimuorphan->Branch("pt", &m_dimuorphan_pt, "pt/F");
   m_dimuorphan->Branch("eta", &m_dimuorphan_eta, "eta/F");
   m_dimuorphan->Branch("phi", &m_dimuorphan_phi, "phi/F");
   m_dimuorphan->Branch("dr", &m_dimuorphan_dr, "dr/F");
   m_dimuorphan->Branch("vprob", &m_dimuorphan_vprob, "vprob/F");
   m_dimuorphan->Branch("lxy", &m_dimuorphan_lxy, "lxy/F");
   m_dimuorphan->Branch("lxyz", &m_dimuorphan_lxyz, "lxyz/F");
   m_dimuorphan->Branch("iso", &m_dimuorphan_iso, "iso/F");

   m_dimudimu = tFile->make<TTree>("dimudimu", "dimudimu");
   m_dimudimu->Branch("qscale", &m_qscale, "qscale/F");
   m_dimudimu->Branch("genLevel", &m_genLevel, "genLevel/I");
   m_dimudimu->Branch("run", &m_run, "run/I");
   m_dimudimu->Branch("event", &m_event, "event/I");
   m_dimudimu->Branch("trigger", &m_trigger, "trigger/I");
   m_dimudimu->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
   m_dimudimu->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
   m_dimudimu->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
   m_dimudimu->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
   m_dimudimu->Branch("muon3pt", &m_muon3pt, "muon3pt/F");
   m_dimudimu->Branch("muon3eta", &m_muon3eta, "muon3eta/F");
   m_dimudimu->Branch("muon4pt", &m_muon4pt, "muon4pt/F");
   m_dimudimu->Branch("muon4eta", &m_muon4eta, "muon4eta/F");
   m_dimudimu->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
   m_dimudimu->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
   m_dimudimu->Branch("orphans", &m_dimudimu_orphans, "orphans/I");
   m_dimudimu->Branch("wholemass", &m_dimudimu_wholemass, "wholemass/F");
   m_dimudimu->Branch("wholept", &m_dimudimu_wholept, "wholept/F");
   m_dimudimu->Branch("deltaphi", &m_dimudimu_deltaphi, "deltaphi/F");
   m_dimudimu->Branch("containstrig", &m_dimudimu_containstrig, "containstrig/I");
   m_dimudimu->Branch("massC", &m_dimudimu_massC, "massC/F");
   m_dimudimu->Branch("ptC", &m_dimudimu_ptC, "ptC/F");
   m_dimudimu->Branch("etaC", &m_dimudimu_etaC, "etaC/F");
   m_dimudimu->Branch("phiC", &m_dimudimu_phiC, "phiC/F");
   m_dimudimu->Branch("drC", &m_dimudimu_drC, "drC/F");
   m_dimudimu->Branch("vprobC", &m_dimudimu_vprobC, "vprobC/F");
   m_dimudimu->Branch("lxyC", &m_dimudimu_lxyC, "lxyC/F");
   m_dimudimu->Branch("lxyzC", &m_dimudimu_lxyzC, "lxyzC/F");
   m_dimudimu->Branch("isoC", &m_dimudimu_isoC, "isoC/F");
   m_dimudimu->Branch("massF", &m_dimudimu_massF, "massF/F");
   m_dimudimu->Branch("ptF", &m_dimudimu_ptF, "ptF/F");
   m_dimudimu->Branch("etaF", &m_dimudimu_etaF, "etaF/F");
   m_dimudimu->Branch("phiF", &m_dimudimu_phiF, "phiF/F");
   m_dimudimu->Branch("drF", &m_dimudimu_drF, "drF/F");
   m_dimudimu->Branch("vprobF", &m_dimudimu_vprobF, "vprobF/F");
   m_dimudimu->Branch("lxyF", &m_dimudimu_lxyF, "lxyF/F");
   m_dimudimu->Branch("lxyzF", &m_dimudimu_lxyzF, "lxyzF/F");
   m_dimudimu->Branch("isoF", &m_dimudimu_isoF, "isoF/F");

   m_twomujets = tFile->make<TTree>("twomujets", "twomujets");
   m_twomujets->Branch("qscale", &m_qscale, "qscale/F");
   m_twomujets->Branch("genLevel", &m_genLevel, "genLevel/I");
   m_twomujets->Branch("run", &m_run, "run/I");
   m_twomujets->Branch("event", &m_event, "event/I");
   m_twomujets->Branch("trigger", &m_trigger, "trigger/I");
   m_twomujets->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
   m_twomujets->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
   m_twomujets->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
   m_twomujets->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
   m_twomujets->Branch("muon3pt", &m_muon3pt, "muon3pt/F");
   m_twomujets->Branch("muon3eta", &m_muon3eta, "muon3eta/F");
   m_twomujets->Branch("muon4pt", &m_muon4pt, "muon4pt/F");
   m_twomujets->Branch("muon4eta", &m_muon4eta, "muon4eta/F");
   m_twomujets->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
   m_twomujets->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
   m_twomujets->Branch("orphans", &m_dimudimu_orphans, "orphans/I");
   m_twomujets->Branch("wholemass", &m_dimudimu_wholemass, "wholemass/F");
   m_twomujets->Branch("wholept", &m_dimudimu_wholept, "wholept/F");
   m_twomujets->Branch("deltaphi", &m_dimudimu_deltaphi, "deltaphi/F");
   m_twomujets->Branch("containstrig", &m_dimudimu_containstrig, "containstrig/I");
   m_twomujets->Branch("massC", &m_dimudimu_massC, "massC/F");
   m_twomujets->Branch("ptC", &m_dimudimu_ptC, "ptC/F");
   m_twomujets->Branch("etaC", &m_dimudimu_etaC, "etaC/F");
   m_twomujets->Branch("phiC", &m_dimudimu_phiC, "phiC/F");
   m_twomujets->Branch("drC", &m_dimudimu_drC, "drC/F");
   m_twomujets->Branch("vprobC", &m_dimudimu_vprobC, "vprobC/F");
   m_twomujets->Branch("lxyC", &m_dimudimu_lxyC, "lxyC/F");
   m_twomujets->Branch("lxyzC", &m_dimudimu_lxyzC, "lxyzC/F");
   m_twomujets->Branch("isoC", &m_dimudimu_isoC, "isoC/F");
   m_twomujets->Branch("massF", &m_dimudimu_massF, "massF/F");
   m_twomujets->Branch("ptF", &m_dimudimu_ptF, "ptF/F");
   m_twomujets->Branch("etaF", &m_dimudimu_etaF, "etaF/F");
   m_twomujets->Branch("phiF", &m_dimudimu_phiF, "phiF/F");
   m_twomujets->Branch("drF", &m_dimudimu_drF, "drF/F");
   m_twomujets->Branch("vprobF", &m_dimudimu_vprobF, "vprobF/F");
   m_twomujets->Branch("lxyF", &m_dimudimu_lxyF, "lxyF/F");
   m_twomujets->Branch("lxyzF", &m_dimudimu_lxyzF, "lxyzF/F");
   m_twomujets->Branch("isoF", &m_dimudimu_isoF, "isoF/F");
   m_twomujets->Branch("muonsC", &m_twomujets_muonsC, "muonsC/I");
   m_twomujets->Branch("muonsF", &m_twomujets_muonsF, "muonsF/I");
   m_twomujets->Branch("massCa", &m_twomujets_massCa, "massCa/F");
   m_twomujets->Branch("massCb", &m_twomujets_massCb, "massCb/F");
   m_twomujets->Branch("massCc", &m_twomujets_massCc, "massCc/F");
   m_twomujets->Branch("massFa", &m_twomujets_massFa, "massFa/F");
   m_twomujets->Branch("massFb", &m_twomujets_massFb, "massFb/F");
   m_twomujets->Branch("massFc", &m_twomujets_massFc, "massFc/F");

   m_moremujets = tFile->make<TTree>("moremujets", "moremujets");
   m_moremujets->Branch("qscale", &m_qscale, "qscale/F");
   m_moremujets->Branch("genLevel", &m_genLevel, "genLevel/I");
   m_moremujets->Branch("run", &m_run, "run/I");
   m_moremujets->Branch("event", &m_event, "event/I");
   m_moremujets->Branch("trigger", &m_trigger, "trigger/I");
   m_moremujets->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
   m_moremujets->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
   m_moremujets->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
   m_moremujets->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
   m_moremujets->Branch("muon3pt", &m_muon3pt, "muon3pt/F");
   m_moremujets->Branch("muon3eta", &m_muon3eta, "muon3eta/F");
   m_moremujets->Branch("muon4pt", &m_muon4pt, "muon4pt/F");
   m_moremujets->Branch("muon4eta", &m_muon4eta, "muon4eta/F");
   m_moremujets->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
   m_moremujets->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
   m_moremujets->Branch("mujets", &m_moremujets_mujets, "mujets/I");
   m_moremujets->Branch("orphans", &m_moremujets_orphans, "orphans/I");
   m_moremujets->Branch("wholemass", &m_moremujets_wholemass, "wholemass/F");
   m_moremujets->Branch("wholept", &m_moremujets_wholept, "wholept/F");
   m_moremujets->Branch("massa", &m_moremujets_massa, "massa/F");
   m_moremujets->Branch("massb", &m_moremujets_massb, "massb/F");
   m_moremujets->Branch("massc", &m_moremujets_massc, "massc/F");
   m_moremujets->Branch("massd", &m_moremujets_massd, "massd/F");
   m_moremujets->Branch("masse", &m_moremujets_masse, "masse/F");
   m_moremujets->Branch("massf", &m_moremujets_massf, "massf/F");
   m_moremujets->Branch("muons1", &m_moremujets_muons1, "muons1/I");
   m_moremujets->Branch("containstrig1", &m_moremujets_containstrig1, "containstrig1/I");
   m_moremujets->Branch("mass1", &m_moremujets_mass1, "mass1/F");
   m_moremujets->Branch("pt1", &m_moremujets_pt1, "pt1/F");
   m_moremujets->Branch("iso1", &m_moremujets_iso1, "iso1/F");
   m_moremujets->Branch("muons2", &m_moremujets_muons2, "muons2/I");
   m_moremujets->Branch("containstrig2", &m_moremujets_containstrig2, "containstrig2/I");
   m_moremujets->Branch("mass2", &m_moremujets_mass2, "mass2/F");
   m_moremujets->Branch("pt2", &m_moremujets_pt2, "pt2/F");
   m_moremujets->Branch("iso2", &m_moremujets_iso2, "iso2/F");
   m_moremujets->Branch("muons3", &m_moremujets_muons3, "muons3/I");
   m_moremujets->Branch("containstrig3", &m_moremujets_containstrig3, "containstrig3/I");
   m_moremujets->Branch("mass3", &m_moremujets_mass3, "mass3/F");
   m_moremujets->Branch("pt3", &m_moremujets_pt3, "pt3/F");
   m_moremujets->Branch("iso3", &m_moremujets_iso3, "iso3/F");
   m_moremujets->Branch("muons4", &m_moremujets_muons4, "muons4/I");
   m_moremujets->Branch("containstrig4", &m_moremujets_containstrig4, "containstrig4/I");
   m_moremujets->Branch("mass4", &m_moremujets_mass4, "mass4/F");
   m_moremujets->Branch("pt4", &m_moremujets_pt4, "pt4/F");
   m_moremujets->Branch("iso4", &m_moremujets_iso4, "iso4/F");
   m_moremujets->Branch("muons5", &m_moremujets_muons5, "muons5/I");
   m_moremujets->Branch("containstrig5", &m_moremujets_containstrig5, "containstrig5/I");
   m_moremujets->Branch("mass5", &m_moremujets_mass5, "mass5/F");
   m_moremujets->Branch("pt5", &m_moremujets_pt5, "pt5/F");
   m_moremujets->Branch("iso5", &m_moremujets_iso5, "iso5/F");
   m_moremujets->Branch("muons6", &m_moremujets_muons6, "muons6/I");
   m_moremujets->Branch("containstrig6", &m_moremujets_containstrig6, "containstrig6/I");
   m_moremujets->Branch("mass6", &m_moremujets_mass6, "mass6/F");
   m_moremujets->Branch("pt6", &m_moremujets_pt6, "pt6/F");
   m_moremujets->Branch("iso6", &m_moremujets_iso6, "iso6/F");

   m_displaced = tFile->make<TTree>("displaced", "displaced");
   m_displaced->Branch("qscale", &m_qscale, "qscale/F");
   m_displaced->Branch("genLevel", &m_genLevel, "genLevel/I");
   m_displaced->Branch("run", &m_run, "run/I");
   m_displaced->Branch("event", &m_event, "event/I");
   m_displaced->Branch("trigger", &m_trigger, "trigger/I");
   m_displaced->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
   m_displaced->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
   m_displaced->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
   m_displaced->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
   m_displaced->Branch("muon3pt", &m_muon3pt, "muon3pt/F");
   m_displaced->Branch("muon3eta", &m_muon3eta, "muon3eta/F");
   m_displaced->Branch("muon4pt", &m_muon4pt, "muon4pt/F");
   m_displaced->Branch("muon4eta", &m_muon4eta, "muon4eta/F");
   m_displaced->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
   m_displaced->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
   m_displaced->Branch("mujets", &m_displaced_mujets, "mujets/I");
   m_displaced->Branch("orphans", &m_displaced_orphans, "orphans/I");
   m_displaced->Branch("muons", &m_displaced_muons, "muons/I");
   m_displaced->Branch("mass", &m_displaced_mass, "mass/F");
   m_displaced->Branch("pt", &m_displaced_pt, "pt/F");
   m_displaced->Branch("eta", &m_displaced_eta, "eta/F");
   m_displaced->Branch("phi", &m_displaced_phi, "phi/F");
   m_displaced->Branch("vx", &m_displaced_vx, "vx/F");
   m_displaced->Branch("vy", &m_displaced_vy, "vy/F");
   m_displaced->Branch("vz", &m_displaced_vz, "vz/F");
   m_displaced->Branch("zclosest", &m_displaced_zclosest, "zclosest/F");
   m_displaced->Branch("zcloseness", &m_displaced_zcloseness, "zcloseness/F");
   m_displaced->Branch("lxy", &m_displaced_lxy, "lxy/F");
   m_displaced->Branch("lxyz", &m_displaced_lxyz, "lxyz/F");
   m_displaced->Branch("iso", &m_displaced_iso, "iso/F");
   m_displaced->Branch("vprob", &m_displaced_vprob, "vprob/F");
   m_displaced->Branch("mupt1", &m_displaced_mupt1, "mupt1/F");
   m_displaced->Branch("mupt2", &m_displaced_mupt2, "mupt2/F");
   m_displaced->Branch("plusmatches", &m_displaced_plusmatches, "plusmatches/I");
   m_displaced->Branch("plusst1x", &m_displaced_plusst1x, "plusst1x/F");
   m_displaced->Branch("plusst1xsig", &m_displaced_plusst1xsig, "plusst1xsig/F");
   m_displaced->Branch("plusst1dxdz", &m_displaced_plusst1dxdz, "plusst1dxdz/F");
   m_displaced->Branch("plusst1dxdzsig", &m_displaced_plusst1dxdzsig, "plusst1dxdzsig/F");
   m_displaced->Branch("plushits", &m_displaced_plushits, "plushits/I");
   m_displaced->Branch("plusnormchi2", &m_displaced_plusnormchi2, "plusnormchi2/F");
   m_displaced->Branch("minusmatches", &m_displaced_minusmatches, "minusmatches/I");
   m_displaced->Branch("minusst1x", &m_displaced_minusst1x, "minusst1x/F");
   m_displaced->Branch("minusst1xsig", &m_displaced_minusst1xsig, "minusst1xsig/F");
   m_displaced->Branch("minusst1dxdz", &m_displaced_minusst1dxdz, "minusst1dxdz/F");
   m_displaced->Branch("minusst1dxdzsig", &m_displaced_minusst1dxdzsig, "minusst1dxdzsig/F");
   m_displaced->Branch("minushits", &m_displaced_minushits, "minushits/I");
   m_displaced->Branch("minusnormchi2", &m_displaced_minusnormchi2, "minusnormchi2/F");

   m_displacedFake = tFile->make<TTree>("displacedFake", "displacedFake");
   m_displacedFake->Branch("qscale", &m_qscale, "qscale/F");
   m_displacedFake->Branch("genLevel", &m_genLevel, "genLevel/I");
   m_displacedFake->Branch("run", &m_run, "run/I");
   m_displacedFake->Branch("event", &m_event, "event/I");
   m_displacedFake->Branch("trigger", &m_trigger, "trigger/I");
   m_displacedFake->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
   m_displacedFake->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
   m_displacedFake->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
   m_displacedFake->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
   m_displacedFake->Branch("muon3pt", &m_muon3pt, "muon3pt/F");
   m_displacedFake->Branch("muon3eta", &m_muon3eta, "muon3eta/F");
   m_displacedFake->Branch("muon4pt", &m_muon4pt, "muon4pt/F");
   m_displacedFake->Branch("muon4eta", &m_muon4eta, "muon4eta/F");
   m_displacedFake->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
   m_displacedFake->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
   m_displacedFake->Branch("mujets", &m_displaced_mujets, "mujets/I");
   m_displacedFake->Branch("orphans", &m_displaced_orphans, "orphans/I");
   m_displacedFake->Branch("muons", &m_displaced_muons, "muons/I");
   m_displacedFake->Branch("mass", &m_displaced_mass, "mass/F");
   m_displacedFake->Branch("pt", &m_displaced_pt, "pt/F");
   m_displacedFake->Branch("eta", &m_displaced_eta, "eta/F");
   m_displacedFake->Branch("phi", &m_displaced_phi, "phi/F");
   m_displacedFake->Branch("vx", &m_displaced_vx, "vx/F");
   m_displacedFake->Branch("vy", &m_displaced_vy, "vy/F");
   m_displacedFake->Branch("vz", &m_displaced_vz, "vz/F");
   m_displacedFake->Branch("zclosest", &m_displaced_zclosest, "zclosest/F");
   m_displacedFake->Branch("zcloseness", &m_displaced_zcloseness, "zcloseness/F");
   m_displacedFake->Branch("lxy", &m_displaced_lxy, "lxy/F");
   m_displacedFake->Branch("lxyz", &m_displaced_lxyz, "lxyz/F");
   m_displacedFake->Branch("iso", &m_displaced_iso, "iso/F");
   m_displacedFake->Branch("vprob", &m_displaced_vprob, "vprob/F");
   m_displacedFake->Branch("mupt1", &m_displaced_mupt1, "mupt1/F");
   m_displacedFake->Branch("mupt2", &m_displaced_mupt2, "mupt2/F");
   m_displacedFake->Branch("plusmatches", &m_displaced_plusmatches, "plusmatches/I");
   m_displacedFake->Branch("plusst1x", &m_displaced_plusst1x, "plusst1x/F");
   m_displacedFake->Branch("plusst1xsig", &m_displaced_plusst1xsig, "plusst1xsig/F");
   m_displacedFake->Branch("plusst1dxdz", &m_displaced_plusst1dxdz, "plusst1dxdz/F");
   m_displacedFake->Branch("plusst1dxdzsig", &m_displaced_plusst1dxdzsig, "plusst1dxdzsig/F");
   m_displacedFake->Branch("plushits", &m_displaced_plushits, "plushits/I");
   m_displacedFake->Branch("plusnormchi2", &m_displaced_plusnormchi2, "plusnormchi2/F");
   m_displacedFake->Branch("minusmatches", &m_displaced_minusmatches, "minusmatches/I");
   m_displacedFake->Branch("minusst1x", &m_displaced_minusst1x, "minusst1x/F");
   m_displacedFake->Branch("minusst1xsig", &m_displaced_minusst1xsig, "minusst1xsig/F");
   m_displacedFake->Branch("minusst1dxdz", &m_displaced_minusst1dxdz, "minusst1dxdz/F");
   m_displacedFake->Branch("minusst1dxdzsig", &m_displaced_minusst1dxdzsig, "minusst1dxdzsig/F");
   m_displacedFake->Branch("minushits", &m_displaced_minushits, "minushits/I");
   m_displacedFake->Branch("minusnormchi2", &m_displaced_minusnormchi2, "minusnormchi2/F");
   m_displacedFake->Branch("fakes", &m_displacedFake_fakes, "fakes/I");
}


FitNtuple::~FitNtuple()
{
 
   // do anything here that needs to be done at desctruction time
   // (e.g. close files, deallocate resources etc.)

}


//
// member functions
//

// ------------ method called to for each event  ------------
void
FitNtuple::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{
   // the QCD scale only exists for certain kinds of Monte Carlo, but if it's available, get it
   m_qscale = 0.;
   if (m_getQscale) {
      edm::Handle<GenEventInfoProduct> genEventInfoProduct;
      iEvent.getByLabel("generator", genEventInfoProduct);
      m_qscale = genEventInfoProduct->qScale();
   }

   // the pair-gun Monte Carlos have a (never used) feature called
   // alteration; the real events are the ones in which particleNumber == 0
   bool alternating = true;
   if (m_getAlternating) {
      edm::Handle<unsigned int> particleNumber;
      iEvent.getByLabel("generator", "particleNumber", particleNumber);
      alternating = (*particleNumber == 0);
   }
   if (!alternating) return;

   // get the run and event number
   m_run = iEvent.id().run();
   m_event = iEvent.id().event();

   // mu-jets (muons grouped by mass and vertex compatibility)
   edm::Handle<pat::MultiMuonCollection> muJets;
   iEvent.getByLabel("MuJetProducer", muJets);

   // orphans (muons not found in any group)
   edm::Handle<pat::MuonCollection> orphans;
   iEvent.getByLabel("MuJetProducer", "Orphans", orphans);

   // mu-jets plus at most two tracks as fake muons
   edm::Handle<pat::MultiMuonCollection> muJetPlusTracks;
   iEvent.getByLabel("MuJetPlusTracks", muJetPlusTracks);

   edm::Handle<pat::MultiMuonCollection> displacedMuJet;
   iEvent.getByLabel("DisplacedMuJet", displacedMuJet);

   edm::Handle<pat::MultiMuonCollection> displacedMuJetFake;
   iEvent.getByLabel("DisplacedMuJetFake", displacedMuJetFake);

   // find the top four muons in the event (-1000. if not found)
   edm::Handle<pat::MuonCollection> muons;
   iEvent.getByLabel("cleanPatMuonsTriggerMatch", muons);
   m_muon1pt = -1000.; m_muon1eta = -1000.;
   m_muon2pt = -1000.; m_muon2eta = -1000.;
   m_muon3pt = -1000.; m_muon3eta = -1000.;
   m_muon4pt = -1000.; m_muon4eta = -1000.;
   m_muontrigpt = -1000.; m_muontrigeta = -1000.;
   std::vector<pat::MuonCollection::const_iterator> muontrig;

   for (pat::MuonCollection::const_iterator muon = muons->begin();  muon != muons->end();  ++muon) {
      if (fabs(muon->eta()) < 2.4  &&  muon->isTrackerMuon() && muon->numberOfMatches(reco::Muon::SegmentAndTrackArbitration) >= 2 && muon->innerTrack()->numberOfValidHits() >= 8 && muon->innerTrack()->normalizedChi2() < 4.) {
         if (muon->pt() > m_muon1pt) {
            m_muon4pt = m_muon3pt; m_muon4eta = m_muon3eta;
            m_muon3pt = m_muon2pt; m_muon3eta = m_muon2eta;
            m_muon2pt = m_muon1pt; m_muon2eta = m_muon1eta;
            m_muon1pt = muon->pt();  m_muon1eta = muon->eta();
         }
         else if (muon->pt() > m_muon2pt) {
            m_muon4pt = m_muon3pt; m_muon4eta = m_muon3eta;
            m_muon3pt = m_muon2pt; m_muon3eta = m_muon2eta;
            m_muon2pt = muon->pt(); m_muon2eta = muon->eta();
         }
         else if (muon->pt() > m_muon3pt) {
            m_muon4pt = m_muon3pt; m_muon4eta = m_muon3eta;
            m_muon3pt = muon->pt(); m_muon3eta = muon->eta();
         }
         else if (muon->pt() > m_muon4pt) {
            m_muon4pt = muon->pt(); m_muon4eta = muon->eta();
         }

         // special muon within a more limited |eta| range, to guarantee the trigger
         if (fabs(muon->eta()) < 0.9) {
            if (muon->pt() > m_muontrigpt) {
               m_muontrigpt = muon->pt(); m_muontrigeta = muon->eta();
            }
         }

         // all muons with pT > 15 GeV/c and |eta| < 0.9, matched to a trigger object
         if (muon->pt() > 15.  &&  fabs(muon->eta()) < 0.9) {
            const pat::TriggerObjectStandAlone *mu3 = muon->triggerObjectMatchByPath("HLT_Mu3");
            const pat::TriggerObjectStandAlone *mu5 = muon->triggerObjectMatchByPath("HLT_Mu5");
            const pat::TriggerObjectStandAlone *mu9 = muon->triggerObjectMatchByPath("HLT_Mu9");
            const pat::TriggerObjectStandAlone *mu11 = muon->triggerObjectMatchByPath("HLT_Mu11");
            const pat::TriggerObjectStandAlone *mu15 = muon->triggerObjectMatchByPath("HLT_Mu15_v1");

            if ((mu3 != NULL  &&  mu3->collection() == std::string("hltL3MuonCandidates::HLT")  &&  mu3->pt() > 15.)  ||
                (mu5 != NULL  &&  mu5->collection() == std::string("hltL3MuonCandidates::HLT")  &&  mu5->pt() > 15.)  ||
                (mu9 != NULL  &&  mu9->collection() == std::string("hltL3MuonCandidates::HLT")  &&  mu9->pt() > 15.)  ||
                (mu11 != NULL  &&  mu11->collection() == std::string("hltL3MuonCandidates::HLT")  &&  mu11->pt() > 15.)  ||
                (mu15 != NULL  &&  mu15->collection() == std::string("hltL3MuonCandidates::HLT")  &&  mu15->pt() > 15.)) {
               muontrig.push_back(muon);
            }
         }
      }
   }

//    // all tracker-tracks
//    edm::Handle<reco::TrackCollection> tracks;
//    iEvent.getByLabel("generalTracks", tracks);

//    // generator-level 4-vectors
//    edm::Handle<reco::GenParticleCollection> genParticles;
//    iEvent.getByLabel("genParticles", genParticles);

   // find the closest primary vertex (in Z) to the first muJet with a valid vertex
   edm::Handle<reco::VertexCollection> primaryVertices;
   iEvent.getByLabel("offlinePrimaryVertices", primaryVertices);
   reco::VertexCollection::const_iterator closestPrimaryVertex = primaryVertices->end();
   if (muJets->size() > 0) {
      pat::MultiMuonCollection::const_iterator muJet0 = muJets->end();
      for (pat::MultiMuonCollection::const_iterator muJet = muJets->begin();  muJet != muJets->end();  ++muJet) {
         if (muJet->vertexValid()) {
            muJet0 = muJet;
            break;
         }
      }

      if (muJet0 != muJets->end()) {
         for (reco::VertexCollection::const_iterator vertex = primaryVertices->begin();  vertex != primaryVertices->end();  ++vertex) {
            if (vertex->isValid()  &&  !vertex->isFake()  &&  vertex->tracksSize() > 3  &&  fabs(vertex->z()) < 24.) {
               if (closestPrimaryVertex == primaryVertices->end()  ||  fabs(vertex->z() - muJet0->vertexPoint().z()) < fabs(closestPrimaryVertex->z() - muJet0->vertexPoint().z())) {
                  closestPrimaryVertex = vertex;
               }
            } // end vertex quality cuts
         } // end loop over primary vertices
      } // end if muJet0 exists
   } // end if muJets->size > 0

   // find out which trigger bits were fired
   edm::Handle<pat::TriggerEvent> triggerEvent;
   iEvent.getByLabel("patTriggerEvent", triggerEvent);
   m_trigger = 0;
   if (triggerEvent->path("HLT_Mu5")  &&  triggerEvent->path("HLT_Mu5")->wasAccept()) m_trigger += 1;
   if (triggerEvent->path("HLT_Mu9")  &&  triggerEvent->path("HLT_Mu9")->wasAccept()) m_trigger += 2;
   if (triggerEvent->path("HLT_Mu11")  &&  triggerEvent->path("HLT_Mu11")->wasAccept()) m_trigger += 4;
   if (triggerEvent->path("HLT_Mu15_v1")  &&  triggerEvent->path("HLT_Mu15_v1")->wasAccept()) m_trigger += 8;
   if (triggerEvent->path("HLT_Mu20_NoVertex")  &&  triggerEvent->path("HLT_Mu20_NoVertex")->wasAccept()) m_trigger += 16;

   // get some gen-level information
   m_genLevel = -1;
   if (m_fillGenLevel) {
      bool any_missing_mcmatch = false;
      bool any_decays_in_flight = false;
      unsigned int bjets = 0;
      unsigned int resonances = 0;
      for (pat::MultiMuonCollection::const_iterator muJet = muJets->begin();  muJet != muJets->end();  ++muJet) {
         for (unsigned int i = 0;  i < muJet->numberOfDaughters();  i++) {
            const pat::Muon *muon = dynamic_cast<const pat::Muon*>(muJet->daughter(i));

            if (muon->genParticlesSize() == 0) any_missing_mcmatch = true;

            else if (muon->genParticle()->mother()) {
               int mother = muon->genParticle()->mother()->pdgId();
               if (abs(mother) == 13  &&  muon->genParticle()->mother()->mother()) muon->genParticle()->mother()->mother()->pdgId();  // don't let muons be their own mothers

               if (abs(mother) == 211 /* charged pion */  ||
                   abs(mother) == 321 /* charged kaon */  ||
                   (3000 <= abs(mother)  &&  abs(mother) < 4000) /* strange baryon */ ) any_decays_in_flight = true;
            }
         }

         if (muJet->genParticlesSize() == 1) {
            for (const reco::Candidate *particle = muJet->genParticle();  particle != NULL;  particle = particle->mother()) {
               int pdgId = particle->pdgId();
               if (abs(pdgId) == 511  ||
                   abs(pdgId) == 521  ||
                   abs(pdgId) == 513  ||
                   abs(pdgId) == 523  ||
                   abs(pdgId) == 515  ||
                   abs(pdgId) == 525  ||
                   abs(pdgId) == 531  ||
                   abs(pdgId) == 533  ||
                   abs(pdgId) == 535  ||
                   abs(pdgId) == 541  ||
                   abs(pdgId) == 543  ||
                   abs(pdgId) == 545  ||
                   (5000 <= abs(pdgId)  &&  abs(pdgId) < 6000)) {
                  bjets++;
                  break;
               }
            }

            if (muJet->genParticle()->pdgId() == 113 /* rho */  ||
                muJet->genParticle()->pdgId() == 221 /* eta */  ||
                muJet->genParticle()->pdgId() == 223 /* omega(782) */  ||
                muJet->genParticle()->pdgId() == 333 /* phi(1020) */  ||
                muJet->genParticle()->pdgId() == 443 /* J/psi */  ||
                muJet->genParticle()->pdgId() == 100443 /* psi(2S) */ ) {
               resonances++;
            }
         }
      }

      m_genLevel = 0;
      if (bjets > 0  &&  bjets == muJets->size()) m_genLevel += 1;
      if (resonances > 0  &&  resonances == muJets->size()) m_genLevel += 2;
      if (any_decays_in_flight) m_genLevel += 4;
      if (any_missing_mcmatch) m_genLevel += 8;
   }

   ////////////////////////////////////////////////////////// alldimuons (special case for scaling to the Z)
   if (muons->size() == 2) {
      double total_energy = (*muons)[0].energy() + (*muons)[1].energy();
      double total_px = (*muons)[0].px() + (*muons)[1].px();
      double total_py = (*muons)[0].py() + (*muons)[1].py();
      double total_pz = (*muons)[0].pz() + (*muons)[1].pz();
      m_alldimuons_mass = sqrt(total_energy*total_energy - total_px*total_px - total_py*total_py - total_pz*total_pz);
      if (muJets->size() == 1  &&  (*muJets)[0].numberOfDaughters() == 2) {
         m_alldimuons_iso = (*muJets)[0].centralTrackIsolation();
      }
      else {
         m_alldimuons_iso = -1.;
      }
      m_alldimuons->Fill();
   }

   ////////////////////////////////////////////////////////// lowdimuon and highdimuon:
   if (muJets->size() == 1  &&  (*muJets)[0].numberOfDaughters() == 2  &&  orphans->size() == 0) {
      pat::MultiMuonCollection::const_iterator muJet = muJets->begin();

      m_lowdimuon_containstrig = 0;
      for (std::vector<pat::MuonCollection::const_iterator>::const_iterator iter = muontrig.begin();  iter != muontrig.end();  ++iter) {
         if (muJet->muon(0)->innerTrack().isAvailable()  &&  (*iter)->innerTrack().isAvailable()  &&
             muJet->sameTrack(&*(muJet->muon(0)->innerTrack()), &*((*iter)->innerTrack()))) {
            m_lowdimuon_containstrig = 1;
         }
         if (muJet->muon(1)->innerTrack().isAvailable()  &&  (*iter)->innerTrack().isAvailable()  &&
             muJet->sameTrack(&*(muJet->muon(1)->innerTrack()), &*((*iter)->innerTrack()))) {
            m_lowdimuon_containstrig = 1;
         }
      }

      // generator-level mass using matched genParticles (for resolution of fit peak)
      m_lowdimuon_genmass = -1000.;
      if (muJet->muon(0)->genParticlesSize() == 1  &&  muJet->muon(1)->genParticlesSize() == 1) {
         const reco::GenParticle *mu0 = muJet->muon(0)->genParticle();
         const reco::GenParticle *mu1 = muJet->muon(1)->genParticle();

         double total_energy = mu0->energy() + mu1->energy();
         double total_px = mu0->px() + mu1->px();
         double total_py = mu0->py() + mu1->py();
         double total_pz = mu0->pz() + mu1->pz();
         m_lowdimuon_genmass = sqrt(total_energy*total_energy - total_px*total_px - total_py*total_py - total_pz*total_pz);
      }

      m_lowdimuon_mass = muJet->mass();
      m_lowdimuon_pt = muJet->pt();
      m_lowdimuon_eta = muJet->eta();
      m_lowdimuon_phi = muJet->phi();
      m_lowdimuon_dr = muJet->dRmax();
      if (muJet->daughter(0)->charge() > 0) {
         m_lowdimuon_pluspx = muJet->daughter(0)->px();
         m_lowdimuon_pluspy = muJet->daughter(0)->py();
         m_lowdimuon_pluspz = muJet->daughter(0)->pz();
         m_lowdimuon_minuspx = muJet->daughter(1)->px();
         m_lowdimuon_minuspy = muJet->daughter(1)->py();
         m_lowdimuon_minuspz = muJet->daughter(1)->pz();
      }
      else {
         m_lowdimuon_pluspx = muJet->daughter(1)->px();
         m_lowdimuon_pluspy = muJet->daughter(1)->py();
         m_lowdimuon_pluspz = muJet->daughter(1)->pz();
         m_lowdimuon_minuspx = muJet->daughter(0)->px();
         m_lowdimuon_minuspy = muJet->daughter(0)->py();
         m_lowdimuon_minuspz = muJet->daughter(0)->pz();
      }
      m_lowdimuon_vprob = -1000.;
      m_lowdimuon_vx = -1000.;
      m_lowdimuon_vy = -1000.;
      m_lowdimuon_vz = -1000.;

      GlobalVector complus;
      if (muJet->daughter(0)->charge() > 0) complus = muJet->daughterCOM(0);
      else complus = muJet->daughterCOM(1);
      m_lowdimuon_compluspx = complus.x();
      m_lowdimuon_compluspy = complus.y();
      m_lowdimuon_compluspz = complus.z();

      if (muJet->daughter(0)->charge() > 0) complus = muJet->daughterCOMrot(0);
      else complus = muJet->daughterCOMrot(1);
      m_lowdimuon_comrotpluspx = complus.x();
      m_lowdimuon_comrotpluspy = complus.y();
      m_lowdimuon_comrotpluspz = complus.z();

      // replace all values with vertex-updated values if vertex fitting succeeded
      if (muJet->vertexValid()) {
         m_lowdimuon_mass = muJet->vertexMass();
         m_lowdimuon_pt = muJet->vertexMomentum().perp();
         m_lowdimuon_eta = muJet->vertexMomentum().eta();
         m_lowdimuon_phi = muJet->vertexMomentum().phi();
         m_lowdimuon_dr = muJet->dRmax(true);
         if (muJet->daughter(0)->charge() > 0) {
            m_lowdimuon_pluspx = muJet->vertexMomentum(0).x();
            m_lowdimuon_pluspy = muJet->vertexMomentum(0).y();
            m_lowdimuon_pluspz = muJet->vertexMomentum(0).z();
            m_lowdimuon_minuspx = muJet->vertexMomentum(1).x();
            m_lowdimuon_minuspy = muJet->vertexMomentum(1).y();
            m_lowdimuon_minuspz = muJet->vertexMomentum(1).z();
         }
         else {
            m_lowdimuon_pluspx = muJet->vertexMomentum(1).x();
            m_lowdimuon_pluspy = muJet->vertexMomentum(1).y();
            m_lowdimuon_pluspz = muJet->vertexMomentum(1).z();
            m_lowdimuon_minuspx = muJet->vertexMomentum(0).x();
            m_lowdimuon_minuspy = muJet->vertexMomentum(0).y();
            m_lowdimuon_minuspz = muJet->vertexMomentum(0).z();
         }
         m_lowdimuon_vprob = muJet->vertexProb();

         if (closestPrimaryVertex != primaryVertices->end()) {
            m_lowdimuon_vx = muJet->vertexPoint().x() - closestPrimaryVertex->x();
            m_lowdimuon_vy = muJet->vertexPoint().y() - closestPrimaryVertex->y();
            m_lowdimuon_vz = muJet->vertexPoint().z() - closestPrimaryVertex->z();
         }

         GlobalVector complus;
         if (muJet->daughter(0)->charge() > 0) complus = muJet->daughterCOM(0, true);
         else complus = muJet->daughterCOM(1, true);
         m_lowdimuon_compluspx = complus.x();
         m_lowdimuon_compluspy = complus.y();
         m_lowdimuon_compluspz = complus.z();

         if (muJet->daughter(0)->charge() > 0) complus = muJet->daughterCOMrot(0, true);
         else complus = muJet->daughterCOMrot(1, true);
         m_lowdimuon_comrotpluspx = complus.x();
         m_lowdimuon_comrotpluspy = complus.y();
         m_lowdimuon_comrotpluspz = complus.z();

      } // end of replacements with quantities measured at the vertex

      m_lowdimuon_iso = muJet->centralTrackIsolation();

      std::vector<reco::MuonChamberMatch> plusmatches, minusmatches;
      if (muJet->daughter(0)->charge() > 0) {
         m_lowdimuon_plusmatches = muJet->muon(0)->numberOfMatches(reco::Muon::SegmentAndTrackArbitration);
         m_lowdimuon_minusmatches = muJet->muon(1)->numberOfMatches(reco::Muon::SegmentAndTrackArbitration);
         plusmatches = muJet->muon(0)->matches();
         minusmatches = muJet->muon(1)->matches();
         m_lowdimuon_plushits = muJet->muon(0)->innerTrack()->numberOfValidHits();
         m_lowdimuon_minushits = muJet->muon(1)->innerTrack()->numberOfValidHits();
         m_lowdimuon_plusnormchi2 = muJet->muon(0)->innerTrack()->normalizedChi2();
         m_lowdimuon_minusnormchi2 = muJet->muon(1)->innerTrack()->normalizedChi2();
      }
      else {
         m_lowdimuon_plusmatches = muJet->muon(1)->numberOfMatches(reco::Muon::SegmentAndTrackArbitration);
         m_lowdimuon_minusmatches = muJet->muon(0)->numberOfMatches(reco::Muon::SegmentAndTrackArbitration);
         plusmatches = muJet->muon(1)->matches();
         minusmatches = muJet->muon(0)->matches();
         m_lowdimuon_plushits = muJet->muon(1)->innerTrack()->numberOfValidHits();
         m_lowdimuon_minushits = muJet->muon(0)->innerTrack()->numberOfValidHits();
         m_lowdimuon_plusnormchi2 = muJet->muon(1)->innerTrack()->normalizedChi2();
         m_lowdimuon_minusnormchi2 = muJet->muon(0)->innerTrack()->normalizedChi2();
      }
 
      m_lowdimuon_plusst1x = -1000.;
      m_lowdimuon_plusst1xsig = -1000.;
      m_lowdimuon_plusst1dxdz = -1000.;
      m_lowdimuon_plusst1dxdzsig = -1000.;
      for (std::vector<reco::MuonChamberMatch>::const_iterator chamberMatch = plusmatches.begin();  chamberMatch != plusmatches.end();  ++chamberMatch) {
         for (std::vector<reco::MuonSegmentMatch>::const_iterator segmentMatch = chamberMatch->segmentMatches.begin();  segmentMatch != chamberMatch->segmentMatches.end();  ++segmentMatch) {
            if (segmentMatch->isMask(reco::MuonSegmentMatch::BestInChamberByDR) &&
                segmentMatch->isMask(reco::MuonSegmentMatch::BelongsToTrackByDR)) {
               if (chamberMatch->station() == 1) {
                  m_lowdimuon_plusst1x = (segmentMatch->x - chamberMatch->x);
                  m_lowdimuon_plusst1xsig = (segmentMatch->x - chamberMatch->x)/sqrt(pow(segmentMatch->xErr, 2) + pow(chamberMatch->xErr, 2));
                  m_lowdimuon_plusst1dxdz = (segmentMatch->dXdZ - chamberMatch->dXdZ);
                  m_lowdimuon_plusst1dxdzsig = (segmentMatch->dXdZ - chamberMatch->dXdZ)/sqrt(pow(segmentMatch->dXdZErr, 2) + pow(chamberMatch->dXdZErr, 2));
               }
            }
         }
      }

      m_lowdimuon_minusst1x = -1000.;
      m_lowdimuon_minusst1xsig = -1000.;
      m_lowdimuon_minusst1dxdz = -1000.;
      m_lowdimuon_minusst1dxdzsig = -1000.;
      for (std::vector<reco::MuonChamberMatch>::const_iterator chamberMatch = minusmatches.begin();  chamberMatch != minusmatches.end();  ++chamberMatch) {
         for (std::vector<reco::MuonSegmentMatch>::const_iterator segmentMatch = chamberMatch->segmentMatches.begin();  segmentMatch != chamberMatch->segmentMatches.end();  ++segmentMatch) {
            if (segmentMatch->isMask(reco::MuonSegmentMatch::BestInChamberByDR) &&
                segmentMatch->isMask(reco::MuonSegmentMatch::BelongsToTrackByDR)) {
               if (chamberMatch->station() == 1) {
                  m_lowdimuon_minusst1x = (segmentMatch->x - chamberMatch->x);
                  m_lowdimuon_minusst1xsig = (segmentMatch->x - chamberMatch->x)/sqrt(pow(segmentMatch->xErr, 2) + pow(chamberMatch->xErr, 2));
                  m_lowdimuon_minusst1dxdz = (segmentMatch->dXdZ - chamberMatch->dXdZ);
                  m_lowdimuon_minusst1dxdzsig = (segmentMatch->dXdZ - chamberMatch->dXdZ)/sqrt(pow(segmentMatch->dXdZErr, 2) + pow(chamberMatch->dXdZErr, 2));
               }
            }
         }
      }

      if (m_lowdimuon_pt > 80.) m_highdimuon->Fill();
      else m_lowdimuon->Fill();
   }

   ////////////////////////////////////////////////////////// quadmuon
   else if (muJets->size() == 1  &&  (*muJets)[0].numberOfDaughters() >= 4) {
      pat::MultiMuonCollection::const_iterator muJet = muJets->begin();

      m_quadmuon_containstrig = 0;
      for (std::vector<pat::MuonCollection::const_iterator>::const_iterator iter = muontrig.begin();  iter != muontrig.end();  ++iter) {
         for (unsigned int i = 0;  i < muJet->numberOfDaughters();  i++) {
            if (muJet->muon(i)->innerTrack().isAvailable()  &&  (*iter)->innerTrack().isAvailable()  &&
                muJet->sameTrack(&*(muJet->muon(i)->innerTrack()), &*((*iter)->innerTrack()))) {
               m_quadmuon_containstrig = 1;
            }
         }
      }

      m_quadmuon_muons = muJet->numberOfDaughters();
      m_quadmuon_orphans = orphans->size();
      m_quadmuon_charge = muJet->charge();
      m_quadmuon_mass = muJet->mass();
      m_quadmuon_massa = -1.;
      m_quadmuon_massb = -1.;
      m_quadmuon_massc = -1.;
      m_quadmuon_massd = -1.;
      m_quadmuon_masse = -1.;
      m_quadmuon_massf = -1.;
      std::vector<double> masses = muJet->consistentPairMasses();
      std::random_shuffle(masses.begin(), masses.end());
      if (masses.size() >= 1) m_quadmuon_massa = masses[0];
      if (masses.size() >= 2) m_quadmuon_massb = masses[1];
      if (masses.size() >= 3) m_quadmuon_massc = masses[2];
      if (masses.size() >= 4) m_quadmuon_massd = masses[3];
      if (masses.size() >= 5) m_quadmuon_masse = masses[4];
      if (masses.size() >= 6) m_quadmuon_massf = masses[5];
      m_quadmuon_pt = muJet->pt();
      m_quadmuon_eta = muJet->eta();
      m_quadmuon_phi = muJet->phi();
      m_quadmuon_dr = muJet->dRmax();
      m_quadmuon_vchi2 = -1.;
      m_quadmuon_vndof = -1;
      m_quadmuon_vprob = -1.;
      m_quadmuon_lxy = -1000.;
      m_quadmuon_lxyz = -1000.;
      m_quadmuon_iso = muJet->centralTrackIsolation();

      if (muJet->vertexValid()) {
         m_quadmuon_mass = muJet->vertexMass();
         m_quadmuon_massa = -1.;
         m_quadmuon_massb = -1.;
         m_quadmuon_massc = -1.;
         m_quadmuon_massd = -1.;
         m_quadmuon_masse = -1.;
         m_quadmuon_massf = -1.;
         std::vector<double> masses = muJet->consistentPairMasses(true);
         std::random_shuffle(masses.begin(), masses.end());
         if (masses.size() >= 1) m_quadmuon_massa = masses[0];
         if (masses.size() >= 2) m_quadmuon_massb = masses[1];
         if (masses.size() >= 3) m_quadmuon_massc = masses[2];
         if (masses.size() >= 4) m_quadmuon_massd = masses[3];
         if (masses.size() >= 5) m_quadmuon_masse = masses[4];
         if (masses.size() >= 6) m_quadmuon_massf = masses[5];
         m_quadmuon_pt = muJet->vertexMomentum().perp();
         m_quadmuon_eta = muJet->vertexMomentum().eta();
         m_quadmuon_phi = muJet->vertexMomentum().phi();
         m_quadmuon_dr = muJet->dRmax(true);
         m_quadmuon_vchi2 = muJet->vertexChi2();
         m_quadmuon_vndof = muJet->vertexNdof();
         m_quadmuon_vprob = muJet->vertexProb();
         if (closestPrimaryVertex != primaryVertices->end()) {
            m_quadmuon_lxy = muJet->lxy(GlobalPoint(closestPrimaryVertex->x(), closestPrimaryVertex->y(), closestPrimaryVertex->z()));
            m_quadmuon_lxyz = muJet->lxyz(GlobalPoint(closestPrimaryVertex->x(), closestPrimaryVertex->y(), closestPrimaryVertex->z()));
         }
      }

      m_quadmuon->Fill();
   }

   ////////////////////////////////////////////////////////// dimuorphan
   else if (muJets->size() == 1  &&  (*muJets)[0].numberOfDaughters() == 2  &&  orphans->size() == 1) {
      pat::MultiMuonCollection::const_iterator muJet = muJets->begin();
      pat::MuonCollection::const_iterator orphan = orphans->begin();

      m_dimuorphan_deltaphi = muJet->phi() - orphan->phi();
      while (m_dimuorphan_deltaphi > M_PI) m_dimuorphan_deltaphi -= 2.*M_PI;
      while (m_dimuorphan_deltaphi < -M_PI) m_dimuorphan_deltaphi += 2.*M_PI;

      m_dimuorphan_orphanpt = orphan->pt();
      m_dimuorphan_orphaneta = orphan->eta();
      m_dimuorphan_orphanphi = orphan->phi();
      m_dimuorphan_orphanisglobal = (orphan->isGlobalMuon() ? 1 : 0);
      m_dimuorphan_orphanmatches = orphan->numberOfMatches(reco::Muon::SegmentAndTrackArbitration);
      m_dimuorphan_orphanhits = (orphan->innerTrack().isAvailable() ? orphan->innerTrack()->numberOfValidHits(): -1);
      m_dimuorphan_orphanchi2 = (orphan->innerTrack().isAvailable() ? orphan->innerTrack()->normalizedChi2(): -1.);
      m_dimuorphan_orphanst1x = -1000.;
      m_dimuorphan_orphanst1xsig = -1000.;
      m_dimuorphan_orphanst1dxdz = -1000.;
      m_dimuorphan_orphanst1dxdzsig = -1000.;
      for (std::vector<reco::MuonChamberMatch>::const_iterator chamberMatch = orphan->matches().begin();  chamberMatch != orphan->matches().end();  ++chamberMatch) {
         for (std::vector<reco::MuonSegmentMatch>::const_iterator segmentMatch = chamberMatch->segmentMatches.begin();  segmentMatch != chamberMatch->segmentMatches.end();  ++segmentMatch) {
            if (segmentMatch->isMask(reco::MuonSegmentMatch::BestInChamberByDR) &&
                segmentMatch->isMask(reco::MuonSegmentMatch::BelongsToTrackByDR)) {
               if (chamberMatch->station() == 1) {
                  m_dimuorphan_orphanst1x = (segmentMatch->x - chamberMatch->x);
                  m_dimuorphan_orphanst1xsig = (segmentMatch->x - chamberMatch->x)/sqrt(pow(segmentMatch->xErr, 2) + pow(chamberMatch->xErr, 2));
                  m_dimuorphan_orphanst1dxdz = (segmentMatch->dXdZ - chamberMatch->dXdZ);
                  m_dimuorphan_orphanst1dxdzsig = (segmentMatch->dXdZ - chamberMatch->dXdZ)/sqrt(pow(segmentMatch->dXdZErr, 2) + pow(chamberMatch->dXdZErr, 2));
               }
            }
         }
      }

      m_dimuorphan_containstrig = 0;
      for (std::vector<pat::MuonCollection::const_iterator>::const_iterator iter = muontrig.begin();  iter != muontrig.end();  ++iter) {
         if (orphan->innerTrack().isAvailable()  &&  (*iter)->innerTrack().isAvailable()  &&
             muJet->sameTrack(&*(orphan->innerTrack()), &*((*iter)->innerTrack()))) {
            m_dimuorphan_containstrig = 1;
         }
      }
         
      m_dimuorphan_mass = muJet->mass();
      m_dimuorphan_pt = muJet->pt();
      m_dimuorphan_eta = muJet->eta();
      m_dimuorphan_phi = muJet->phi();
      m_dimuorphan_dr = muJet->dRmax();
      m_dimuorphan_vprob = -1.;
      m_dimuorphan_lxy = -1000.;
      m_dimuorphan_lxyz = -1000.;
      m_dimuorphan_iso = muJet->centralTrackIsolation();

      if (muJet->vertexValid()) {
         m_dimuorphan_deltaphi = muJet->vertexMomentum().phi() - orphan->phi();
         while (m_dimuorphan_deltaphi > M_PI) m_dimuorphan_deltaphi -= 2.*M_PI;
         while (m_dimuorphan_deltaphi < -M_PI) m_dimuorphan_deltaphi += 2.*M_PI;

         m_dimuorphan_mass = muJet->vertexMass();
         m_dimuorphan_pt = muJet->vertexMomentum().perp();
         m_dimuorphan_eta = muJet->vertexMomentum().eta();
         m_dimuorphan_phi = muJet->vertexMomentum().phi();
         m_dimuorphan_dr = muJet->dRmax(true);
         m_dimuorphan_vprob = muJet->vertexProb();
         if (closestPrimaryVertex != primaryVertices->end()) {
            m_dimuorphan_lxy = muJet->lxy(GlobalPoint(closestPrimaryVertex->x(), closestPrimaryVertex->y(), closestPrimaryVertex->z()));
            m_dimuorphan_lxyz = muJet->lxyz(GlobalPoint(closestPrimaryVertex->x(), closestPrimaryVertex->y(), closestPrimaryVertex->z()));
         }
      }

      m_dimuorphan->Fill();
   }

   ////////////////////////////////////////////////////////// dimudimu
   else if (muJets->size() == 2) {
      const pat::MultiMuon *muJet0 = &((*muJets)[0]);
      const pat::MultiMuon *muJet1 = &((*muJets)[1]);

      bool muJet0_canBeC = false;
      bool muJet1_canBeC = false;
      for (std::vector<pat::MuonCollection::const_iterator>::const_iterator iter = muontrig.begin();  iter != muontrig.end();  ++iter) {
         if (muJet0->muon(0)->innerTrack().isAvailable()  &&  (*iter)->innerTrack().isAvailable()  &&
             muJet0->sameTrack(&*(muJet0->muon(0)->innerTrack()), &*((*iter)->innerTrack()))) {
            muJet0_canBeC = true;
         }
         if (muJet0->muon(1)->innerTrack().isAvailable()  &&  (*iter)->innerTrack().isAvailable()  &&
             muJet0->sameTrack(&*(muJet0->muon(1)->innerTrack()), &*((*iter)->innerTrack()))) {
            muJet0_canBeC = true;
         }

         if (muJet1->muon(0)->innerTrack().isAvailable()  &&  (*iter)->innerTrack().isAvailable()  &&
             muJet1->sameTrack(&*(muJet1->muon(0)->innerTrack()), &*((*iter)->innerTrack()))) {
            muJet1_canBeC = true;
         }
         if (muJet1->muon(1)->innerTrack().isAvailable()  &&  (*iter)->innerTrack().isAvailable()  &&
             muJet1->sameTrack(&*(muJet1->muon(1)->innerTrack()), &*((*iter)->innerTrack()))) {
            muJet1_canBeC = true;
         }
      }

      const pat::MultiMuon *muJetC = NULL;
      const pat::MultiMuon *muJetF = NULL;
      if (muJet0_canBeC  &&  muJet1_canBeC) {
         if (m_trandom3.Integer(2) == 0) {
            muJetC = muJet0;
            muJetF = muJet1;
         }
         else {
            muJetC = muJet1;
            muJetF = muJet0;
         }
      }
      else if (muJet0_canBeC) {
         muJetC = muJet0;
         muJetF = muJet1;
      }
      else if (muJet1_canBeC) {
         muJetC = muJet1;
         muJetF = muJet0;
      }

      if (muJetC != NULL  &&  muJetF != NULL) {
         m_dimudimu_orphans = orphans->size();

         double total_energy = muJetC->energy() + muJetF->energy();
         double total_px = muJetC->px() + muJetF->px();
         double total_py = muJetC->py() + muJetF->py();
         double total_pz = muJetC->pz() + muJetF->pz();

         m_dimudimu_wholemass = sqrt(total_energy*total_energy - total_px*total_px - total_py*total_py - total_pz*total_pz);
         m_dimudimu_wholept = sqrt(total_px*total_px + total_py*total_py);

         m_dimudimu_deltaphi = muJetC->phi() - muJetF->phi();
         while (m_dimudimu_deltaphi > M_PI) m_dimudimu_deltaphi -= 2.*M_PI;
         while (m_dimudimu_deltaphi < -M_PI) m_dimudimu_deltaphi += 2.*M_PI;

         m_dimudimu_massC = muJetC->mass();
         m_dimudimu_ptC = muJetC->pt();
         m_dimudimu_etaC = muJetC->eta();
         m_dimudimu_phiC = muJetC->phi();
         m_dimudimu_drC = muJetC->dRmax();
         m_dimudimu_vprobC = -1.;
         m_dimudimu_lxyC = -1000.;
         m_dimudimu_lxyzC = -1000.;
         m_dimudimu_isoC = muJetC->centralTrackIsolation();

         m_dimudimu_massF = muJetF->mass();
         m_dimudimu_ptF = muJetF->pt();
         m_dimudimu_etaF = muJetF->eta();
         m_dimudimu_phiF = muJetF->phi();
         m_dimudimu_drF = muJetF->dRmax();
         m_dimudimu_vprobF = -1.;
         m_dimudimu_lxyF = -1000.;
         m_dimudimu_lxyzF = -1000.;
         m_dimudimu_isoF = muJetF->centralTrackIsolation();

         if (muJetC->vertexValid() && muJetF->vertexValid()) {
            total_energy = muJetC->vertexP4().energy() + muJetF->vertexP4().energy();
            total_px = muJetC->vertexMomentum().x() + muJetF->vertexMomentum().x();
            total_py = muJetC->vertexMomentum().y() + muJetF->vertexMomentum().y();
            total_pz = muJetC->vertexMomentum().z() + muJetF->vertexMomentum().z();

            m_dimudimu_wholemass = sqrt(total_energy*total_energy - total_px*total_px - total_py*total_py - total_pz*total_pz);
            m_dimudimu_wholept = sqrt(total_px*total_px + total_py*total_py);

            m_dimudimu_deltaphi = muJetC->vertexMomentum().phi() - muJetF->vertexMomentum().phi();
            while (m_dimudimu_deltaphi > M_PI) m_dimudimu_deltaphi -= 2.*M_PI;
            while (m_dimudimu_deltaphi < -M_PI) m_dimudimu_deltaphi += 2.*M_PI;

            m_dimudimu_massC = muJetC->vertexMass();
            m_dimudimu_ptC = muJetC->vertexMomentum().perp();
            m_dimudimu_etaC = muJetC->vertexMomentum().eta();
            m_dimudimu_phiC = muJetC->vertexMomentum().phi();
            m_dimudimu_drC = muJetC->dRmax(true);
            m_dimudimu_vprobC = muJetC->vertexProb();
            if (closestPrimaryVertex != primaryVertices->end()) {
               m_dimudimu_lxyC = muJetC->lxy(GlobalPoint(closestPrimaryVertex->x(), closestPrimaryVertex->y(), closestPrimaryVertex->z()));
               m_dimudimu_lxyzC = muJetC->lxyz(GlobalPoint(closestPrimaryVertex->x(), closestPrimaryVertex->y(), closestPrimaryVertex->z()));
            }

            m_dimudimu_massF = muJetF->vertexMass();
            m_dimudimu_ptF = muJetF->vertexMomentum().perp();
            m_dimudimu_etaF = muJetF->vertexMomentum().eta();
            m_dimudimu_phiF = muJetF->vertexMomentum().phi();
            m_dimudimu_drF = muJetF->dRmax(true);
            m_dimudimu_vprobF = muJetF->vertexProb();
            if (closestPrimaryVertex != primaryVertices->end()) {
               m_dimudimu_lxyF = muJetF->lxy(GlobalPoint(closestPrimaryVertex->x(), closestPrimaryVertex->y(), closestPrimaryVertex->z()));
               m_dimudimu_lxyzF = muJetF->lxyz(GlobalPoint(closestPrimaryVertex->x(), closestPrimaryVertex->y(), closestPrimaryVertex->z()));
            }
         }

         if (muJetC->numberOfDaughters() == 2  &&  muJetF->numberOfDaughters() == 2) {
            m_dimudimu->Fill();
         }
         else {
            m_twomujets_muonsC = muJetC->numberOfDaughters();
            m_twomujets_muonsF = muJetF->numberOfDaughters();

            m_twomujets_massCa = -1.;
            m_twomujets_massCb = -1.;
            m_twomujets_massCc = -1.;
            m_twomujets_massFa = -1.;
            m_twomujets_massFb = -1.;
            m_twomujets_massFc = -1.;
            std::vector<double> massesC = muJetC->consistentPairMasses(muJetC->vertexValid() && muJetF->vertexValid());
            std::random_shuffle(massesC.begin(), massesC.end());
            if (massesC.size() >= 1) m_twomujets_massCa = massesC[0];
            if (massesC.size() >= 2) m_twomujets_massCb = massesC[1];
            if (massesC.size() >= 3) m_twomujets_massCc = massesC[2];
            std::vector<double> massesF = muJetF->consistentPairMasses(muJetC->vertexValid() && muJetF->vertexValid());
            std::random_shuffle(massesF.begin(), massesF.end());
            if (massesF.size() >= 1) m_twomujets_massFa = massesF[0];
            if (massesF.size() >= 2) m_twomujets_massFb = massesF[1];
            if (massesF.size() >= 3) m_twomujets_massFc = massesF[2];

            m_twomujets->Fill();
         }
      }
   }

   ////////////////////////////////////////////////////////// moremujets

   else if (muJets->size() > 2) {
      m_moremujets_mujets = muJets->size();
      m_moremujets_orphans = orphans->size();

      double total_energy = 0.;
      double total_px = 0.;
      double total_py = 0.;
      double total_pz = 0.;
      for (pat::MultiMuonCollection::const_iterator muJet = muJets->begin();  muJet != muJets->end();  ++muJet) {
         total_energy += muJet->energy();
         total_px += muJet->px();
         total_py += muJet->py();
         total_pz += muJet->pz();
      }
      m_moremujets_wholemass = sqrt(total_energy*total_energy - total_px*total_px - total_py*total_py - total_pz*total_pz);
      m_moremujets_wholept = sqrt(total_px*total_px + total_py*total_py);

      m_moremujets_massa = -1.;
      m_moremujets_massb = -1.;
      m_moremujets_massc = -1.;
      m_moremujets_massd = -1.;
      m_moremujets_masse = -1.;
      m_moremujets_massf = -1.;
      std::vector<double> all_masses;
      for (pat::MultiMuonCollection::const_iterator muJet = muJets->begin();  muJet != muJets->end();  ++muJet) {
         std::vector<double> masses = muJet->consistentPairMasses();
         std::random_shuffle(masses.begin(), masses.end());
         for (std::vector<double>::const_iterator mass = masses.begin();  mass != masses.end();  ++mass) {
            all_masses.push_back(*mass);
         }
      }
      if (all_masses.size() >= 1) m_moremujets_massa = all_masses[0];
      if (all_masses.size() >= 2) m_moremujets_massb = all_masses[1];
      if (all_masses.size() >= 3) m_moremujets_massc = all_masses[2];
      if (all_masses.size() >= 4) m_moremujets_massd = all_masses[3];
      if (all_masses.size() >= 5) m_moremujets_masse = all_masses[4];
      if (all_masses.size() >= 6) m_moremujets_massf = all_masses[5];

      m_moremujets_muons1 = -1;
      m_moremujets_containstrig1 = -1;
      m_moremujets_mass1 = -1.;
      m_moremujets_pt1 = -1.;
      m_moremujets_iso1 = -1.;
      m_moremujets_muons2 = -1;
      m_moremujets_containstrig2 = -1;
      m_moremujets_mass2 = -1.;
      m_moremujets_pt2 = -1.;
      m_moremujets_iso2 = -1.;
      m_moremujets_muons3 = -1;
      m_moremujets_containstrig3 = -1;
      m_moremujets_mass3 = -1.;
      m_moremujets_pt3 = -1.;
      m_moremujets_iso3 = -1.;
      m_moremujets_muons4 = -1;
      m_moremujets_containstrig4 = -1;
      m_moremujets_mass4 = -1.;
      m_moremujets_pt4 = -1.;
      m_moremujets_iso4 = -1.;
      m_moremujets_muons5 = -1;
      m_moremujets_containstrig5 = -1;
      m_moremujets_mass5 = -1.;
      m_moremujets_pt5 = -1.;
      m_moremujets_iso5 = -1.;
      m_moremujets_muons6 = -1;
      m_moremujets_containstrig6 = -1;
      m_moremujets_mass6 = -1.;
      m_moremujets_pt6 = -1.;
      m_moremujets_iso6 = -1.;

      if (muJets->size() >= 1) {
         const pat::MultiMuon *muJet = &((*muJets)[0]);
         m_moremujets_muons1 = muJet->numberOfDaughters();
         m_moremujets_containstrig1 = 0;
         for (std::vector<pat::MuonCollection::const_iterator>::const_iterator iter = muontrig.begin();  iter != muontrig.end();  ++iter) {
            for (unsigned int i = 0;  i < muJet->numberOfDaughters();  i++) {
               if (muJet->muon(i)->innerTrack().isAvailable()  &&  (*iter)->innerTrack().isAvailable()  &&
                   muJet->sameTrack(&*(muJet->muon(i)->innerTrack()), &*((*iter)->innerTrack()))) {
                  m_moremujets_containstrig1 = 1;
               }
            }
         }
         m_moremujets_mass1 = muJet->mass();
         m_moremujets_pt1 = muJet->pt();
         m_moremujets_iso1 = muJet->centralTrackIsolation();
      }

      if (muJets->size() >= 2) {
         const pat::MultiMuon *muJet = &((*muJets)[1]);
         m_moremujets_muons2 = muJet->numberOfDaughters();
         m_moremujets_containstrig2 = 0;
         for (std::vector<pat::MuonCollection::const_iterator>::const_iterator iter = muontrig.begin();  iter != muontrig.end();  ++iter) {
            for (unsigned int i = 0;  i < muJet->numberOfDaughters();  i++) {
               if (muJet->muon(i)->innerTrack().isAvailable()  &&  (*iter)->innerTrack().isAvailable()  &&
                   muJet->sameTrack(&*(muJet->muon(i)->innerTrack()), &*((*iter)->innerTrack()))) {
                  m_moremujets_containstrig2 = 1;
               }
            }
         }
         m_moremujets_mass2 = muJet->mass();
         m_moremujets_pt2 = muJet->pt();
         m_moremujets_iso2 = muJet->centralTrackIsolation();
      }

      if (muJets->size() >= 3) {
         const pat::MultiMuon *muJet = &((*muJets)[2]);
         m_moremujets_muons3 = muJet->numberOfDaughters();
         m_moremujets_containstrig3 = 0;
         for (std::vector<pat::MuonCollection::const_iterator>::const_iterator iter = muontrig.begin();  iter != muontrig.end();  ++iter) {
            for (unsigned int i = 0;  i < muJet->numberOfDaughters();  i++) {
               if (muJet->muon(i)->innerTrack().isAvailable()  &&  (*iter)->innerTrack().isAvailable()  &&
                   muJet->sameTrack(&*(muJet->muon(i)->innerTrack()), &*((*iter)->innerTrack()))) {
                  m_moremujets_containstrig3 = 1;
               }
            }
         }
         m_moremujets_mass3 = muJet->mass();
         m_moremujets_pt3 = muJet->pt();
         m_moremujets_iso3 = muJet->centralTrackIsolation();
      }

      if (muJets->size() >= 4) {
         const pat::MultiMuon *muJet = &((*muJets)[3]);
         m_moremujets_muons4 = muJet->numberOfDaughters();
         m_moremujets_containstrig4 = 0;
         for (std::vector<pat::MuonCollection::const_iterator>::const_iterator iter = muontrig.begin();  iter != muontrig.end();  ++iter) {
            for (unsigned int i = 0;  i < muJet->numberOfDaughters();  i++) {
               if (muJet->muon(i)->innerTrack().isAvailable()  &&  (*iter)->innerTrack().isAvailable()  &&
                   muJet->sameTrack(&*(muJet->muon(i)->innerTrack()), &*((*iter)->innerTrack()))) {
                  m_moremujets_containstrig4 = 1;
               }
            }
         }
         m_moremujets_mass4 = muJet->mass();
         m_moremujets_pt4 = muJet->pt();
         m_moremujets_iso4 = muJet->centralTrackIsolation();
      }

      if (muJets->size() >= 5) {
         const pat::MultiMuon *muJet = &((*muJets)[4]);
         m_moremujets_muons5 = muJet->numberOfDaughters();
         m_moremujets_containstrig5 = 0;
         for (std::vector<pat::MuonCollection::const_iterator>::const_iterator iter = muontrig.begin();  iter != muontrig.end();  ++iter) {
            for (unsigned int i = 0;  i < muJet->numberOfDaughters();  i++) {
               if (muJet->muon(i)->innerTrack().isAvailable()  &&  (*iter)->innerTrack().isAvailable()  &&
                   muJet->sameTrack(&*(muJet->muon(i)->innerTrack()), &*((*iter)->innerTrack()))) {
                  m_moremujets_containstrig5 = 1;
               }
            }
         }
         m_moremujets_mass5 = muJet->mass();
         m_moremujets_pt5 = muJet->pt();
         m_moremujets_iso5 = muJet->centralTrackIsolation();
      }

      if (muJets->size() >= 6) {
         const pat::MultiMuon *muJet = &((*muJets)[5]);
         m_moremujets_muons6 = muJet->numberOfDaughters();
         m_moremujets_containstrig6 = 0;
         for (std::vector<pat::MuonCollection::const_iterator>::const_iterator iter = muontrig.begin();  iter != muontrig.end();  ++iter) {
            for (unsigned int i = 0;  i < muJet->numberOfDaughters();  i++) {
               if (muJet->muon(i)->innerTrack().isAvailable()  &&  (*iter)->innerTrack().isAvailable()  &&
                   muJet->sameTrack(&*(muJet->muon(i)->innerTrack()), &*((*iter)->innerTrack()))) {
                  m_moremujets_containstrig6 = 1;
               }
            }
         }
         m_moremujets_mass6 = muJet->mass();
         m_moremujets_pt6 = muJet->pt();
         m_moremujets_iso6 = muJet->centralTrackIsolation();
      }

      m_moremujets->Fill();
   }

   ////////////////////////////////////////////////////////// mujetplustracks
   if (muJetPlusTracks->size() == 1  &&  ((*muJetPlusTracks)[0].userInt("numReal") == 2  ||  (*muJetPlusTracks)[0].userInt("numReal") == 3)  &&  (*muJetPlusTracks)[0].userInt("numFakes") > 0) {
      const pat::MultiMuon *muJet = &((*muJetPlusTracks)[0]);

      m_mujetplustracks_containstrig = 0;
      for (std::vector<pat::MuonCollection::const_iterator>::const_iterator iter = muontrig.begin();  iter != muontrig.end();  ++iter) {
         if (muJet->muon(0)->innerTrack().isAvailable()  &&  (*iter)->innerTrack().isAvailable()  &&
             muJet->sameTrack(&*(muJet->muon(0)->innerTrack()), &*((*iter)->innerTrack()))) {
            m_mujetplustracks_containstrig = 1;
         }
         if (muJet->muon(1)->innerTrack().isAvailable()  &&  (*iter)->innerTrack().isAvailable()  &&
             muJet->sameTrack(&*(muJet->muon(1)->innerTrack()), &*((*iter)->innerTrack()))) {
            m_mujetplustracks_containstrig = 1;
         }
      }

      m_mujetplustracks_muons = muJet->userInt("numReal");
      m_mujetplustracks_fakes = muJet->userInt("numFakes");
      m_mujetplustracks_mass = muJet->mass();
      m_mujetplustracks_charge = muJet->charge();
      m_mujetplustracks_massa = -1.;
      m_mujetplustracks_massb = -1.;
      std::vector<double> masses = muJet->consistentPairMasses();
      std::random_shuffle(masses.begin(), masses.end());
      if (masses.size() >= 1) m_mujetplustracks_massa = masses[0];
      if (masses.size() >= 2) m_mujetplustracks_massb = masses[1];

      m_mujetplustracks_pt = muJet->pt();
      m_mujetplustracks_eta = muJet->eta();
      m_mujetplustracks_phi = muJet->phi();
      m_mujetplustracks_dr = muJet->dRmax();
      m_mujetplustracks_vprob = -1.;
      m_mujetplustracks_lxy = -1000.;
      m_mujetplustracks_lxyz = -1000.;
      m_mujetplustracks_iso = muJet->centralTrackIsolation();

      m_mujetplustracks_mupt1 = -1.;
      m_mujetplustracks_mupt2 = -1.;
      m_mujetplustracks_mupt3 = -1.;
      m_mujetplustracks_fakept1 = -1.;
      m_mujetplustracks_fakept2 = -1.;
      for (unsigned int i = 0;  i < muJet->numberOfDaughters();  i++) {
         if (muJet->muon(i)->matches().size() > 0) {
            if (muJet->muon(i)->pt() > m_mujetplustracks_mupt1) {
               m_mujetplustracks_mupt3 = m_mujetplustracks_mupt2;
               m_mujetplustracks_mupt2 = m_mujetplustracks_mupt1;
               m_mujetplustracks_mupt1 = muJet->muon(i)->pt();
            }
            else if (muJet->muon(i)->pt() > m_mujetplustracks_mupt2) {
               m_mujetplustracks_mupt3 = m_mujetplustracks_mupt2;
               m_mujetplustracks_mupt2 = muJet->muon(i)->pt();
            }
            else if (muJet->muon(i)->pt() > m_mujetplustracks_mupt3) {
               m_mujetplustracks_mupt3 = muJet->muon(i)->pt();
            }
         }
         else {
            if (muJet->muon(i)->pt() > m_mujetplustracks_fakept1) {
               m_mujetplustracks_fakept2 = m_mujetplustracks_fakept1;
               m_mujetplustracks_fakept1 = muJet->muon(i)->pt();
            }
            else if (muJet->muon(i)->pt() > m_mujetplustracks_fakept1) {
               m_mujetplustracks_fakept2 = muJet->muon(i)->pt();
            }
         }
      }

      if (muJet->vertexValid()) {
         m_mujetplustracks_mass = muJet->vertexMass();
         m_mujetplustracks_massa = -1.;
         m_mujetplustracks_massb = -1.;
         std::vector<double> masses = muJet->consistentPairMasses(true);
         std::random_shuffle(masses.begin(), masses.end());
         if (masses.size() >= 1) m_mujetplustracks_massa = masses[0];
         if (masses.size() >= 2) m_mujetplustracks_massb = masses[1];

         m_mujetplustracks_pt = muJet->vertexMomentum().perp();
         m_mujetplustracks_eta = muJet->vertexMomentum().eta();
         m_mujetplustracks_phi = muJet->vertexMomentum().phi();
         m_mujetplustracks_dr = muJet->dRmax(true);
         m_mujetplustracks_vprob = muJet->vertexProb();
         if (closestPrimaryVertex != primaryVertices->end()) {
            m_mujetplustracks_lxy = muJet->lxy(GlobalPoint(closestPrimaryVertex->x(), closestPrimaryVertex->y(), closestPrimaryVertex->z()));
            m_mujetplustracks_lxyz = muJet->lxyz(GlobalPoint(closestPrimaryVertex->x(), closestPrimaryVertex->y(), closestPrimaryVertex->z()));
         }

         m_mujetplustracks_mupt1 = -1.;
         m_mujetplustracks_mupt2 = -1.;
         m_mujetplustracks_mupt3 = -1.;
         m_mujetplustracks_fakept1 = -1.;
         m_mujetplustracks_fakept2 = -1.;
         for (unsigned int i = 0;  i < muJet->numberOfDaughters();  i++) {
            if (muJet->muon(i)->matches().size() > 0) {
               if (muJet->vertexMomentum(i).perp() > m_mujetplustracks_mupt1) {
                  m_mujetplustracks_mupt3 = m_mujetplustracks_mupt2;
                  m_mujetplustracks_mupt2 = m_mujetplustracks_mupt1;
                  m_mujetplustracks_mupt1 = muJet->vertexMomentum(i).perp();
               }
               else if (muJet->vertexMomentum(i).perp() > m_mujetplustracks_mupt2) {
                  m_mujetplustracks_mupt3 = m_mujetplustracks_mupt2;
                  m_mujetplustracks_mupt2 = muJet->vertexMomentum(i).perp();
               }
               else if (muJet->vertexMomentum(i).perp() > m_mujetplustracks_mupt3) {
                  m_mujetplustracks_mupt3 = muJet->vertexMomentum(i).perp();
               }
            }
            else {
               if (muJet->vertexMomentum(i).perp() > m_mujetplustracks_fakept1) {
                  m_mujetplustracks_fakept2 = m_mujetplustracks_fakept1;
                  m_mujetplustracks_fakept1 = muJet->vertexMomentum(i).perp();
               }
               else if (muJet->vertexMomentum(i).perp() > m_mujetplustracks_fakept1) {
                  m_mujetplustracks_fakept2 = muJet->vertexMomentum(i).perp();
               }
            }
         }
      }

      m_mujetplustracks->Fill();
   }

   ////////////////////////////////////////////////////////// displaced
   if (displacedMuJet->size() > 0) {
      m_displaced_mujets = displacedMuJet->size();
      m_displaced_orphans = orphans->size();

      pat::MultiMuonCollection::const_iterator muJet = displacedMuJet->begin();
      if (muJet->vertexValid()) {
         m_displaced_muons = muJet->numberOfDaughters();
         m_displaced_mass = muJet->mass();
         m_displaced_pt = muJet->pt();
         m_displaced_eta = muJet->eta();
         m_displaced_phi = muJet->phi();
         m_displaced_vx = muJet->vertexPoint().x();
         m_displaced_vy = muJet->vertexPoint().y();
         m_displaced_vz = muJet->vertexPoint().z();

         double px = muJet->vertexMomentum().x();
         double py = muJet->vertexMomentum().y();
         double pz = muJet->vertexMomentum().z();
         m_displaced_zclosest = m_displaced_vz - pz * (px*m_displaced_vx + py*m_displaced_vy) / pow(px + py, 2);

         reco::VertexCollection::const_iterator closest = primaryVertices->end();
         for (reco::VertexCollection::const_iterator vertex = primaryVertices->begin();  vertex != primaryVertices->end();  ++vertex) {
            if (vertex->isValid()  &&  !vertex->isFake()  &&  vertex->tracksSize() > 3  &&  fabs(vertex->z()) < 24.) {
               if (closest == primaryVertices->end()  ||  fabs(vertex->z() - m_displaced_zclosest) < fabs(closest->z() - m_displaced_zclosest)) {
                  closest = vertex;
               }
            } // end vertex quality cuts
         } // end loop over primary vertices
      
         m_displaced_zcloseness = -1000.;
         m_displaced_lxy = -1000.;
         m_displaced_lxyz = -1000.;
         if (closest != primaryVertices->end()) {
            m_displaced_zcloseness = m_displaced_zclosest - closest->z();
            m_displaced_lxy = muJet->lxy(GlobalPoint(closest->x(), closest->y(), closest->z()));
            m_displaced_lxyz = muJet->lxyz(GlobalPoint(closest->x(), closest->y(), closest->z()));
         }
         else if (m_getAlternating) {
            m_displaced_zcloseness = m_displaced_zclosest - 0.;
            m_displaced_lxy = muJet->lxy(GlobalPoint(0., 0., 0.));
            m_displaced_lxyz = muJet->lxyz(GlobalPoint(0., 0., 0.));
         }

         m_displaced_iso = muJet->centralTrackIsolation();
         m_displaced_vprob = muJet->vertexProb();
         
         const pat::Muon *plus = NULL;
         const pat::Muon *minus = NULL;
         m_displaced_mupt1 = 0.;
         m_displaced_mupt2 = 0.;
         for (unsigned int i = 0;  i < muJet->numberOfDaughters();  i++) {
            const pat::Muon *muon = muJet->muon(i);
            if (muon->pt() > m_displaced_mupt1) {
               m_displaced_mupt2 = m_displaced_mupt1;
               m_displaced_mupt1 = muon->pt();
            }
            else if (muon->pt() > m_displaced_mupt2) {
               m_displaced_mupt2 = muon->pt();
            }

            if (muon->charge() > 0  &&  muon->innerTrack().isAvailable()  &&  (plus == NULL  ||  muon->pt() > plus->pt())) {
               plus = muon;
            }
            if (muon->charge() < 0  &&  muon->innerTrack().isAvailable()  &&  (minus == NULL  ||  muon->pt() > minus->pt())) {
               minus = muon;
            }
         }

         m_displaced_plusmatches = -1;
         m_displaced_plusst1x = -1000.;
         m_displaced_plusst1xsig = -1000.;
         m_displaced_plusst1dxdz = -1000.;
         m_displaced_plusst1dxdzsig = -1000.;
         m_displaced_plushits = -1;
         m_displaced_plusnormchi2 = -1.;
         m_displaced_minusmatches = -1;
         m_displaced_minusst1x = -1000.;
         m_displaced_minusst1xsig = -1000.;
         m_displaced_minusst1dxdz = -1000.;
         m_displaced_minusst1dxdzsig = -1000.;
         m_displaced_minushits = -1;
         m_displaced_minusnormchi2 = -1.;

         if (plus != NULL  &&  minus != NULL) {
            m_displaced_plusmatches = plus->numberOfMatches(reco::Muon::SegmentAndTrackArbitration);
            m_displaced_plushits = plus->innerTrack()->numberOfValidHits();
            m_displaced_plusnormchi2 = plus->innerTrack()->normalizedChi2();

            std::vector<reco::MuonChamberMatch> plusmatches = plus->matches();
            for (std::vector<reco::MuonChamberMatch>::const_iterator chamberMatch = plusmatches.begin();  chamberMatch != plusmatches.end();  ++chamberMatch) {
               for (std::vector<reco::MuonSegmentMatch>::const_iterator segmentMatch = chamberMatch->segmentMatches.begin();  segmentMatch != chamberMatch->segmentMatches.end();  ++segmentMatch) {
                  if (segmentMatch->isMask(reco::MuonSegmentMatch::BestInChamberByDR) &&
                      segmentMatch->isMask(reco::MuonSegmentMatch::BelongsToTrackByDR)) {
                     if (chamberMatch->station() == 1) {
                        m_displaced_plusst1x = (segmentMatch->x - chamberMatch->x);
                        m_displaced_plusst1xsig = (segmentMatch->x - chamberMatch->x)/sqrt(pow(segmentMatch->xErr, 2) + pow(chamberMatch->xErr, 2));
                        m_displaced_plusst1dxdz = (segmentMatch->dXdZ - chamberMatch->dXdZ);
                        m_displaced_plusst1dxdzsig = (segmentMatch->dXdZ - chamberMatch->dXdZ)/sqrt(pow(segmentMatch->dXdZErr, 2) + pow(chamberMatch->dXdZErr, 2));
                     }
                  }
               }
            }

            m_displaced_minusmatches = minus->numberOfMatches(reco::Muon::SegmentAndTrackArbitration);
            m_displaced_minushits = minus->innerTrack()->numberOfValidHits();
            m_displaced_minusnormchi2 = minus->innerTrack()->normalizedChi2();

            std::vector<reco::MuonChamberMatch> minusmatches = minus->matches();
            for (std::vector<reco::MuonChamberMatch>::const_iterator chamberMatch = minusmatches.begin();  chamberMatch != minusmatches.end();  ++chamberMatch) {
               for (std::vector<reco::MuonSegmentMatch>::const_iterator segmentMatch = chamberMatch->segmentMatches.begin();  segmentMatch != chamberMatch->segmentMatches.end();  ++segmentMatch) {
                  if (segmentMatch->isMask(reco::MuonSegmentMatch::BestInChamberByDR) &&
                      segmentMatch->isMask(reco::MuonSegmentMatch::BelongsToTrackByDR)) {
                     if (chamberMatch->station() == 1) {
                        m_displaced_minusst1x = (segmentMatch->x - chamberMatch->x);
                        m_displaced_minusst1xsig = (segmentMatch->x - chamberMatch->x)/sqrt(pow(segmentMatch->xErr, 2) + pow(chamberMatch->xErr, 2));
                        m_displaced_minusst1dxdz = (segmentMatch->dXdZ - chamberMatch->dXdZ);
                        m_displaced_minusst1dxdzsig = (segmentMatch->dXdZ - chamberMatch->dXdZ)/sqrt(pow(segmentMatch->dXdZErr, 2) + pow(chamberMatch->dXdZErr, 2));
                     }
                  }
               }
            }
         }

         if (sqrt(pow(m_displaced_vx, 2) + pow(m_displaced_vy, 2)) > 0.2) {
            m_displaced->Fill();
         }
      }
   }

   ////////////////////////////////////////////////////////// displacedFake
   if (displacedMuJetFake->size() > 0) {
      m_displaced_mujets = displacedMuJetFake->size();
      m_displaced_orphans = orphans->size();

      pat::MultiMuonCollection::const_iterator muJet = displacedMuJetFake->begin();
      if (muJet->vertexValid()) {
         m_displaced_muons = muJet->numberOfDaughters();
         m_displaced_mass = muJet->mass();
         m_displaced_pt = muJet->pt();
         m_displaced_eta = muJet->eta();
         m_displaced_phi = muJet->phi();
         m_displaced_vx = muJet->vertexPoint().x();
         m_displaced_vy = muJet->vertexPoint().y();
         m_displaced_vz = muJet->vertexPoint().z();

         double px = muJet->vertexMomentum().x();
         double py = muJet->vertexMomentum().y();
         double pz = muJet->vertexMomentum().z();
         m_displaced_zclosest = m_displaced_vz - pz * (px*m_displaced_vx + py*m_displaced_vy) / pow(px + py, 2);
         reco::VertexCollection::const_iterator closest = primaryVertices->end();
         for (reco::VertexCollection::const_iterator vertex = primaryVertices->begin();  vertex != primaryVertices->end();  ++vertex) {
            if (vertex->isValid()  &&  !vertex->isFake()  &&  vertex->tracksSize() > 3  &&  fabs(vertex->z()) < 24.) {
               if (closest == primaryVertices->end()  ||  fabs(vertex->z() - m_displaced_zclosest) < fabs(closest->z() - m_displaced_zclosest)) {
                  closest = vertex;
               }
            } // end vertex quality cuts
         } // end loop over primary vertices
      
         m_displaced_zcloseness = -1000.;
         m_displaced_lxy = -1000.;
         m_displaced_lxyz = -1000.;
         if (closest != primaryVertices->end()) {
            m_displaced_zcloseness = m_displaced_zclosest - closest->z();
            m_displaced_lxy = muJet->lxy(GlobalPoint(closest->x(), closest->y(), closest->z()));
            m_displaced_lxyz = muJet->lxyz(GlobalPoint(closest->x(), closest->y(), closest->z()));
         }
         else if (m_getAlternating) {
            m_displaced_zcloseness = m_displaced_zclosest - 0.;
            m_displaced_lxy = muJet->lxy(GlobalPoint(0., 0., 0.));
            m_displaced_lxyz = muJet->lxyz(GlobalPoint(0., 0., 0.));
         }

         m_displaced_iso = muJet->centralTrackIsolation();
         m_displaced_vprob = muJet->vertexProb();
         
         const pat::Muon *plus = NULL;
         const pat::Muon *minus = NULL;
         m_displaced_mupt1 = 0.;
         m_displaced_mupt2 = 0.;
         for (unsigned int i = 0;  i < muJet->numberOfDaughters();  i++) {
            const pat::Muon *muon = muJet->muon(i);
            if (muon->pt() > m_displaced_mupt1) {
               m_displaced_mupt2 = m_displaced_mupt1;
               m_displaced_mupt1 = muon->pt();
            }
            else if (muon->pt() > m_displaced_mupt2) {
               m_displaced_mupt2 = muon->pt();
            }

            if (muon->charge() > 0  &&  muon->innerTrack().isAvailable()  &&  (plus == NULL  ||  muon->pt() > plus->pt())) {
               plus = muon;
            }
            if (muon->charge() < 0  &&  muon->innerTrack().isAvailable()  &&  (minus == NULL  ||  muon->pt() > minus->pt())) {
               minus = muon;
            }
         }

         m_displaced_plusmatches = -1;
         m_displaced_plusst1x = -1000.;
         m_displaced_plusst1xsig = -1000.;
         m_displaced_plusst1dxdz = -1000.;
         m_displaced_plusst1dxdzsig = -1000.;
         m_displaced_plushits = -1;
         m_displaced_plusnormchi2 = -1.;
         m_displaced_minusmatches = -1;
         m_displaced_minusst1x = -1000.;
         m_displaced_minusst1xsig = -1000.;
         m_displaced_minusst1dxdz = -1000.;
         m_displaced_minusst1dxdzsig = -1000.;
         m_displaced_minushits = -1;
         m_displaced_minusnormchi2 = -1.;

         if (plus != NULL  &&  minus != NULL) {
            m_displaced_plusmatches = plus->numberOfMatches(reco::Muon::SegmentAndTrackArbitration);
            m_displaced_plushits = plus->innerTrack()->numberOfValidHits();
            m_displaced_plusnormchi2 = plus->innerTrack()->normalizedChi2();

            std::vector<reco::MuonChamberMatch> plusmatches = plus->matches();
            for (std::vector<reco::MuonChamberMatch>::const_iterator chamberMatch = plusmatches.begin();  chamberMatch != plusmatches.end();  ++chamberMatch) {
               for (std::vector<reco::MuonSegmentMatch>::const_iterator segmentMatch = chamberMatch->segmentMatches.begin();  segmentMatch != chamberMatch->segmentMatches.end();  ++segmentMatch) {
                  if (segmentMatch->isMask(reco::MuonSegmentMatch::BestInChamberByDR) &&
                      segmentMatch->isMask(reco::MuonSegmentMatch::BelongsToTrackByDR)) {
                     if (chamberMatch->station() == 1) {
                        m_displaced_plusst1x = (segmentMatch->x - chamberMatch->x);
                        m_displaced_plusst1xsig = (segmentMatch->x - chamberMatch->x)/sqrt(pow(segmentMatch->xErr, 2) + pow(chamberMatch->xErr, 2));
                        m_displaced_plusst1dxdz = (segmentMatch->dXdZ - chamberMatch->dXdZ);
                        m_displaced_plusst1dxdzsig = (segmentMatch->dXdZ - chamberMatch->dXdZ)/sqrt(pow(segmentMatch->dXdZErr, 2) + pow(chamberMatch->dXdZErr, 2));
                     }
                  }
               }
            }

            m_displaced_minusmatches = minus->numberOfMatches(reco::Muon::SegmentAndTrackArbitration);
            m_displaced_minushits = minus->innerTrack()->numberOfValidHits();
            m_displaced_minusnormchi2 = minus->innerTrack()->normalizedChi2();

            std::vector<reco::MuonChamberMatch> minusmatches = minus->matches();
            for (std::vector<reco::MuonChamberMatch>::const_iterator chamberMatch = minusmatches.begin();  chamberMatch != minusmatches.end();  ++chamberMatch) {
               for (std::vector<reco::MuonSegmentMatch>::const_iterator segmentMatch = chamberMatch->segmentMatches.begin();  segmentMatch != chamberMatch->segmentMatches.end();  ++segmentMatch) {
                  if (segmentMatch->isMask(reco::MuonSegmentMatch::BestInChamberByDR) &&
                      segmentMatch->isMask(reco::MuonSegmentMatch::BelongsToTrackByDR)) {
                     if (chamberMatch->station() == 1) {
                        m_displaced_minusst1x = (segmentMatch->x - chamberMatch->x);
                        m_displaced_minusst1xsig = (segmentMatch->x - chamberMatch->x)/sqrt(pow(segmentMatch->xErr, 2) + pow(chamberMatch->xErr, 2));
                        m_displaced_minusst1dxdz = (segmentMatch->dXdZ - chamberMatch->dXdZ);
                        m_displaced_minusst1dxdzsig = (segmentMatch->dXdZ - chamberMatch->dXdZ)/sqrt(pow(segmentMatch->dXdZErr, 2) + pow(chamberMatch->dXdZErr, 2));
                     }
                  }
               }
            }
         }

         m_displacedFake_fakes = 0;
         for (unsigned int i = 0;  i < muJet->numberOfDaughters();  i++) {
            if (muJet->muon(i)->numberOfMatches(reco::Muon::SegmentAndTrackArbitration) == 0) {
               m_displacedFake_fakes++;
            }
         }

         if (m_displacedFake_fakes > 0  &&  sqrt(pow(m_displaced_vx, 2) + pow(m_displaced_vy, 2)) > 0.2) {
            m_displacedFake->Fill();
         }
      }
   }
}


// ------------ method called once each job just before starting event loop  ------------
void 
FitNtuple::beginJob()
{
}

// ------------ method called once each job just after ending the event loop  ------------
void 
FitNtuple::endJob() {
}

//define this as a plug-in
DEFINE_FWK_MODULE(FitNtuple);
Revision:	1.14
Committed:	Tue Dec 21 01:16:31 2010 UTC (14 years, 4 months ago) by pivarski
Content type:	text/plain
Branch:	MAIN
CVS Tags:	JP-2010-12-20-a
Changes since 1.13:	+419 -2 lines
Log Message:	add TTrees of displaced mu-jets