//#include "GaudiKernel/AlgFactory.h" //#include "GaudiKernel/IToolSvc.h" //#include "GaudiKernel/AlgTool.h" #include "AnalysisTools/IAnalysisTools.h" #include "StoreGate/StoreGateSvc.h" #include "StoreGate/DataHandle.h" #include "GaudiKernel/ITHistSvc.h" #include "TH1F.h" #include "TH2F.h" #include "egammaEvent/ElectronContainer.h" #include "egammaEvent/PhotonContainer.h" #include "muonEvent/Muon.h" #include "muonEvent/MuonContainer.h" #include "Particle/TrackParticle.h" #include "Particle/TrackParticleContainer.h" #include #include "TLorentzVector.h" #include "TVector3.h" #include "FourMomUtils/P4Helpers.h" #include "CxxUtils/SealCommon.h" #include "CxxUtils/SealDebug.h" #include "JetEvent/JetCollection.h" #include "JetTagEvent/TrackAssociation.h" #include "JetTagInfo/TruthInfo.h" #include "JetTagInfo/IPInfoBase.h" #include "JetTagInfo/IPInfoPlus.h" #include "JetTagInfo/IPTrackInfo.h" #include "JetTagInfo/SVInfoPlus.h" #include "MissingETEvent/MissingEtTruth.h" #include "ParticleEvent/ParticleBaseContainer.h" #include "OSUserHiggs/OSHbb.h" #include "OSSelectors/OSSelectorsEnums.h" #include "EventInfo/EventInfo.h" #include "EventInfo/EventID.h" #include "EventInfo/TriggerInfo.h" #include "EventInfo/EventType.h" #include "InDetBeamSpotService/IBeamCondSvc.h" #include "boost/foreach.hpp" #include #include #include #include #include #include "OSUserHiggs/stringify.h" using std::sort; namespace Analysis { std::string OSHbb::m_btaggerName; OSHbb::OSHbb(const std::string& name, ISvcLocator* pSvcLocator) : OSHBase(name, pSvcLocator), m_selectJetlowpt("Analysis::OSSelectJet"), m_selectJetformet("Analysis::OSSelectJet"), m_kinFitter("Analysis::OSKinFitter"), m_folder(""), m_event_counter(m_numEC), m_EC_names(m_numEC), m_event_weight(1), m_passCutsTrigger(true), m_passCutThirdJetWeight(true), m_primaryVertex(0), m_idxNN(0), m_bookedHistos(0) { declareProperty( "ElecSelector", m_selectElec); declareProperty( "PhotonSelector", m_selectPhot); declareProperty( "MuonSelector", m_selectMuon); declareProperty( "JetSelector", m_selectJet); declareProperty( "JetLowPtSelector", m_selectJetlowpt); declareProperty( "JetForMETSelector", m_selectJetformet); declareProperty( "TauSelector", m_selectTau); declareProperty( "TruthSelector", m_selectTruth); declareProperty( "WeightSvc", m_osWeightSvc); declareProperty("TrackParticleContainer", m_trackParticleContainerName="TrackParticleCandidate"); declareProperty("VectorBoson", m_VectorBosonOpS = "eZ"); // declareProperty("VectorBoson2", m_VectorBoson2OpS = "eZbb"); // "JetFitterCOMBNN" declareProperty("JetTagger", m_btaggerName = ""); //Default combination declareProperty("TruthLater", m_truthLater= false); // Set to true if you want to unpack truth later after some rec cuts and so speed things up declareProperty("OSHistSvc", m_osHistSvc); declareProperty("beamCondSvc", m_beamCondSvc); declareProperty("Isolation", m_doIsol = true); declareProperty("DoPartonJets", m_doPartonJets = false); declareProperty("RemoveHadOverlapWithWZDecays", m_removeHadOverlapWZDecays = false); // Not needed if running the WZ truth filter declareProperty("LLnoMFake", m_llnomFake = false); declareProperty( "KinFitter", m_kinFitter); } OSHbb::~OSHbb() {} StatusCode OSHbb::initialize() { msg(MSG::INFO) << "Initializing OSHbb" << endreq; // base class OSHBase::initialize(); // Get tools and services if( FetchToolsAndServices().isFailure() ) { msg(MSG::ERROR) << "OSHbb::Initialize failed FetchToolsAndServices" << endreq; return StatusCode::FAILURE; } // OSHistSvc if (m_osHistSvc.retrieve().isFailure()) { msg(MSG::ERROR) << "OSHistSvc could NOT bo loaded." << endreq; return StatusCode::FAILURE; } // NN Ntuple m_osHistSvc->bookTree("/file1/NNVariables", "NN Ntuple"); m_osHistSvc->addBranch("/file1/NNVariables", "NN", &fNNVariables, "Jet0Pt:Jet0Eta:Jet0Phi:Jet0M:Jet0Weight:Jet1Pt:Jet1Eta:Jet1Phi:Jet1M:Jet1Weight:Jet2Pt:Jet2Eta:Jet2Phi:Jet2M:Jet2Weight:Lep0Pt:Lep0Eta:Lep0Phi:Lep0M:Lep1Pt:Lep1Eta:Lep1Phi:Lep1M:ZPt:ZM:DeltaEtaJJ:DeltaPhiJJ:DeltaEtaLL:DeltaPhiLL:DeltaPhiJJZ:Ht:EtMiss:NJets:MJJ:CosThetaStar:TopMass:SEtJets:SEtObjs:EvWeight"); fTHistStreamName="/file1/";// must have format /name/ and match the THistSvc.Output joboptions steering // Read the string for the type of selection and set accordingly. This could be an enum on the py side m_VectorBosonOp = eZ; //This is Z->e or mu if(0 == m_VectorBosonOpS.compare("eZnu")) m_VectorBosonOp = eZnu; else if(0 == m_VectorBosonOpS.compare("eZ")) m_VectorBosonOp = eZ; else if(0 == m_VectorBosonOpS.compare("eW")) m_VectorBosonOp = eW; // Type of second vector boson for H->ZZ analysis m_VectorBoson2Op = eZbb; //Warning order of comparison important here if (0==m_VectorBoson2OpS.compare("eZnu"))m_VectorBoson2Op = eZnu; // This is Z->nu else if(0==m_VectorBoson2OpS.compare("eZ"))m_VectorBoson2Op = eZ; // This is Z->e or mu std::string VectorBoson("No Boson Selected"); if (m_VectorBosonOp==eZ) VectorBoson=" Selecting Z->ee, Z->mumu (eZ)"; else if (m_VectorBosonOp==eW) VectorBoson=" Selecting W->enu, W->munu (eW)"; else if (m_VectorBosonOp==eZnu) VectorBoson=" Selecting Z->nunu (eZnu)"; msg(MSG::INFO) << VectorBoson << endreq; VectorBoson=("No Boson Selected"); if (m_VectorBoson2Op==eZ) VectorBoson="2nd Vector Boson Selecting Z->ee, Z->mumu (eZ)"; else if (m_VectorBoson2Op==eZbb) VectorBoson="2nd Vector Boson Selecting Z->bb (eZbb)"; else if (m_VectorBoson2Op==eZnu) VectorBoson="2nd Vector Boson Selecting Z->nunu (eZnu)"; msg(MSG::INFO) << VectorBoson << endreq; // Stream/Trigger selection - done in OSHBase // Different directories for histograms m_histDirCuts = fTHistStreamName+"Cuts/"; m_histDirTruth = fTHistStreamName+"Gen/"; m_histDirFake = fTHistStreamName+"Fake/"; std::cout << "Cuts Histo directory " << m_histDirCuts << std::endl; std::cout << "Truth Histo directory " << m_histDirTruth << std::endl; for (int i=0; i< m_numEC; i++) m_event_counter[i]=0; // Cuts m_wmasslow = 40*GeV; //ZMass ll used in note // m_zmasslow = 79*GeV; // m_zmasshigh = 103*GeV; //ZMass +-15 due to poorer resolution in data m_zmasslow = 76*GeV; m_zmasshigh = 106*GeV; m_zmasslowsb = 60*GeV; //We can't go lower cause of cut on PYTHIA m_zmasshighsb = 150*GeV; m_passCutZMass =false; //true if event has is within tight mass limits //ZMass bb m_zbmasslowsb = 40*GeV; m_zbmasshighsb = 150*GeV; m_zbmasslow = 70*GeV; // m_zbmasshigh = 110*GeV; m_zbmasshigh = 105*GeV; // MET m_etmisscut_eW = 30*GeV; m_etmisscut_eZ = 50*GeV; // m_etmisscut_eZ = 30*GeV; m_etmisscut_eZnu = (m_VectorBoson2Op==eZnu)? 100*GeV : 150*GeV; m_passCutEtmiss=true; //true if event passes etmiss cut // lepton m_ptcutelec_eW = 20*GeV; m_nelecmin_eW = 1; // minimum tight m_nelecmax_eW = 1; // max good m_etcutelec0_eZ = 20*GeV; // m_etcutelec1_eZ = 15*GeV; m_etcutelec1_eZ = 20*GeV; m_nelecmin_eZ = 1; // minimum tight, always at least 2 good m_nelecmax_eZ = 99; // max good m_nelecmax_eZnu = 0; // m_selectMuon->setProperty("MinPt", "10000"); m_ptcutmu_eW = 15*GeV; m_nmumin_eW = 1; // minimum tight m_nmumax_eW = 1; // max good m_ptcutmu0_eZ = 20*GeV; m_ptcutmu1_eZ = 20*GeV; m_nmumin_eZ = 1; // minimum tight, always at least 2 good m_nmumax_eZ = 99; // max good m_nmumax_eZnu = 0; // max good m_taucut = 0; // no taus alllowed m_deltaPhiJJCut = 1.570796; m_deltaPhiLLCut = 1.570796; m_PtZCut = 50*GeV; m_dibmassscale = 1.05; //Replaces scale in OSSelectJet //m_RwEtaBError = 0.7; //m_RwEtBError = 0.3; m_RwEtaBError = 0.6; m_RwEtBError = 0.5; m_RwEtBCen = 50.*GeV; // central Et reweight value in GeV m_RwEtaBCen =1; // central eta reweight value m_maxetabjet =2.5; //Maximum eta for b tagged jets if(m_btaggerName.find("JetFitter")!=std::string::npos) { if (m_VectorBosonOp==eZnu||m_VectorBosonOp==eZ) { m_jetWeight0Cut=0; m_jetWeight1Cut=0; m_jetWeight2Cut=0; //This is an anti-b cut } else { m_jetWeight0Cut=3; m_jetWeight1Cut=3; m_jetWeight2Cut=3; //This is an anti-b cut } } else if(m_btaggerName.find("SV0")!=std::string::npos) { m_jetWeight0Cut=5.85; m_jetWeight1Cut=5.85; m_jetWeight2Cut=5.85; //This is an anti-b cut } else { if (m_VectorBosonOp==eZnu||m_VectorBosonOp==eZ) { m_jetWeight0Cut=1; m_jetWeight1Cut=1; m_jetWeight2Cut=1; //This is an anti-b cut } else { if(m_VectorBosonOp==eW) { // looser cuts for WH m_jetWeight0Cut=0; m_jetWeight1Cut=0; } else { m_jetWeight0Cut=3; m_jetWeight1Cut=3; } m_jetWeight2Cut=0; //This is an anti-b cut } } m_passCutsTrigger=true; // TEST //m_jetWeight0Cut=-10000; //m_jetWeight1Cut=-10000; //m_jetWeight2Cut=-10000; m_sysNames[0]="SysJetEUp"; m_sysNames[1]="SysJetEDo"; m_sysNames[2]="SysJetEResol"; m_sysNames[3]="SysElecEUp"; m_sysNames[4]="SysElecEDo"; m_sysNames[5]="SysElecEResol"; m_sysNames[6]="SysMuonEUp"; m_sysNames[7]="SysMuonEDo"; m_sysNames[8]="SysMuonEResolIDUp"; m_sysNames[9]="SysMuonEResolIDDo"; m_sysNames[10]="SysMuonEResolMSUp"; m_sysNames[11]="SysMuonEResolMSDo"; m_sysNames[12]="SysMuonEResolNoSmear"; m_sysNames[13]="SysMETEUp"; m_sysNames[14]="SysMETEDo"; m_sysNames[15]="SysMETEResol"; // Event counter names m_EC_names[eECAll]= " All "; m_EC_names[eEC1GoodE]= ">=1 good electrons "; m_EC_names[eEC2GoodE]= ">=2 good electrons "; m_EC_names[eEC1TightE]= ">=1 medium electrons "; m_EC_names[eEC1GoodMu]= ">=1 good muons "; m_EC_names[eEC2GoodMu]= ">=2 good muons "; m_EC_names[eEC1TightMu]=">=1 medium muons "; m_EC_names[eECAllElectronCuts]=" all electron cuts "; m_EC_names[eECAllMuonCuts]=" all muon cuts "; m_EC_names[eECOnly1Lep]="exactly 1 lepton "; m_EC_names[eECNoTau] ="no taus "; m_EC_names[eECMET] ="MET cuts "; m_EC_names[eECLepMass]= "lep. mass requirements"; m_EC_names[eECNumJetsLT4]="<4 Jets "; m_EC_names[eECJetWeight2]="Jet 2 Weight Cut "; m_EC_names[eECNumJetsGT1]=">1 Jet "; m_EC_names[eECNumBJetsGT1]=">1 bJet "; m_EC_names[eECBJetMass] = stringify(m_zbmasslow)+" < M(bb) < "+ stringify(m_zbmasshigh); m_EC_names[eECDeltaPhiCut] = " delta phi cut.(|dphi|<"+stringify(m_deltaPhiJJCut)+")"; PrintCuts("OSHbb"); // Histograms BookHistos(); return StatusCode::SUCCESS; } void OSHbb::PrintCuts(std::string algname) { std::cout << "============================ " << std::endl; std::cout << " Cuts for " << algname << " " << this->name() << std::endl; std::cout << "============================ " << std::endl; if (m_VectorBosonOp==eZ) { std::cout << " eZ Cuts " << std::endl; std::cout << " Always require at least 2 good leptons " << std::endl; std::cout << " Et miss cut : " << m_etmisscut_eZ << std::endl; std::cout << " Min num (medium) elec: " << m_nelecmin_eZ << std::endl; std::cout << " Max num (good) elec: " << m_nelecmax_eZ << std::endl; std::cout << " Pt cut elec1: " << m_etcutelec0_eZ << std::endl; std::cout << " Pt cut elec2: " << m_etcutelec1_eZ << std::endl; std::cout << " Min num (medium) mu: " << m_nmumin_eZ << std::endl; std::cout << " Max num (good) mu: " << m_nmumax_eZ << std::endl; std::cout << " Pt cut mu1: " << m_ptcutmu0_eZ << std::endl; std::cout << " Pt cut mu2: " << m_ptcutmu1_eZ << std::endl; } else if (m_VectorBosonOp==eW) { std::cout << " eW Cuts " << std::endl; std::cout << " Et miss cut : " << m_etmisscut_eW << std::endl; std::cout << " Min num (medium) elec: " << m_nelecmin_eW << std::endl; std::cout << " Max num (good) elec: " << m_nelecmax_eW << std::endl; std::cout << " Pt cut elec1: " << m_ptcutelec_eW << std::endl; std::cout << " Min num (medium) mu: " << m_nmumin_eW << std::endl; std::cout << " Max num (good) mu: " << m_nmumax_eW << std::endl; std::cout << " Pt cut mu1: " << m_ptcutmu_eW << std::endl; } else if (m_VectorBosonOp==eZnu) { std::cout << " eZnu Cuts " << std::endl; std::cout << " Et miss cut : " << m_etmisscut_eZnu << std::endl; std::cout << " Max num elec: " << m_nelecmax_eZnu << std::endl; std::cout << " Max num mu: " << m_nmumax_eZnu << std::endl; } std::cout << " No Taus allowed " << m_taucut << std::endl; std::cout << " Jet Weight Cut 0: " << m_jetWeight0Cut << std::endl; std::cout << " Jet Weight Cut 1: " << m_jetWeight1Cut << std::endl; return; } void OSHbb::BookHistos() { BookGenHistos(); BookMETHistos(); BookZHistos(); BookJetHistos(); BookJetCorrectionHistos(); BookElectronHistos(); BookMuonHistos(); BookHighProbHistos(); BookJetTagInfoHistos(); BookCutFlowHistos(); BookZSBHistos(); BookIsolHistos(); //Default histos for final fits BookFinalHistos(""); //Systematics histos for final fits for (int i=0; i< nSys; i++){ BookFinalHistos(m_sysNames[i]); } m_bookedHistos = 1; return; } void OSHbb::BookMuonHistos() { // Muon histogram booking m_histDirMuon = fTHistStreamName+"Muon/"; m_osHistSvc->bookFlavHist(m_histDirMuon+"MuonPt","p_{t} #mu",50,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMuon+"MuonEta","#eta #mu",70,-3.5,3.5); m_osHistSvc->bookFlavHist(m_histDirMuon+"MuonChi2","chi2 mu",200,0.0,1000.0); m_osHistSvc->bookFlavHist(m_histDirMuon+"MuonQuality","Quality #mu", 4, 0., 4.); m_osHistSvc->bookFlavHist(m_histDirMuon+"MuonD0","D0 #mu", 1000, -10., 10.); m_osHistSvc->bookFlavHist(m_histDirMuon+"MuonDZ0ZVertex","DZ0ZVertex #mu", 310, -10., 300.); m_osHistSvc->bookFlavHist(m_histDirMuon+"MuonIPSig","IPSig #mu", 200, -20., 20.); m_osHistSvc->bookFlavHist(m_histDirMuon+"MuonD0MZ","D0 #mu", 1000, -10., 10.); m_osHistSvc->bookFlavHist(m_histDirMuon+"MuonDZ0ZVertexMZ","DZ0ZVertex #mu", 310, -10., 300.); m_osHistSvc->bookFlavHist(m_histDirMuon+"MuonIPSigMZ","IPSig #mu", 200, -20., 20.); m_osHistSvc->bookFlavHist(m_histDirMuon+"MuonPtAfterTrig","p_{t} #mu",50,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMuon+"MuonPtBothAfterTrig","p_{t} #mu",50,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMuon+"MuonEtaAfterTrig","#eta #mu",70,-3.5,3.5); m_osHistSvc->bookFlavHist(m_histDirMuon+"MuonEtaBothAfterTrig","#eta #mu",70,-3.5,3.5); return; } void OSHbb::BookElectronHistos() { // Electron histogram booking m_histDirElec = fTHistStreamName+"Electron/"; m_osHistSvc->bookFlavHist(m_histDirElec+"ElecPt","pt el",50,0.,250.); m_osHistSvc->bookFlavHist(m_histDirElec+"ElecEta","eta el",70,-3.5,3.5); m_osHistSvc->bookFlavHist(m_histDirElec+"ElecPtAfterTrig","pt el",50,0.,250.); m_osHistSvc->bookFlavHist(m_histDirElec+"ElecPtBothAfterTrig","pt el",50,0.,250.); m_osHistSvc->bookFlavHist(m_histDirElec+"ElecEtaAfterTrig","eta el",70,-3.5,3.5); m_osHistSvc->bookFlavHist(m_histDirElec+"ElecEtaBothAfterTrig","eta el",70,-3.5,3.5); m_osHistSvc->bookFlavHist(m_histDirElec+"ElecDZ0ZVertex","DZ0ZVertex e", 310, -10., 300.); m_osHistSvc->bookFlavHist(m_histDirElec+"ElecDZ0ZVertexMZ","DZ0ZVertex e", 310, -10., 300.); return; } void OSHbb::BookZHistos() { // jet BHist histograms m_histDirZ = fTHistStreamName+"Z/"; m_osHistSvc->bookFlavHist(m_histDirZ+"PtZ","pt Z",50,0.,400.); m_osHistSvc->bookFlavHist(m_histDirZ+"MZ","ZMass",150,0.,150.); m_osHistSvc->bookFlavHist(m_histDirZ+"MZee","ZMass ee ",150,0.,150.); m_osHistSvc->bookFlavHist(m_histDirZ+"MZmumu","ZMass mumu",150,0.,150.); m_osHistSvc->bookFlavHist(m_histDirZ+"MZeeOrig","ZMass ee orig",150,0.,150.); m_osHistSvc->bookFlavHist(m_histDirZ+"MZmumuOrig","ZMass mumu orig",150,0.,150.); m_osHistSvc->bookFlavHist(m_histDirZ+"MZAllMass","ZMass",50,0.,150.); m_osHistSvc->bookFlavHist(m_histDirZ+"MZeeAllMass","ZMass ee ",50,0.,150.); m_osHistSvc->bookFlavHist(m_histDirZ+"MZmumuAllMass","ZMass mumu",50,0.,150.); m_osHistSvc->bookFlavHist(m_histDirZ+"MZee_e10_loose","ZMass ee ",150,0.,150.); m_osHistSvc->bookFlavHist(m_histDirZ+"MZee_e10_medium","ZMass ee ",150,0.,150.); m_osHistSvc->bookFlavHist(m_histDirZ+"MZee_e15_loose","ZMass ee ",150,0.,150.); m_osHistSvc->bookFlavHist(m_histDirZ+"MZmumu_mu10","ZMass mumu",150,0.,150.); m_osHistSvc->bookFlavHist(m_histDirZ+"MZmumu_mu13","ZMass mumu",150,0.,150.); m_osHistSvc->bookFlavHist(m_histDirZ+"MZmumu_mu15","ZMass mumu",150,0.,150.); m_osHistSvc->bookFlavHist(m_histDirZ+"MZmumu_mu20","ZMass mumu",150,0.,150.); m_osHistSvc->bookFlavHist(m_histDirZ+"MZeeSameQ","ZMass ee SameQ",150,0.,150.); m_osHistSvc->bookFlavHist(m_histDirZ+"MZmumuSameQ","ZMass mumu SameQ",150,0.,150.); m_osHistSvc->bookFlavHist(m_histDirZ+"MZeePosCrack","ZMass ee Crack",150,0.,150.); m_osHistSvc->bookFlavHist(m_histDirZ+"MZeeNegCrack","ZMass ee Crack",150,0.,150.); m_osHistSvc->bookFlavHist(m_histDirZ+"MZeeNoCrack","ZMass ee No Crack",150,0.,150.); m_osHistSvc->bookFlavHist(m_histDirZ+"MZeeHighMET","ZMass ee high MET",150,0.,150.); m_osHistSvc->bookFlavHist(m_histDirZ+"MZmumuHighMET","ZMass mumu high MET",150,0.,150.); m_osHistSvc->bookFlavHist(m_histDirZ+"PtZDijet","pt Z (2 jets)",50,0.,400.); m_osHistSvc->bookFlavHist(m_histDirZ+"MZDijet","M Z (2 jets)",150,0.,150.); m_osHistSvc->bookFlavHist(m_histDirZ+"MZeeDijet","ZMass ee (2 jets)",150,0.,150.); m_osHistSvc->bookFlavHist(m_histDirZ+"MZmumuDijet","ZMass mumu (2 jets)",150,0.,150.); m_osHistSvc->bookFlavHist(m_histDirZ+"MZeeDijetWC","ZMass ee (2 jets)",150,0.,150.); m_osHistSvc->bookFlavHist(m_histDirZ+"MZmumuDijetWC","ZMass mumu (2 jets)",150,0.,150.); /* m_osHistSvc->bookFlavHist(m_histDirZ+"MZFake","ZMassFake",50,0.,150.); m_osHistSvc->bookFlavHist(m_histDirZ+"MZDijetFake","ZMassFake (2 jets)",50,0.,150.); m_osHistSvc->bookFlavHist(m_histDirZ+"MZDijetWCFake","ZMassFake (2 tagged jets)",50,0.,150.); m_osHistSvc->bookFlavHist(m_histDirZ+"MZeeFake","ZMassFake",50,0.,150.); m_osHistSvc->bookFlavHist(m_histDirZ+"MZeeDijetFake","ZMassFake (2 jets)",50,0.,150.); m_osHistSvc->bookFlavHist(m_histDirZ+"MZeeDijetWCFake","ZMassFake (2 tagged jets)",50,0.,150.); */ // trigger m_osHistSvc->bookFlavHist(m_histDirZ+"PtZDijetNoTrig","pt Z (2 jets)",50,0.,400.); m_osHistSvc->bookFlavHist(m_histDirZ+"PtZDijetTrig","pt Z (2 jets)",50,0.,400.); m_osHistSvc->bookFlavHist(m_histDirZ+"PtZDijetWC","pt Z (2 tagged jets)",50,0.,400.); m_osHistSvc->bookFlavHist(m_histDirZ+"MZDijetWC","M Z (2 tagged jets)",150,0.,150.); m_osHistSvc->bookFlavHist(m_histDirZ+"MZDijet1WC","M Z (2 tagged jets)",150,0.,150.); m_osHistSvc->bookFlavHist(m_histDirZ+"AngleEE","",180,0.0,180.0); m_osHistSvc->bookFlavHist(m_histDirZ+"AngleMuMu","",180,0.0,180.0); m_osHistSvc->bookFlavHist(m_histDirZ+"AngleEEMZ","",180,0.0,180.0); m_osHistSvc->bookFlavHist(m_histDirZ+"AngleMuMuMZ","",180,0.0,180.0); return; } void OSHbb::BookJetHistos() { // jet BHist histograms m_histDirJets = fTHistStreamName+"Jet/"; m_histDirNN = fTHistStreamName+"NN/"; // number of jets m_osHistSvc->bookFlavHist(m_histDirJets+"NJet","N^{jet}",10,0.,10); m_osHistSvc->bookFlavHist(m_histDirJets+"ElecNJet","N^{jet}",10,0.,10); m_osHistSvc->bookFlavHist(m_histDirJets+"MuonNJet","N^{jet}",10,0.,10); m_osHistSvc->bookFlavHist(m_histDirJets+"NBJet","N^{jet}",10,0.,10); m_osHistSvc->bookFlavHist(m_histDirJets+"ElecNBJet","N^{jet}",10,0.,10); m_osHistSvc->bookFlavHist(m_histDirJets+"MuonNBJet","N^{jet}",10,0.,10); m_osHistSvc->bookFlavHist(m_histDirJets+"NBJetDijet","N^{jet}",10,0.,10); m_osHistSvc->bookFlavHist(m_histDirJets+"ElecNBJetDijet","N^{jet}",10,0.,10); m_osHistSvc->bookFlavHist(m_histDirJets+"NBJetNJet1","N^{jet}",10,0.,10); m_osHistSvc->bookFlavHist(m_histDirJets+"ElecNBJetNJet1","N^{jet}",10,0.,10); m_osHistSvc->bookFlavHist(m_histDirJets+"MuonNBJetNJet1","N^{jet}",10,0.,10); m_osHistSvc->bookFlavHist(m_histDirJets+"MuonNBJetDijet","N^{jet}",10,0.,10); m_osHistSvc->bookFlavHist(m_histDirJets+"NBJetNJet2","N^{jet}",10,0.,10); m_osHistSvc->bookFlavHist(m_histDirJets+"ElecNBJetNJet2","N^{jet}",10,0.,10); m_osHistSvc->bookFlavHist(m_histDirJets+"MuonNBJetNJet2","N^{jet}",10,0.,10); m_osHistSvc->bookFlavHist(m_histDirJets+"NBJetNJet3","N^{jet}",10,0.,10); m_osHistSvc->bookFlavHist(m_histDirJets+"ElecNBJetNJet3","N^{jet}",10,0.,10); m_osHistSvc->bookFlavHist(m_histDirJets+"MuonNBJetNJet3","N^{jet}",10,0.,10); m_osHistSvc->bookFlavHist(m_histDirJets+"NJetBTag","N^{jet}",10,0.,10); // Ht m_osHistSvc->bookFlavHist(m_histDirJets+"Ht", " ",100,0.0,1000); m_osHistSvc->bookFlavHist(m_histDirJets+"HtDijet", " ",100,0.0,1000); m_osHistSvc->bookFlavHist(m_histDirJets+"HtDijetNBJet1", " ",100,0.0,1000); m_osHistSvc->bookFlavHist(m_histDirJets+"HtDijetNBJet2", " ",100,0.0,1000); m_osHistSvc->bookFlavHist(m_histDirJets+"HtNJet2NBJet1", " ",100,0.0,1000); m_osHistSvc->bookFlavHist(m_histDirJets+"HtNJet2NBJet2", " ",100,0.0,1000); m_osHistSvc->bookFlavHist(m_histDirJets+"HtNJet3NBJet1", " ",100,0.0,1000); m_osHistSvc->bookFlavHist(m_histDirJets+"HtNJet3NBJet2", " ",100,0.0,1000); // Et m_osHistSvc->bookFlavHist(m_histDirJets+"JetEt","E_{T}^{jet}",50,0.,250.); m_osHistSvc->bookFlavHist(m_histDirJets+"JetEtSysEUp","E_{T}^{jet}",50,0.,250.); m_osHistSvc->bookFlavHist(m_histDirJets+"JetEtSysEDo","E_{T}^{jet}",50,0.,250.); /* //For fakes m_osHistSvc->bookFlavHist(m_histDirFake+"PtJet","P_{T}^{jet}",50,0.,250.); m_osHistSvc->bookFlavHist(m_histDirFake+"PtJetPassE","P_{T}^{jet}",50,0.,250.); m_osHistSvc->bookFlavHist(m_histDirFake+"PtEPtJetRat","",100,0.0,2.0); m_osHistSvc->bookFlavHist(m_histDirFake+"PtJetPassMu","P_{T}^{jet}",50,0.,250.); m_osHistSvc->bookFlavHist(m_histDirFake+"PtMuPtJetRat","",100,0.0,2.0); // trigger efficiency m_osHistSvc->bookFlavHist(m_histDirFake+"PtJetTrig","P_{T}^{jet}",50,0.,250.); m_osHistSvc->bookFlavHist(m_histDirFake+"PtJetTrigMon","P_{T}^{jet}",50,0.,250.); m_osHistSvc->bookFlavHist(m_histDirFake+"PtJetTrigActual","P_{T}^{jet}",50,0.,250.); m_osHistSvc->bookFlavHist(m_histDirFake+"PtJetTrigRaw","P_{T}^{jet}",50,0.,250.); */ // model error m_osHistSvc->bookFlavHist(m_histDirJets+"JetEtSysRwEtUp","E_{T}^{jet}",50,0.,250.); m_osHistSvc->bookFlavHist(m_histDirJets+"JetEtSysRwEtDo","E_{T}^{jet}",50,0.,250.); // eta m_osHistSvc->bookFlavHist(m_histDirJets+"JetEta","{#eta}^{jet}",70,-3.5,3.5); m_osHistSvc->bookFlavHist(m_histDirJets+"JetEtaSysRwEtaUp","{#eta}^{jet}",70,-3.5,3.5); m_osHistSvc->bookFlavHist(m_histDirJets+"JetEtaSysRwEtaDo","{#eta}^{jet}",70,-3.5,3.5); // mass m_osHistSvc->bookFlavHist(m_histDirJets+"DijetMassPtZ50","M_{j#bar{j}} (P_{T}^{W}>50 GeV)",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetMassPtZ100","M_{j#bar{j}} (P_{T}^{W}>100 GeV)",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetMassNVtx1","M_{j#bar{j}}",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetMassNVtxGt1","M_{j#bar{j}}",50,0.,250); //Masses for lost HF events in the jet selection m_osHistSvc->bookFlavHist(m_histDirJets+"Lost02DijetMass02","M_{j#bar{j}}",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+"Lost0Dijet","M_{j#bar{j}}",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+"Lost12DijetMass12","M_{j#bar{j}}",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+"Lost1DijetMass","M_{j#bar{j}}",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+"Lost2DijetMass","M_{j#bar{j}}",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+"LostDijetMass","M_{j#bar{j}}",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetMassWCNoCal","M_{b#bar{b}}",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetMassWCPtZ50","M_{b#bar{b}}",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetMassWCPtZ100","M_{b#bar{b}}",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassPtZ","M_{H}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassPtZWC","M_{H}",400,5,1005); //Higgs Mass for only LL LB or LC combinations m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassWCL","M_{H}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassDPhiCutsWCL","M_{H}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetMassWCL","M_{b#bar{b}}",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassNNCutLoose","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassNNCut","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassNNCutTight","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassNNloMHCutLoose","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassNNloMHCutLooseWC","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassNNloMHCut","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassNNloMHCutWC","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassNNloMHCutTight","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassNNloMHCutTightWC","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassNNhiMHCutLoose","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassNNhiMHCutLooseWC","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassNNhiMHCut","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassNNhiMHCutWC","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassNNhiMHCutTight","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassNNhiMHCutTightWC","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassSBLow","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassSBLowPtJ50","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassSBLowDPhiCuts","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassSBLowPtJ50DPhiCuts","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassSBLowDPhiCutsFake","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassSBLowWC","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassSBLowDPhiCutsWC","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassSBLowNNCut","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassSBLowNNCutWC","M_{j#bar{j}}",400,5,1005); /* m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassSBLowFake","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassSBLowWCFake","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassSBLowDPhiCutsWCFake","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassSBHighFake","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassSBHighDPhiCutsFake","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassSBHighWCFake","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassSBHighDPhiCutsWCFake","M_{j#bar{j}}",400,5,1005); */ m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassSBHigh","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassSBHighPtJ50","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassSBHighDPhiCuts","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassSBHighPtJ50DPhiCuts","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassSBHighWC","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassSBHighDPhiCutsWC","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassSBHighNNCut","M_{j#bar{j}}",400,5,1005); m_osHistSvc->bookFlavHist(m_histDirJets+"DijetZMassSBHighNNCutWC","M_{j#bar{j}}",400,5,1005); //NN inputs m_osHistSvc->bookFlavHist(m_histDirNN+"NNptZ"," ",100,0.0,1000.0); m_osHistSvc->bookFlavHist(m_histDirNN+"NNjetEt0Sc"," ",100,0.0,4); m_osHistSvc->bookFlavHist(m_histDirNN+"NNjetEt1Sc"," ",100,0.0,4); m_osHistSvc->bookFlavHist(m_histDirNN+"NNjetEt0"," ",100,0.0,1000); m_osHistSvc->bookFlavHist(m_histDirNN+"NNjetEt1"," ",100,0.0,1000); m_osHistSvc->bookFlavHist(m_histDirNN+"NNjjDeltaR"," ",100,0.0,4); m_osHistSvc->bookFlavHist(m_histDirNN+"NNjjZDeltaPhi"," ",100,0.0,4); m_osHistSvc->bookFlavHist(m_histDirNN+"NNHtSc", " ",100,0.0,10); m_osHistSvc->bookFlavHist(m_histDirNN+"NNHt", " ",100,0.0,1000); m_osHistSvc->bookFlavHist(m_histDirNN+"NNmjj"," ",100,0.0,400); m_osHistSvc->bookFlavHist(m_histDirNN+"NNEtmiss"," ",100,0.0,400); m_osHistSvc->bookFlavHist(m_histDirNN+"NNnumjets"," ",10,-0.5,9.5); m_osHistSvc->bookFlavHist(m_histDirNN+"NNmZ"," ",100,0.0,200); m_osHistSvc->bookFlavHist(m_histDirNN+"NNtrPt0"," ",100,0.0,200); m_osHistSvc->bookFlavHist(m_histDirNN+"NNtrPt1"," ",100,0.0,200); m_osHistSvc->bookFlavHist(m_histDirNN+"NNjetEt2"," ",100,0.0,1000); m_osHistSvc->bookFlavHist(m_histDirNN+"NNSEtjets"," ",100,0.0,2000); m_osHistSvc->bookFlavHist(m_histDirNN+"NNSEtnobjets"," ",100,0.0,2000); m_osHistSvc->bookFlavHist(m_histDirNN+"NNtopmassAll"," ",100,0.0,1000); m_osHistSvc->bookFlavHist(m_histDirNN+"NNtopmassNear"," ",100,0.0,1000); m_osHistSvc->bookFlavHist(m_histDirNN+"NNptZWC"," ",100,0.0,1000.0); m_osHistSvc->bookFlavHist(m_histDirNN+"NNjetEt0ScWC"," ",100,0.0,4); m_osHistSvc->bookFlavHist(m_histDirNN+"NNjetEt1ScWC"," ",100,0.0,4); m_osHistSvc->bookFlavHist(m_histDirNN+"NNjetEt0WC"," ",100,0.0,1000); m_osHistSvc->bookFlavHist(m_histDirNN+"NNjetEt1WC"," ",100,0.0,1000); m_osHistSvc->bookFlavHist(m_histDirNN+"NNjjDeltaRWC"," ",100,0.0,4); m_osHistSvc->bookFlavHist(m_histDirNN+"NNjjZDeltaPhiWC"," ",100,0.0,4); m_osHistSvc->bookFlavHist(m_histDirNN+"NNHtScWC", " ",100,0.0,10); m_osHistSvc->bookFlavHist(m_histDirNN+"NNHtWC", " ",100,0.0,1000); m_osHistSvc->bookFlavHist(m_histDirNN+"NNmjjWC"," ",100,0.0,400); m_osHistSvc->bookFlavHist(m_histDirNN+"NNEtmissWC"," ",100,0.0,400); m_osHistSvc->bookFlavHist(m_histDirNN+"NNnumjetsWC"," ",10,-0.5,9.5); m_osHistSvc->bookFlavHist(m_histDirNN+"NNmZWC"," ",100,0.0,10000); m_osHistSvc->bookFlavHist(m_histDirNN+"NNtrPt0WC"," ",100,0.0,200); m_osHistSvc->bookFlavHist(m_histDirNN+"NNtrPt1WC"," ",100,0.0,200); m_osHistSvc->bookFlavHist(m_histDirNN+"NNjetEt2WC"," ",100,0.0,1000); m_osHistSvc->bookFlavHist(m_histDirNN+"NNSEtjetsWC"," ",100,0.0,2000); m_osHistSvc->bookFlavHist(m_histDirNN+"NNSEtnobjetsWC"," ",100,0.0,2000); m_osHistSvc->bookFlavHist(m_histDirNN+"NNtopmassAllWC"," ",100,0.0,1000); m_osHistSvc->bookFlavHist(m_histDirNN+"NNtopmassNearWC"," ",100,0.0,1000); //Make a grid of NNcuts the first index is the cut on Top the second on W+Jets for (int i=0; i<10; i++){ for (int j=0; j<10; j++){ TString name="DijetMassNNCutG"; if(i<10) name+="0"; name+=i; if(j<10) name+="0"; name+=j; TH1F* hist = new TH1F(m_histDirNN+name,"M_{j#bar{j}} (P_{T}^{W}>100 GeV) NNCut",50,0.,250); m_osHistSvc->bookFlavHist(hist); TString name2="DijetMassNNCutWCG"; if(i<10) name2+="0"; name2+=i; if(j<10) name2+="0"; name2+=j; TH1F* hist2 = new TH1F(m_histDirNN+name2,"M_{j#bar{j}} (P_{T}^{W}>100 GeV) NNCut",50,0.,250); m_osHistSvc->bookFlavHist(hist2); } } float weightmin=-10; float weightmax=25; int nweightbins=35; m_osHistSvc->bookFlavHist(m_histDirJets+"JetWeight","jet weight",nweightbins,weightmin,weightmax); m_osHistSvc->bookFlavHist(m_histDirJets+"JetWeightSysEUp","weight b truth",nweightbins,weightmin,weightmax); m_osHistSvc->bookFlavHist(m_histDirJets+"JetWeightSysEDo","weight b truth",nweightbins,weightmin,weightmax); m_osHistSvc->bookFlavHist(m_histDirJets+"JetWeightSysRwEtUp","weight b truth",nweightbins,weightmin,weightmax); m_osHistSvc->bookFlavHist(m_histDirJets+"JetWeightSysRwEtDo","weight b truth",nweightbins,weightmin,weightmax); m_osHistSvc->bookFlavHist(m_histDirJets+"JetWeightSysRwEtaUp","weight b truth",nweightbins,weightmin,weightmax); m_osHistSvc->bookFlavHist(m_histDirJets+"JetWeightSysRwEtaDo","weight b truth",nweightbins,weightmin,weightmax); m_osHistSvc->bookFlavHist(m_histDirJets+"JetWeightIP3D","IP3D weight",nweightbins,weightmin,weightmax); m_osHistSvc->bookFlavHist(m_histDirJets+"JetWeightIP2D","IP2D weight",nweightbins,weightmin,weightmax); m_osHistSvc->bookFlavHist(m_histDirJets+"JetWeightSV1","SV1 weight",nweightbins,weightmin,weightmax); m_osHistSvc->bookFlavHist(m_histDirJets+"JetWeightSV0","SV0 weight",30,-10.0,50.0); // for njet=1,2 for W+b jets m_osHistSvc->bookFlavHist(m_histDirJets+"JetWeightNJet1","jet weight",nweightbins,weightmin,weightmax); m_osHistSvc->bookFlavHist(m_histDirJets+"JetWeightIP3DNJet1","IP3D weight",nweightbins,weightmin,weightmax); m_osHistSvc->bookFlavHist(m_histDirJets+"JetWeightSV1NJet1","SV1 weight",nweightbins,weightmin,weightmax); m_osHistSvc->bookFlavHist(m_histDirJets+"JetWeightSV0NJet1","SV0 weight",30,-10.0,50.0); m_osHistSvc->bookFlavHist(m_histDirJets+"JetWeightNJet2","jet weight",nweightbins,weightmin,weightmax); m_osHistSvc->bookFlavHist(m_histDirJets+"JetWeightIP3DNJet2","IP3D weight",nweightbins,weightmin,weightmax); m_osHistSvc->bookFlavHist(m_histDirJets+"JetWeightSV1NJet2","SV1 weight",nweightbins,weightmin,weightmax); m_osHistSvc->bookFlavHist(m_histDirJets+"JetWeightSV0NJet2","SV0 weight",30,-10.0,50.0); // m_osHistSvc->bookFlavHist(m_histDirJets+"DPhiJJ","#Delta #phi JJ",100,0.0,TMath::Pi()); m_osHistSvc->bookFlavHist(m_histDirJets+"DPhiJJBeforeMCut","#Delta #phi JJ",100,0.0,TMath::Pi()); m_osHistSvc->bookFlavHist(m_histDirJets+"DPhiJJWC","#Delta #phi JJ WC",100,0.0,TMath::Pi()); m_osHistSvc->bookFlavHist(m_histDirZ+"DPhiLL","#Delta #phi ll",100,0.0,TMath::Pi()); m_osHistSvc->bookFlavHist(m_histDirZ+"DPhiLLBeforeMCut","#Delta #phi ll",100,0.0,TMath::Pi()); m_osHistSvc->bookFlavHist(m_histDirZ+"DPhiLLWC","#Delta #phi ll WC",100,0.0,TMath::Pi()); m_osHistSvc->bookFlavHist(m_histDirZ+"DPhiLLDPhiJJCut","#Delta #phi ll",100,0.0,TMath::Pi()); m_osHistSvc->bookFlavHist(m_histDirZ+"DPhiLLDPhiJJCutWC","#Delta #phi ll WC",100,0.0,TMath::Pi()); m_osHistSvc->bookFlavHist(m_histDirJets+"JetFitterWeight","jet weight",nweightbins,weightmin,weightmax); m_osHistSvc->bookFlavHist(m_histDirJets+"JetFitterCombWeight","jet weight",nweightbins,weightmin,weightmax); m_osHistSvc->bookFlavHist(m_histDirJets+"JetFitterTagNNWeight","jet weight",nweightbins,weightmin,weightmax); m_osHistSvc->bookFlavHist(m_histDirJets+"JetFitterCOMBNNWeight","jet weight",nweightbins,weightmin,weightmax); m_osHistSvc->bookFlavHist(m_histDirJets+"JetFitterWeightNJet1","jet weight",nweightbins,weightmin,weightmax); m_osHistSvc->bookFlavHist(m_histDirJets+"JetFitterWeightNJet2","jet weight",nweightbins,weightmin,weightmax); //Jet Vertex Fraction m_osHistSvc->bookFlavHist(m_histDirJets+"JetVertexFraction","JVF",200,-2.0,2.0); m_osHistSvc->bookFlavHist(m_histDirJets+"JetVertexFractionCentral","JVF",200,-2.0,2.0); m_osHistSvc->bookFlavHist(m_histDirJets+"JetVertexFractionForward","JVF",200,-2.0,2.0); m_osHistSvc->bookFlavHist(m_histDirJets+"JetVertexFractionBarrel","JVF",200,-2.0,2.0); m_osHistSvc->bookFlavHist(m_histDirJets+"JetVertexFractionDijet","JVF",200,-2.0,2.0); m_osHistSvc->bookFlavHist(m_histDirJets+"JetVertexFractionDijetCentral","JVF",200,-2.0,2.0); m_osHistSvc->bookFlavHist(m_histDirJets+"JetVertexFractionDijetForward","JVF",200,-2.0,2.0); m_osHistSvc->bookFlavHist(m_histDirJets+"JetVertexFractionDijetBarrel","JVF",200,-2.0,2.0); m_osHistSvc->bookFlavHist(m_histDirJets+"JetVertexFractionEarly","JVF",200,-2.0,2.0); m_osHistSvc->bookFlavHist(m_histDirJets+"JetVertexFractionEarlyBarrel","JVF",200,-2.0,2.0); m_osHistSvc->bookFlavHist(m_histDirJets+"JetVertexFractionEarlyCentral","JVF",200,-2.0,2.0); m_osHistSvc->bookFlavHist(m_histDirJets+"JetVertexFractionEarlyForward","JVF",200,-2.0,2.0); return; } void OSHbb::BookFinalHistos(std::string tag) { //Book final histograms with and electron, muon version std::cout << "Booking final histos with tag "<< tag <bookFlavHist(m_histDirZ+pre+"ZMass"+tag+post,"MZ",250,0.,250); m_osHistSvc->bookFlavHist(m_histDirMET+pre+"MissingET"+tag+post,"MET",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetMass"+tag+post,"M_{j#bar{j}}",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetMassWC"+tag+post,"M_{b#bar{b}}",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetMass0WC"+tag+post,"M_{j#bar{j}}",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetMass1WC"+tag+post,"M_{j#bar{j}}",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetMassCutsLoose"+tag+post,"M_{b#bar{b}}",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetMassCutsMedium"+tag+post,"M_{b#bar{b}}",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetMassCutsTight"+tag+post,"M_{b#bar{b}}",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetMassCutsLooseWC"+tag+post,"M_{b#bar{b}}",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetMassCutsMediumWC"+tag+post,"M_{b#bar{b}}",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetMassCutsTightWC"+tag+post,"M_{b#bar{b}}",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetMassNNCutWJ"+tag+post,"M_{j#bar{j}} (P_{T}^{W}>100 GeV) NNCut W+Jets",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetMassNNCutWJWC"+tag+post,"M_{b#bar{b}} (P_{T}^{W}>100 GeV) BTag NNCut W+Jets",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetMassNNCutTT"+tag+post,"M_{j#bar{j}} (P_{T}^{W}>100 GeV) NNCut Ttbar",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetMassNNCutTTWC"+tag+post,"M_{b#bar{b}} (P_{T}^{W}>100 GeV) BTag NNCut Ttbar",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetMassNNCutLooseTT"+tag+post,"M_{j#bar{j}} (P_{T}^{W}>100 GeV) NNCut Ttbar",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetMassNNCutLooseTTWC"+tag+post,"M_{b#bar{b}} (P_{T}^{W}>100 GeV) BTag NNCut Ttbar",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetMassNNCutTightTT"+tag+post,"M_{j#bar{j}} (P_{T}^{W}>100 GeV) NNCut Ttbar",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetMassNNCutTightTTWC"+tag+post,"M_{b#bar{b}} (P_{T}^{W}>100 GeV) BTag NNCut Ttbar",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetMassNNCut"+tag+post,"M_{j#bar{j}} (P_{T}^{W}>100 GeV) NNCut Both",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetMassNNCutWC"+tag+post,"M_{b#bar{b}} (P_{T}^{W}>100 GeV) BTag NNCut Both",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetMassNNCutLoose"+tag+post,"M_{j#bar{j}} (P_{T}^{W}>100 GeV) NNCut Both",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetMassNNCutLooseWC"+tag+post,"M_{b#bar{b}} (P_{T}^{W}>100 GeV) BTag NNCut Both",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetMassNNCutTight"+tag+post,"M_{j#bar{j}} (P_{T}^{W}>100 GeV) NNCut Both",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetMassNNCutTightWC"+tag+post,"M_{b#bar{b}} (P_{T}^{W}>100 GeV) BTag NNCut Both",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetZMassDPhiCuts"+tag+post,"M_H",400,5.0,1005); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetZMassDPhiCutsWC"+tag+post,"M_H",400,5.0,1005); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetZMass"+tag+post,"M_H",400,5.0,1005); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetZMassWC"+tag+post,"M_H",400,5.0,1005); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetZMass0WC"+tag+post,"M_H",400,5.0,1005); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetZMass1WC"+tag+post,"M_H",400,5.0,1005); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetZMassNNCut0WC"+tag+post,"M_H",400,5.0,1005); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetZMassNNCut1WC"+tag+post,"M_H",400,5.0,1005); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetZMassNNCutWC"+tag+post,"M_H",400,5.0,1005); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetZMassNNCutTight0WC"+tag+post,"M_H",400,5.0,1005); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetZMassNNCutTight1WC"+tag+post,"M_H",400,5.0,1005); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetZMassNNCutTightWC"+tag+post,"M_H",400,5.0,1005); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetZMassNNCutLoose0WC"+tag+post,"M_H",400,5.0,1005); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetZMassNNCutLoose1WC"+tag+post,"M_H",400,5.0,1005); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetZMassNNCutLooseWC"+tag+post,"M_H",400,5.0,1005); //histograms with >=4 jets m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetZMassNJ4"+tag+post,"M_H",400,5.0,1005); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetZMassNJ4WC"+tag+post,"M_H",400,5.0,1005); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetZMassDPhiJJNJ4"+tag+post,"M_H",400,5.0,1005); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetZMassDPhiJJNJ4WC"+tag+post,"M_H",400,5.0,1005); //histograms with >=2 jets with pt> 50 GeV m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetZMassPtJ50"+tag+post,"M_H",400,5.0,1005); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetZMassPtJ50WC"+tag+post,"M_H",400,5.0,1005); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetZMassDPhiJJPtJ50"+tag+post,"M_H",400,5.0,1005); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetZMassDPhiJJPtJ50WC"+tag+post,"M_H",400,5.0,1005); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetZMassPtJ500WC"+tag+post,"M_H",400,5.0,1005); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetZMassDPhiJJPtJ500WC"+tag+post,"M_H",400,5.0,1005); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetZMassPtJ501WC"+tag+post,"M_H",400,5.0,1005); m_osHistSvc->bookFlavHist(m_histDirJets+pre+"DijetZMassDPhiJJPtJ501WC"+tag+post,"M_H",400,5.0,1005); } } void OSHbb::BookJetCorrectionHistos() { m_histDirJetCorr = fTHistStreamName+"JetCorrection/"; m_osHistSvc->bookFlavHist(m_histDirJetCorr+"JetEtMinusHad"," ",50,-50.,50.); m_osHistSvc->bookFlavHist(m_histDirJetCorr+"JetEt_vs_Had"," ",50,0.,250.,50,0.,250.); m_osHistSvc->bookFlavHist(m_histDirJetCorr+"JetEtMinusHad_vs_Had"," ",50,0.,250.,50,-50.,50.); return; } void OSHbb::BookMETHistos() { m_histDirMET = fTHistStreamName+"MissingET/"; m_osHistSvc->bookFlavHist(m_histDirMET+"METOrig","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METEle","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METGamma","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METTau","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METJet","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METRefMuon","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METCellOut","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METCryo","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METMuonBoy","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METTruthNonInt","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METTruthIntCent","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METTruthMuon","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METTruthSum","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"MET","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METNonSemiLep","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METSemiLep","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METElec","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METMuon","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METMZ","Missing E_{T} after Z Mass Cut",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METMZee","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METMZmumu","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METMZmumuCombined","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METMZeeExtraMu","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METMZmumuExtraMu","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METMZmumuSysMuonEResol","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"HT", "Scalar Sum Missing E_{T}",125,0.,1000.); m_osHistSvc->bookFlavHist(m_histDirMET+"MET_vs_HT","MET vs H_{T}",50,0.,1000.,50,0.,250.); // for QCD fakes m_osHistSvc->bookFlavHist(m_histDirMET+"MET Fake J0 Elec","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"MET Fake J1 Elec","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"MET Fake J2 Elec","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"MET Fake J3Plus Elec","Missing E_{T}",125,0.,250.); // and per b-jet for different jet multiplicities m_osHistSvc->bookFlavHist(m_histDirMET+"METNBJet0","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METNBJet1","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METNBJet2","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METNBJet3Plus","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METNJet1NBJet0","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METNJet1NBJet1","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METDijetNBJet0","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METDijetNBJet1","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METDijetNBJet2","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METDijetNBJet3Plus","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METNJet2NBJet0","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METNJet2NBJet1","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METNJet2NBJet2","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METNJet3NBJet0","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METNJet3NBJet1","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METNJet3NBJet2","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METNJet3NBJet3","Missing E_{T}",125,0.,250.); // after cuts against top m_osHistSvc->bookFlavHist(m_histDirMET+"METJetCut","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METElecJetCut","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METMuonJetCut","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METDijetCut","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METElecDijetCut","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METMuonDijetCut","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"MET1WC","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METElec1WC","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METMuon1WC","Missing E_{T}",125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METWC","Missing E_{T}",50,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METElecWC","Missing E_{T}",50,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"METMuonWC","Missing E_{T}",50,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"HTJetCut","Scalar Sum Missing E_{T}",125,0.,1000.); m_osHistSvc->bookFlavHist(m_histDirMET+"MET_vs_HTJetCut","MET vs H_{T}",50,0.,1000.,125,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"JetMETPhi","#Delta #phi (jet - MET)",45,0.,180); m_osHistSvc->bookFlavHist(m_histDirMET+"JetMETPhi2D","#Delta #phi (jet - MET) vs MET",45,0.,180,50,0.,250.); m_osHistSvc->bookFlavHist(m_histDirMET+"JetMETPhi2DClosest","#Delta #phi (jet - MET) vs MET",45,0.,180,50,0.,250.); return; } void OSHbb::BookZSBHistos() { m_histDirZSB=fTHistStreamName+"ZSB/"; m_osHistSvc->bookFlavHist(m_histDirZSB+"METZSB","Missing E_{T} in Z Mass side bands",50,0.,250.); m_osHistSvc->bookFlavHist(m_histDirZSB+"METZSBWC","Missing E_{T} in Z Mass side bands",50,0.,250.); m_osHistSvc->bookFlavHist(m_histDirZSB+"DijetMassZSBLowMET","M_{b#bar{b}} ",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirZSB+"DijetMassZSBHighMET","M_{b#bar{b}} ",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirZSB+"DijetMassZSBLowMETWC","M_{b#bar{b}} ",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirZSB+"DijetMassZSBHighMETWC","M_{b#bar{b}} ",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirZSB+"DijetZMassZSBLowMET","M_{ZZ}",400,5.,1005); m_osHistSvc->bookFlavHist(m_histDirZSB+"DijetZMassZSBLowMETWC","M_{ZZ}",400,5.,1005); m_osHistSvc->bookFlavHist(m_histDirZSB+"DijetZMassZSBHighMET","M_{ZZ}",400,5.,1005); m_osHistSvc->bookFlavHist(m_histDirZSB+"DijetZMassZSBHighMETWC","M_{ZZ}",400,5.,1005); /* m_osHistSvc->bookFlavHist(m_histDirZSB+"DijetMassZSBLowMETFake","M_{b#bar{b}} ",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirZSB+"DijetMassZSBHighMETFake","M_{b#bar{b}} ",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirZSB+"DijetMassZSBLowMETWCFake","M_{b#bar{b}} ",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirZSB+"DijetMassZSBHighMETWCFake","M_{b#bar{b}} ",50,0.,250); m_osHistSvc->bookFlavHist(m_histDirZSB+"DijetZMassZSBLowMETFake","M_{ZZ}",400,5.,1005); m_osHistSvc->bookFlavHist(m_histDirZSB+"DijetZMassZSBLowMETWCFake","M_{ZZ}",400,5.,1005); m_osHistSvc->bookFlavHist(m_histDirZSB+"DijetZMassZSBHighMETFake","M_{ZZ}",400,5.,1005); m_osHistSvc->bookFlavHist(m_histDirZSB+"DijetZMassZSBHighMETWCFake","M_{ZZ}",400,5.,1005); */ return; } void OSHbb::BookGenHistos() { // Truth Histograms if (m_doPartonJets) { m_osHistSvc->bookFlavHist(m_histDirTruth+"JetEtGenPart","E_{T}^{jet}",50,0.,250.); m_osHistSvc->bookFlavHist(m_histDirTruth+"JetEtaGenPart","{#eta}^{jet}",70,-3.5,3.5); m_osHistSvc->bookFlavHist(m_histDirTruth+"NPartJet","N^{jet}_{part}",5,0.,5.); m_osHistSvc->bookFlavHist(m_histDirTruth+"Flav1PartJet","N^{jet}_{part}",6,0.,6.); m_osHistSvc->bookFlavHist(m_histDirTruth+"Flav2PartJet","N^{jet}_{part}",6,0.,6.); m_osHistSvc->bookFlavHist(m_histDirTruth+"McPartBJets","",10,0.,10.); } m_osHistSvc->bookFlavHist(m_histDirTruth+"JetEtGen","E_{T}^{jet}",50,0.,250.); m_osHistSvc->bookFlavHist(m_histDirTruth+"JetEtaGen","{#eta}^{jet}",70,-3.5,3.5); m_osHistSvc->bookFlavHist(m_histDirTruth+"NHadJet","N^{jet}_{had}",5,0.,5.); m_osHistSvc->bookFlavHist(m_histDirTruth+"Flav1HadJet","N^{jet}_{had}",6,0.,6.); m_osHistSvc->bookFlavHist(m_histDirTruth+"Flav2HadJet","N^{jet}_{had}",6,0.,6.); m_osHistSvc->bookFlavHist(m_histDirTruth+"McHadBJets","",10,0.,10.); m_osHistSvc->bookFlavHist(m_histDirTruth+"McHadBJetsSysRwEtUp","",10,0.,10.); m_osHistSvc->bookFlavHist(m_histDirTruth+"McHadBJetsSysRwEtDo","",10,0.,10.); m_osHistSvc->bookFlavHist(m_histDirTruth+"McHadBJetsSysRwEtaUp","",10,0.,10.); m_osHistSvc->bookFlavHist(m_histDirTruth+"McHadBJetsSysRwEtaDo","",10,0.,10.); m_osHistSvc->bookFlavHist(m_histDirTruth+"TruthHiggsMass","",400,5.,1005); m_osHistSvc->bookFlavHist(m_histDirTruth+"PtZGen","pt WGen",50,0.,400.); m_osHistSvc->bookFlavHist(m_histDirTruth+"PtLepfromZGen","PtLepfromZGen",50,0.,400.); m_osHistSvc->bookFlavHist(m_histDirTruth+"PtbfromZGen","PtbfromZGen",50,0.,400.); return; } void OSHbb::BookIsolHistos() { // Isolation histograms if (!m_doIsol) return; m_histDirIsol = fTHistStreamName+"Isol/"; m_osHistSvc->bookFlavHist(m_histDirIsol+"MuonCaloIsol","Calo Isol (muon)",50,0.,50.); m_osHistSvc->bookFlavHist(m_histDirIsol+"ElecCaloIsol","Calo Isol (elec)",50,0.,50.); m_osHistSvc->bookFlavHist(m_histDirIsol+"MuonTrackIsol","Track Isol (muon)",50,0.,50.); m_osHistSvc->bookFlavHist(m_histDirIsol+"ElecTrackIsol","Track Isol (elec)",50,0.,50.); m_osHistSvc->bookFlavHist(m_histDirIsol+"MuonCaloIsolRel","Relative Calo Isol (muon)",50,0.,1); m_osHistSvc->bookFlavHist(m_histDirIsol+"ElecCaloIsolRel","Relative Calo Isol (elec)",50,0.,1); m_osHistSvc->bookFlavHist(m_histDirIsol+"MuonTrackIsolRel","Relative Track Isol (muon)",50,0.,1); m_osHistSvc->bookFlavHist(m_histDirIsol+"ElecTrackIsolRel","Relative Track Isol (elec)",50,0.,1); m_osHistSvc->bookFlavHist(m_histDirIsol+"MuonJetDR","Muon-Jet #DeltaR",50,0.,1); m_osHistSvc->bookFlavHist(m_histDirIsol+"ElecJetDR","Elec-Jet #DeltaR",50,0.,1); m_osHistSvc->bookFlavHist(m_histDirIsol+"MuonJetDRWC","Muon-Jet #DeltaR WC",50,0.,1); m_osHistSvc->bookFlavHist(m_histDirIsol+"ElecJetDRWC","Elec-Jet #DeltaR WC",50,0.,1); m_osHistSvc->bookFlavHist(m_histDirIsol+"MuonCaloIsolMZ","Calo Isol (muon)",50,0.,50.); m_osHistSvc->bookFlavHist(m_histDirIsol+"ElecCaloIsolMZ","Calo Isol (elec)",50,0.,50.); m_osHistSvc->bookFlavHist(m_histDirIsol+"MuonTrackIsolMZ","Track Isol (muon)",50,0.,50.); m_osHistSvc->bookFlavHist(m_histDirIsol+"ElecTrackIsolMZ","Track Isol (elec)",50,0.,50.); m_osHistSvc->bookFlavHist(m_histDirIsol+"MuonCaloIsolRelMZ","Relative Calo Isol (muon)",50,0.,1); m_osHistSvc->bookFlavHist(m_histDirIsol+"ElecCaloIsolRelMZ","Relative Calo Isol (elec)",50,0.,1); m_osHistSvc->bookFlavHist(m_histDirIsol+"MuonTrackIsolRelMZ","Relative Track Isol (muon)",50,0.,1); m_osHistSvc->bookFlavHist(m_histDirIsol+"ElecTrackIsolRelMZ","Relative Track Isol (elec)",50,0.,1); return; } void OSHbb::BookCutFlowHistos() { //Book Cut flow histograms static const int ncut =eCutBZMass+1 ; TString cuts[ncut] = {"All","Trig","1e","1Tighte","2e","Only2e","Pte","1mu","1Tightmu","2mu","Only2mu","Ptmu","Zee", "Zmm","Z","emu","Tau","MET","ZMass","CutJetAll","CutJetClean","Dijet","CutExtraJets","Cut3rdJetWC","Cut1stJetWC","DPhiJJ","BDijet","BTag","BZMass"}; //All events BookCutHisto("Cut",ncut, cuts); //Just Zb evnents BookCutHisto("CutZb",ncut, cuts); // Book Cut histos for H->ZZ->llqq static const int nllqqCuts = ellqqCutNum; TString llqqCuts[nllqqCuts] = {"llqqCutAll", "llqqCutTrig", "llqqCutVtx", "llqqCutJetClean", "llqqCut2l", "llqqCutemu", "llqqCutMll", "llqqCutMET", "lqqCutJetVeto", "llqqCutBJetVeto", "llqqCutDiJet", "llqqCutMqq", "llqqCutDPhiqq", "llqqCutDPhill", "llqqCutPt50"}; BookCutHisto("Cutllqq",nllqqCuts, llqqCuts); BookCutHisto("Cuteeqq",nllqqCuts, llqqCuts); BookCutHisto("Cutmumuqq",nllqqCuts, llqqCuts); BookCutHisto("CutllqqNoWeight",nllqqCuts, llqqCuts); BookCutHisto("CuteeqqNoWeight",nllqqCuts, llqqCuts); BookCutHisto("CutmumuqqNoWeight",nllqqCuts, llqqCuts); // Book Cut histos for WH->lnubb static const int nlnubbCuts = elnubbCutNum; TString lnubbCuts[nlnubbCuts] = {"lnubbCutAll", "lnubbCutLepton", "lnubbCutTrig", "lnubbCutVtx", "lnubbCutMETClean", "lnubbCut1eor1mu", "lnubbCutemu", "lnubbCutMT", "lnubbCutMET", "lnubbCutJetVeto", "lnubbCutBJetVeto", "lnubbCutDiJet", "lnubbCutDijetBtag", "lnubbCutDiJet1BJet", "lnubbCutDiJet2BJet", "lnubbCut2Jet", "lnubbCut2Jet1BJet", "lnubbCut2Jet2BJet"}; BookCutHisto("Cutlnubb",nlnubbCuts, lnubbCuts); BookCutHisto("Cutenubb",nlnubbCuts, lnubbCuts); BookCutHisto("Cutmunubb",nlnubbCuts, lnubbCuts); BookCutHisto("CutlnubbNoWeight",nlnubbCuts, lnubbCuts); BookCutHisto("CutenubbNoWeight",nlnubbCuts, lnubbCuts); BookCutHisto("CutmunubbNoWeight",nlnubbCuts, lnubbCuts); // vertex test m_osHistSvc->bookFlavHist(m_histDirCuts+"NTrkVertex","",51,-0.5,50.5); m_osHistSvc->bookFlavHist(m_histDirCuts+"NTrkVertexCut","",51,-0.5,50.5); m_osHistSvc->bookFlavHist(m_histDirCuts+"NVertex2BeforeWeight","",21,-0.5,20.5); m_osHistSvc->bookFlavHist(m_histDirCuts+"NVertex3BeforeWeight","",21,-0.5,20.5); m_osHistSvc->bookFlavHist(m_histDirCuts+"NVertex5BeforeWeight","",21,-0.5,20.5); m_osHistSvc->bookFlavHist(m_histDirCuts+"NVertex2","",21,-0.5,20.5); m_osHistSvc->bookFlavHist(m_histDirCuts+"NVertex3","",21,-0.5,20.5); m_osHistSvc->bookFlavHist(m_histDirCuts+"NVertex5","",21,-0.5,20.5); m_osHistSvc->bookFlavHist(m_histDirCuts+"NVertex2MZBeforeWeight","",21,-0.5,20.5); m_osHistSvc->bookFlavHist(m_histDirCuts+"NVertex3MZBeforeWeight","",21,-0.5,20.5); m_osHistSvc->bookFlavHist(m_histDirCuts+"NVertex5MZBeforeWeight","",21,-0.5,20.5); m_osHistSvc->bookFlavHist(m_histDirCuts+"NVertex2MZ","",21,-0.5,20.5); m_osHistSvc->bookFlavHist(m_histDirCuts+"NVertex3MZ","",21,-0.5,20.5); m_osHistSvc->bookFlavHist(m_histDirCuts+"NVertex3MZDijet","",21,-0.5,20.5); m_osHistSvc->bookFlavHist(m_histDirCuts+"NVertex5MZ","",21,-0.5,20.5); m_osHistSvc->bookFlavHist(m_histDirCuts+"ZVertexBeforeCut","",100,-400.,400.); m_osHistSvc->bookFlavHist(m_histDirCuts+"ZVertexRelBSBeforeCut","",100,-400.,400.); m_osHistSvc->bookFlavHist(m_histDirCuts+"XVertex","",100,-20.,20.); m_osHistSvc->bookFlavHist(m_histDirCuts+"YVertex","",100,-20.,20.); m_osHistSvc->bookFlavHist(m_histDirCuts+"ZVertex","",100,-400.,400.); m_osHistSvc->bookFlavHist(m_histDirCuts+"ZVertexRelBS","",100,-400.,400.); m_osHistSvc->bookFlavHist(m_histDirCuts+"VertexTypeEle","",11,-1.5,9.5); m_osHistSvc->bookFlavHist(m_histDirCuts+"VertexTypeMu" ,"",11,-1.5,9.5); m_osHistSvc->bookFlavHist(m_histDirCuts+"XBeamSpot","",100,-20.,20.); m_osHistSvc->bookFlavHist(m_histDirCuts+"YBeamSpot","",100,-20.,20.); m_osHistSvc->bookFlavHist(m_histDirCuts+"ZBeamSpot","",100,-400.,400.); //Double event histo //m_osHistSvc->bookHist(m_histDirCuts+"EventNumber","Event Number;event number;# Events",10000000,0,10000000); m_triggers[0] = "None"; m_triggers[1] = "EF_e15_loose"; m_triggers[2] = "EF_e10_loose"; m_triggers[3] = "EF_e10_medium"; m_triggers[4] = "EF_mu13"; m_triggers[5] = "EF_2mu6"; m_triggers[6] = "EF_2e5_medium"; m_triggers[7] = "EF_mu6"; m_triggers[8] = "EF_mu10"; m_triggers[9] = "EF_mu13"; TH1F* trigHistElec = new TH1F(m_histDirCuts+TString("TriggersElec"),TString("TriggersElec"),10,0.,10.); TH1F* trigHistMuon = new TH1F(m_histDirCuts+TString("TriggersMuon"),TString("TriggersMuon"),10,0.,10.); for(int i(0); i < 10; i++) { trigHistElec->GetXaxis()->SetBinLabel(i+1,TString(m_triggers[i].c_str())); trigHistMuon->GetXaxis()->SetBinLabel(i+1,TString(m_triggers[i].c_str())); } m_osHistSvc->bookFlavHist(trigHistElec); m_osHistSvc->bookFlavHist(trigHistMuon); return; } void OSHbb::BookCutHisto(TString histoname, int ncut, TString cuts[]){ TH1F* histo = new TH1F(m_histDirCuts+histoname,histoname,ncut,-0.5,ncut-0.5); for(int i=0; iGetXaxis()->SetBinLabel(i+1,cuts[i]); m_osHistSvc->bookFlavHist(histo); return; } void OSHbb::FillllqqCutHisto(int cut, float weight, int flav) { m_osHistSvc->fillTH1Flav(m_histDirCuts+"Cutllqq",cut, weight, flav); m_osHistSvc->fillTH1Flav(m_histDirCuts+"CutllqqNoWeight",cut, 1., flav); if (m_passCutsElectron) { m_osHistSvc->fillTH1Flav(m_histDirCuts+"Cuteeqq",cut, weight, flav); m_osHistSvc->fillTH1Flav(m_histDirCuts+"CuteeqqNoWeight",cut, 1., flav); } if (m_passCutsMuon) { m_osHistSvc->fillTH1Flav(m_histDirCuts+"Cutmumuqq",cut, weight, flav); m_osHistSvc->fillTH1Flav(m_histDirCuts+"CutmumuqqNoWeight",cut, 1., flav); } return; } void OSHbb::FilllnubbCutHisto(int cut, float weight, int flav) { m_osHistSvc->fillTH1Flav(m_histDirCuts+"Cutlnubb",cut, weight, flav); m_osHistSvc->fillTH1Flav(m_histDirCuts+"CutlnubbNoWeight",cut, 1., flav); if (m_passCutsElectron) { m_osHistSvc->fillTH1Flav(m_histDirCuts+"Cutenubb",cut, weight, flav); m_osHistSvc->fillTH1Flav(m_histDirCuts+"CutenubbNoWeight",cut, 1., flav); } if (m_passCutsMuon) { m_osHistSvc->fillTH1Flav(m_histDirCuts+"Cutmunubb",cut, weight, flav); m_osHistSvc->fillTH1Flav(m_histDirCuts+"CutmunubbNoWeight",cut, 1., flav); } return; } StatusCode OSHbb::finalize() { OSHBase::finalize(); ECOutput(); return StatusCode::SUCCESS; } void OSHbb::ECOutput(){ //Output event infomation std::cout << " OSHbb:: Event Counts, Fractions" <0) { for (int i=0; ichronoStart("FillGenHistos"); FillGenHistos(); m_chronoSvc->chronoStop("FillGenHistos"); if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "run number " << m_eventInfo->event_ID()->run_number() << endreq; // Remove duplicate events between streams and fill m_isMC, // if (removeDuplicate()) return StatusCode::SUCCESS; // Reset all event level quantities this->Reset(); // Retrieve verticies const VxContainer* vxContainer(0); if (evtStore()->retrieve(vxContainer,"VxPrimaryCandidate").isFailure() || !vxContainer) { msg(MSG::WARNING) << "Can't retrieve vertex container from StoreGate" << endreq; return StatusCode::SUCCESS; } VxContainer::const_iterator vtxSizeItr(vxContainer->begin()), vtxSizeEnd(vxContainer->end()); int nVtx2(0), nVtx5(0); m_nVtx3=0; m_ntrkVtx=-1; m_primaryVertex = &(* vtxSizeItr)->recVertex(); // the primary vertex itself m_xVtx = m_primaryVertex->position().x(); m_yVtx = m_primaryVertex->position().y(); m_zVtx = m_primaryVertex->position().z(); msg(MSG::DEBUG) << "x vertex " << m_xVtx << endreq; msg(MSG::DEBUG) << "y vertex " << m_yVtx << endreq; msg(MSG::DEBUG) << "z vertex " << m_zVtx << endreq; // loop over and count tracks for (; vtxSizeItr != vtxSizeEnd; ++vtxSizeItr) { const std::vector* tracks = (*vtxSizeItr)->vxTrackAtVertex(); //Is this OK? if (!tracks) continue; if(m_ntrkVtx ==-1) m_ntrkVtx=tracks->size(); if (tracks->size() > 4) nVtx5++; if (tracks->size() > 2) m_nVtx3++; if (tracks->size() > 1) nVtx2++; } // MC Event Weight m_chronoSvc->chronoStart("mcweight"); m_osWeightSvc->setNVertex(m_nVtx3); float event_weight_gen = m_osWeightSvc->getGenWeight(); m_event_weight = m_osWeightSvc->getWeight(); m_chronoSvc->chronoStop("mcweight"); m_event_counter[eECAll] +=m_event_weight; // Before any cuts (but after correct weighting of events) FillllqqCutHisto(ellqqCutAll, event_weight_gen, m_event_flav); FilllnubbCutHisto(elnubbCutAll, event_weight_gen, m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirCuts+"Cut",eCutAll, event_weight_gen,m_event_flav); if(m_zbzlgen) m_osHistSvc->fillTH1Flav(m_histDirCuts+"CutZb", eCutAll, event_weight_gen,m_event_flav); // Find "good" reconstructed particles and do overlap removal. m_chronoSvc->chronoStart("FindReconstructedParticles"); FindReconstructedParticles(); m_chronoSvc->chronoStop("FindReconstructedParticles"); m_passCutEtmiss=true; //Fill Met variables now for Z->nunu and later for others if(m_VectorBosonOp == eZnu) { m_chronoSvc->chronoStart("METCuts"); m_passCutEtmiss = METCuts(true); if (!m_passCutEtmiss) return StatusCode::SUCCESS; m_chronoSvc->chronoStop("METCuts"); } // Apply Electron Cuts and do W/Z mass reconstuction // WRONG - should be after trigger for correct Eta/Pt plots etc. Also, what about getting the weights into here m_chronoSvc->chronoStart("CutsElectron"); m_passCutsElectron = ApplyElectronCuts(true); //We should not use fill muon histos for fake events if(m_llnomFake&&!m_passCutsElectron) return StatusCode::SUCCESS; m_chronoSvc->chronoStop("CutsElectron"); // Apply Muon Cuts and do W/Z mass reconstruction // WRONG - should be after trigger for correct Eta/Pt plots etc. Also, what about getting the weights into here m_chronoSvc->chronoStart("CutsMuon"); m_passCutsMuon = ApplyMuonCuts(true); m_chronoSvc->chronoStop("CutsMuon"); // FillFakeHistos(); int nummu=m_selectMuon->numSelected(); int numele=m_selectElec->numSelected(); if (nummu>0 || numele>0) FilllnubbCutHisto(elnubbCutLepton, m_event_weight, m_event_flav); if (!m_llnomFake) { m_passCutsTrigger=ApplyTriggerCuts(numele,nummu); // also selects on stream } else { // QCD fakes - different triggers for Z/W due to num leptons if (m_VectorBosonOp==eZ) m_passCutsTrigger=ApplyQCDFakeTriggerCuts(1); else m_passCutsTrigger=ApplyQCDFakeTriggerCuts(0); } // Now fakes are done can apply trigger selection for certain stream if (!m_passCutsTrigger) return StatusCode::SUCCESS; // Set correct weights float etamuon[2]; if(nummu>0) etamuon[0]=m_selectMuon->hlvCal(m_mu0).eta(); if(nummu>1) etamuon[1]=m_selectMuon->hlvCal(m_mu1).eta(); float etaelec[2]; if(numele>0) etaelec[0]=m_selectElec->hlvCal(m_e0).eta(); if(numele>1) etaelec[1]=m_selectElec->hlvCal(m_e1).eta(); // m_osWeightSvc->setLeptonTriggerWeight(nummu,etamuon,numele); bool tightele=0; bool combonlymu=0; if (m_VectorBosonOp==eW) { tightele=1; combonlymu=1; } // swiched off in steering m_osWeightSvc->setMuonEfficiencyWeight(nummu, combonlymu); m_osWeightSvc->setElectronEfficiencyWeight(numele, etaelec, tightele); m_event_weight = m_osWeightSvc->getWeight(); FillllqqCutHisto(ellqqCutTrig, m_event_weight, m_event_flav); FilllnubbCutHisto(elnubbCutTrig, m_event_weight, m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirCuts+"Cut",eCutTrig, m_event_weight,m_event_flav); if(m_zbzlgen) m_osHistSvc->fillTH1Flav(m_histDirCuts+"CutZb", eCutTrig, m_event_weight,m_event_flav); // Fill vertex plots HepPoint3D bsPosition = m_beamCondSvc->beamPos(); m_zVtxBS = m_zVtx - bsPosition.z(); m_osHistSvc->fillTH1Flav(m_histDirCuts+"ZBeamSpot",bsPosition.z(), m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirCuts+"XBeamSpot",bsPosition.x(), m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirCuts+"YBeamSpot",bsPosition.y(), m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirCuts+"ZVertexBeforeCut",m_zVtx , m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirCuts+"ZVertexRelBSBeforeCut", m_zVtxBS , m_event_weight,m_event_flav); // Fill N vertex before and after weight (could do this with athena setProperty()) m_osWeightSvc->applyNVertexWeight(false); float eventWeightNoNVtx = m_osWeightSvc->getWeight(); m_osHistSvc->fillTH1Flav(m_histDirCuts+"NVertex2BeforeWeight",nVtx2, eventWeightNoNVtx,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirCuts+"NVertex3BeforeWeight",m_nVtx3, eventWeightNoNVtx,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirCuts+"NVertex5BeforeWeight",nVtx5, eventWeightNoNVtx,m_event_flav); m_osWeightSvc->applyNVertexWeight(true); m_osHistSvc->fillTH1Flav(m_histDirCuts+"NVertex2",nVtx2, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirCuts+"NVertex3",m_nVtx3, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirCuts+"NVertex5",nVtx5, m_event_weight,m_event_flav); // Vertex cut bool vtxcut=1; //if(m_nVtx3<1||m_zVtx<-300.0||m_zVtx>300.0) vtxcut=0; // The z vertex cut kills no events and should not be needed with the beamspot constraint if(m_nVtx3<1) vtxcut=0; if (!vtxcut) return StatusCode::SUCCESS; FillllqqCutHisto(ellqqCutVtx, m_event_weight, m_event_flav); FilllnubbCutHisto(elnubbCutVtx, m_event_weight, m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirCuts+"NTrkVertex",m_ntrkVtx, m_event_weight,m_event_flav); //// Jet cleaning (before fill MET) if(!JetCleaningCut(true)) return StatusCode::SUCCESS; FillllqqCutHisto(ellqqCutJetClean, m_event_weight, m_event_flav); FilllnubbCutHisto(elnubbCutMETClean, m_event_weight, m_event_flav); //Only take events if they have one (W) or two (Z) good electron(s) or muon(s) if (!(m_VectorBosonOp==eZnu)&&!m_passCutsMuon && ! m_passCutsElectron) return StatusCode::SUCCESS; if (m_VectorBosonOp==eZnu&&(!m_passCutsMuon || ! m_passCutsElectron)) return StatusCode::SUCCESS; FillllqqCutHisto(ellqqCut2l, m_event_weight, m_event_flav); FilllnubbCutHisto(elnubbCut1eor1mu, m_event_weight, m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirCuts+"Cut", eCutZ, m_event_weight,m_event_flav); if(m_zbzlgen) m_osHistSvc->fillTH1Flav(m_histDirCuts+"CutZb", eCutZ, m_event_weight,m_event_flav); if(m_passCutsElectron) { m_osHistSvc->fillTH1Flav(m_histDirCuts+"Cut", eCutZee, m_event_weight,m_event_flav); if(m_zbzlgen) m_osHistSvc->fillTH1Flav(m_histDirCuts+"CutZb", eCutZee, m_event_weight,m_event_flav); } if(m_passCutsMuon) { m_osHistSvc->fillTH1Flav(m_histDirCuts+"Cut", eCutZmumu, m_event_weight,m_event_flav); if(m_zbzlgen) m_osHistSvc->fillTH1Flav(m_histDirCuts+"CutZb", eCutZmumu, m_event_weight,m_event_flav); } //For W and Z events reject events if they have leptons of different types if((m_VectorBosonOp==eW||m_VectorBosonOp==eZ)&&m_muons->size()>0&&m_electrons->size()>0) return StatusCode::SUCCESS; //For Z events reject events if they have >=3 leptons of same types if(m_VectorBosonOp==eZ&&(m_muons->size()>2||m_electrons->size()>2)) return StatusCode::SUCCESS; //For Z reject events if they have 2 lepton pairs of different types if(m_VectorBosonOp==eZ&&m_passCutsMuon&&m_passCutsElectron) return StatusCode::SUCCESS; FillllqqCutHisto(ellqqCutemu, m_event_weight, m_event_flav); FilllnubbCutHisto(elnubbCutemu, m_event_weight, m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirCuts+"Cut", eCutemu, m_event_weight,m_event_flav); if(m_zbzlgen) m_osHistSvc->fillTH1Flav(m_histDirCuts+"CutZb", eCutemu, m_event_weight,m_event_flav); m_event_counter[eECOnly1Lep]+= m_event_weight; m_event_counter[eECNoTau]+= m_event_weight; m_osHistSvc->fillTH1Flav(m_histDirCuts+"Cut", eCutTau, m_event_weight,m_event_flav); if(m_zbzlgen) m_osHistSvc->fillTH1Flav(m_histDirCuts+"CutZb", eCutTau, m_event_weight,m_event_flav); //Fill Met variables if(m_VectorBosonOp != eZnu) { m_chronoSvc->chronoStart("METCuts"); m_passCutEtmiss = METCuts(true); //For Z+W analysis we now don't return if event fails missing Et cut so we can look at background from e.g. top +QCD // if (m_VectorBosonOp == eW&&!m_passCutEtmiss) return StatusCode::SUCCESS; m_chronoSvc->chronoStop("METCuts"); } if(m_passCutEtmiss) setFilterPassed(true); float lepangle=180.0*AngleBetweenLeptons()/TMath::Pi(); msg(MSG::DEBUG) << "After call AngleBetweenLeptons()" << endreq; if(lepangle>0) { if(m_passCutsElectron) m_osHistSvc->fillTH1Flav(m_histDirZ+"AngleEE", lepangle, m_event_weight,m_event_flav); if(m_passCutsMuon) m_osHistSvc->fillTH1Flav(m_histDirZ+"AngleMuMu", lepangle, m_event_weight,m_event_flav); } MakeDPhiLL(); // Vertex after MZ ZLeptonMassCut(); // Needs to be called here to calc m_passCutZMass but don't cut on side bands so ignore return if(m_passCutZMass) { m_osHistSvc->fillTH1Flav(m_histDirCuts+"NVertex2MZ", nVtx2, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirCuts+"NVertex3MZ", m_nVtx3, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirCuts+"NVertex5MZ", nVtx5, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirCuts+"NVertex2MZBeforeWeight",nVtx2, eventWeightNoNVtx,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirCuts+"NVertex3MZBeforeWeight",m_nVtx3, eventWeightNoNVtx,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirCuts+"NVertex5MZBeforeWeight",nVtx5, eventWeightNoNVtx,m_event_flav); JetCollection::const_iterator jetItr (m_jets->begin()), jetItrE (m_jets->end()); for (; jetItr != jetItrE; ++jetItr) { Jet* jet = *jetItr; float jvf=jet->getMoment("JVF"); m_osHistSvc->fillTH1Flav(m_histDirJets+"JetVertexFractionEarly", jvf, m_event_weight,m_event_flav); if(jet->eta()<2.5) m_osHistSvc->fillTH1Flav(m_histDirJets+"JetVertexFractionEarlyBarrel", jvf, m_event_weight,m_event_flav); if(jet->eta()<2.1) m_osHistSvc->fillTH1Flav(m_histDirJets+"JetVertexFractionEarlyCentral", jvf, m_event_weight,m_event_flav); else m_osHistSvc->fillTH1Flav(m_histDirJets+"JetVertexFractionEarlyForward", jvf, m_event_weight,m_event_flav); } } for(int iel=0; ielat(iel)->trackParticle()->particleOriginType(); if(vertextype<-1) vertextype=-1; m_osHistSvc->fillTH1Flav(m_histDirCuts+"VertexTypeEle", vertextype, m_event_weight,m_event_flav); } for(int imu=0; imuat(imu)->inDetTrackParticle()->particleOriginType(); if(vertextype<-1) vertextype=-1; m_osHistSvc->fillTH1Flav(m_histDirCuts+"VertexTypeMu", vertextype, m_event_weight,m_event_flav); } //We fill the Systematics here. After the vertex and trigger cuts Fill(); //main filling loop of standard histograms FillSystematics(); //systematics loop return StatusCode::SUCCESS; } void OSHbb::Fill(){ msg(MSG::DEBUG) << "Fill()" << endreq; //Cut on invariant mass of leptons. //This now does not cut out events in the side bands -> only cuts to side band values m_osHistSvc->fillTH1Flav(m_histDirZ+"MZAllMass", m_mZ/GeV, m_event_weight,m_event_flav); if(m_passCutsElectron) m_osHistSvc->fillTH1Flav(m_histDirZ+"MZeeAllMass", m_mZ/GeV, m_event_weight,m_event_flav); if(m_passCutsMuon) m_osHistSvc->fillTH1Flav(m_histDirZ+"MZmumuAllMass", m_mZ/GeV, m_event_weight,m_event_flav); if( !ZLeptonMassCut() ) return; m_osHistSvc->fillTH1Flav(m_histDirCuts+"NTrkVertexCut",m_ntrkVtx, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirCuts+"XVertex",m_xVtx, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirCuts+"YVertex",m_yVtx, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirCuts+"ZVertex",m_zVtx, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirCuts+"ZVertexRelBS", m_zVtxBS , m_event_weight,m_event_flav); //Fill the type of vertex the track is fitted to //enum VertexType { // NoVtx = 0, //!< Dummy vertex, TrackParticle was not used in vertex fit // PriVtx = 1, //!< Primary Vertex // SecVtx = 2, //!< Secondary Vertex // PileUp = 3, //!< Pile Up Vertex // ConvVtx = 4, //!< Converstion Vertex // V0Vtx = 5, //!< Vertex from V0 Decay // KinkVtx = 6, //!< Kink Vertex // V0Lambda = 7, //!< Temporary addition for V0 Lambda // V0LambdaBar = 8, //!< Temporary addition for V0 LambdaBar // V0KShort = 9, //!< Temporary addition for KShort // NotSpecified = -99 //!< this is the default // }; //We call the truth info here after basic Z/W cuts for background processes to speed up code if(m_truthLater) { m_chronoSvc->chronoStart("FindTruthParticles"); FindTruthParticles(); m_chronoSvc->chronoStop("FindTruthParticles"); msg(MSG::DEBUG) << "after call to OSSelecTruth::select()" << endreq; } m_event_counter[eECLepMass]+= m_event_weight; if(m_passCutZMass) FillLeptonHistos(); FillIsolHistos(); if(m_passCutEtmiss&&m_passCutZMass) { float lepangle=180.0*AngleBetweenLeptons()/TMath::Pi(); msg(MSG::DEBUG) << "After call AngleBetweenLeptons()" << endreq; if(lepangle>0&&m_passCutsElectron) m_osHistSvc->fillTH1Flav(m_histDirZ+"AngleEEMZ", lepangle, m_event_weight,m_event_flav); if(m_passCutsMuon) { if(lepangle>0) m_osHistSvc->fillTH1Flav(m_histDirZ+"AngleMuMuMZ", lepangle, m_event_weight,m_event_flav); int nummu=m_selectMuon->numSelected(); for (int i=0; ifillTH1Flav(m_histDirMuon+"MuonD0MZ",m_selectMuon->d0(m_muons->at(i)),m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirMuon+"MuonDZ0ZVertexMZ",m_selectMuon->dZ0ZVertex(m_muons->at(i)),m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirMuon+"MuonIPSigMZ",m_selectMuon->IPSig(m_muons->at(i)),m_event_weight,m_event_flav); } } if(m_passCutsElectron) { int nume=m_selectMuon->numSelected(); for (int i=0; iat(i)->trackParticle()->measuredPerigee()->parameters()[Trk::z0] - m_primaryVertex->position().z()); m_osHistSvc->fillTH1Flav(m_histDirElec+"ElecDZ0ZVertexMZ",dZ0ZVertex,m_event_weight,m_event_flav); } } } if(m_passCutZMass) { m_osHistSvc->fillTH1Flav(m_histDirZ+"PtZ", m_ptZ/GeV, m_event_weight,m_event_flav); FillllqqCutHisto(ellqqCutMll, m_event_weight, m_event_flav); FilllnubbCutHisto(elnubbCutMT, m_event_weight, m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirCuts+"Cut", eCutZMass, m_event_weight,m_event_flav); if(m_zbzlgen) m_osHistSvc->fillTH1Flav(m_histDirCuts+"CutZb", eCutZMass, m_event_weight,m_event_flav); } if(m_passCutEtmiss&&m_passCutZMass) { FillllqqCutHisto(ellqqCutMET, m_event_weight, m_event_flav); FilllnubbCutHisto(elnubbCutMET, m_event_weight, m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirCuts+"Cut", eCutMET, m_event_weight,m_event_flav); if(m_zbzlgen) m_osHistSvc->fillTH1Flav(m_histDirCuts+"CutZb", eCutMET, m_event_weight,m_event_flav); } //Cut on number of jets. This has to be done before FindJets as we //change the jet array //if(!JetCleaningCut(true)) return; if(!JetVertexFractionCut(true)) return; m_chronoSvc->chronoStart("FindJets"); FindJets(); m_chronoSvc->chronoStop("FindJets"); // int numjets= m_selectJet->numSelected(); if(m_passCutsElectron) { if (numjets==0) m_osHistSvc->fillTH1Flav(m_histDirMET+"MET Fake J0 Elec",m_EtMiss/GeV, m_event_weight,m_event_flav); if (numjets==1) m_osHistSvc->fillTH1Flav(m_histDirMET+"MET Fake J1 Elec",m_EtMiss/GeV, m_event_weight,m_event_flav); if (numjets==2) m_osHistSvc->fillTH1Flav(m_histDirMET+"MET Fake J2 Elec",m_EtMiss/GeV, m_event_weight,m_event_flav); if (numjets>=3) m_osHistSvc->fillTH1Flav(m_histDirMET+"MET Fake J3Plus Elec",m_EtMiss/GeV, m_event_weight,m_event_flav); } m_osHistSvc->fillTH1Flav(m_histDirZ+"MZ", m_mZ/GeV, m_event_weight,m_event_flav); if(numjets>1) m_osHistSvc->fillTH1Flav(m_histDirCuts+"NVertex3MZDijet", m_nVtx3, m_event_weight,m_event_flav); if(m_passCutsElectron) { int numele=m_electrons->size(); m_osHistSvc->fillTH1Flav(m_histDirZ+"MZee", m_mZ/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirZ+"MZeeOrig", m_mZOrig/GeV, m_event_weight,m_event_flav); if(numele>1&&m_e0->charge()+m_e1->charge()) m_osHistSvc->fillTH1Flav(m_histDirZ+"MZeeSameQ", m_mZ/GeV, m_event_weight,m_event_flav); //Check to see if either electron is in crack if(numele>0) { float eta0=m_selectElec->hlvCal(m_e0).eta(); bool crackpose=false; bool cracknege=false; if(m_selectElec->isInCrack(eta0)) { if(eta0>0) crackpose=true; else cracknege=true; } if(numele>1) { float eta1=m_selectElec->hlvCal(m_e1).eta(); if(m_selectElec->isInCrack(eta1)) { if(eta1>0) crackpose=true; else cracknege=true; } } if(crackpose) m_osHistSvc->fillTH1Flav(m_histDirZ+"MZeePosCrack", m_mZ/GeV, m_event_weight,m_event_flav); if(cracknege) m_osHistSvc->fillTH1Flav(m_histDirZ+"MZeeNegCrack", m_mZ/GeV, m_event_weight,m_event_flav); if(!cracknege&&!crackpose) m_osHistSvc->fillTH1Flav(m_histDirZ+"MZeeNoCrack", m_mZ/GeV, m_event_weight,m_event_flav); } if (m_trigDec->isPassed("EF_e10_loose")) m_osHistSvc->fillTH1Flav(m_histDirZ+"MZee_e10_loose", m_mZ/GeV, m_event_weight,m_event_flav); if (m_trigDec->isPassed("EF_e10_medium")) m_osHistSvc->fillTH1Flav(m_histDirZ+"MZee_e10_medium", m_mZ/GeV, m_event_weight,m_event_flav); if (m_trigDec->isPassed("EF_e15_loose")) m_osHistSvc->fillTH1Flav(m_histDirZ+"MZee_e15_loose", m_mZ/GeV, m_event_weight,m_event_flav); } if( m_passCutsMuon) { m_osHistSvc->fillTH1Flav(m_histDirZ+"MZmumu", m_mZ/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirZ+"MZmumuOrig", m_mZOrig/GeV, m_event_weight,m_event_flav); if (m_trigDec->isPassed("EF_mu10")) m_osHistSvc->fillTH1Flav(m_histDirZ+"MZmumu_mu10", m_mZ/GeV, m_event_weight,m_event_flav); if (m_trigDec->isPassed("EF_mu13")) m_osHistSvc->fillTH1Flav(m_histDirZ+"MZmumu_mu13", m_mZ/GeV, m_event_weight,m_event_flav); if (m_trigDec->isPassed("EF_mu15")) m_osHistSvc->fillTH1Flav(m_histDirZ+"MZmumu_mu15", m_mZ/GeV, m_event_weight,m_event_flav); if (m_trigDec->isPassed("EF_mu20")) m_osHistSvc->fillTH1Flav(m_histDirZ+"MZmumu_mu20", m_mZ/GeV, m_event_weight,m_event_flav); if(m_muons->size()>1&&m_mu0->charge()+m_mu1->charge()) m_osHistSvc->fillTH1Flav(m_histDirZ+"MZmumuSameQ", m_mZ/GeV, m_event_weight,m_event_flav); } if(!m_passCutEtmiss) { if(m_passCutsElectron) m_osHistSvc->fillTH1Flav(m_histDirZ+"MZeeHighMET", m_mZ/GeV, m_event_weight,m_event_flav); if(m_passCutsMuon) m_osHistSvc->fillTH1Flav(m_histDirZ+"MZmumuHighMET", m_mZ/GeV, m_event_weight,m_event_flav); } double sumet = m_sumEtMiss; //Inclusive Missing Et if (m_passCutsMuon) { m_osHistSvc->fillTH1Flav(m_histDirMET+"METMuon",m_EtMiss/GeV, m_event_weight,m_event_flav); } else if (m_passCutsElectron) { m_osHistSvc->fillTH1Flav(m_histDirMET+"METElec",m_EtMiss/GeV, m_event_weight,m_event_flav); } if(m_passCutZMass) { m_osHistSvc->fillTH1Flav(m_histDirMET+"METMZ",m_EtMiss/GeV, m_event_weight,m_event_flav); int nummu=m_selectMuon->numSelected(); if(m_passCutsElectron) { m_osHistSvc->fillTH1Flav(m_histDirMET+"METMZee", m_EtMiss/GeV, m_event_weight); if(nummu>0) m_osHistSvc->fillTH1Flav(m_histDirMET+"METMZeeExtraMu", m_EtMiss/GeV, m_event_weight); } if(m_passCutsMuon) { m_osHistSvc->fillTH1Flav(m_histDirMET+"METMZmumu", m_EtMiss/GeV, m_event_weight); // watch out for W events - only 1 mu if(nummu>1 && static_cast(m_mu0)->isCombinedMuon() && static_cast(m_mu1)->isCombinedMuon()) m_osHistSvc->fillTH1Flav(m_histDirMET+"METMZmumuCombined", m_EtMiss/GeV, m_event_weight); if(nummu>2) m_osHistSvc->fillTH1Flav(m_histDirMET+"METMZmumuExtraMu", m_EtMiss/GeV, m_event_weight); } } // Missing ET histos for events with at least 1 jet // int numjets= m_selectJet->numSelected(); if ( numjets > 0 ) { m_osHistSvc->fillTH1Flav(m_histDirMET+"HTJetCut", sumet/GeV, m_event_weight); m_osHistSvc->fillTH1Flav(m_histDirMET+"METJetCut", m_EtMiss/GeV, m_event_weight); m_osHistSvc->fillTH2Flav(m_histDirMET+"MET_vs_HTJetCut", sumet/GeV, m_EtMiss, m_event_weight); if (m_passCutsMuon) { m_osHistSvc->fillTH1Flav(m_histDirMET+"METMuonJetCut",m_EtMiss/GeV, m_event_weight,m_event_flav); } else if (m_passCutsElectron) { m_osHistSvc->fillTH1Flav(m_histDirMET+"METElecJetCut",m_EtMiss/GeV, m_event_weight,m_event_flav); } } if ( numjets > 1 ) { m_osHistSvc->fillTH1Flav(m_histDirMET+"METDijetCut", m_EtMiss/GeV, m_event_weight); if (m_passCutsMuon) { m_osHistSvc->fillTH1Flav(m_histDirMET+"METMuonDijetCut",m_EtMiss/GeV, m_event_weight,m_event_flav); } else if (m_passCutsElectron) { m_osHistSvc->fillTH1Flav(m_histDirMET+"METElecDijetCut",m_EtMiss/GeV, m_event_weight,m_event_flav); } } if ( NBJets()>0 && m_passCutThirdJetWeight) { m_osHistSvc->fillTH1Flav(m_histDirMET+"MET1WC", m_EtMiss/GeV, m_event_weight); if (m_passCutsMuon) { m_osHistSvc->fillTH1Flav(m_histDirMET+"METMuon1WC",m_EtMiss/GeV, m_event_weight,m_event_flav); } else if (m_passCutsElectron) { m_osHistSvc->fillTH1Flav(m_histDirMET+"METElec1WC",m_EtMiss/GeV, m_event_weight,m_event_flav); } } if ( NBJets()>1 && m_passCutThirdJetWeight) { m_osHistSvc->fillTH1Flav(m_histDirMET+"METWC", m_EtMiss/GeV, m_event_weight); if (m_passCutsMuon) { m_osHistSvc->fillTH1Flav(m_histDirMET+"METMuonWC",m_EtMiss/GeV, m_event_weight,m_event_flav); } else if (m_passCutsElectron) { m_osHistSvc->fillTH1Flav(m_histDirMET+"METElecWC",m_EtMiss/GeV, m_event_weight,m_event_flav); } } if(m_passCutEtmiss&&m_passCutZMass) { m_chronoSvc->chronoStart("FillJetHistos"); FillHt(); FillJetHistos(); FillJetSystematicsHistos() ; m_chronoSvc->chronoStop("FillJetHistos"); } if(m_passCutZMass && m_passCutsElectron) { FillMETJetHistos(); // for fake rate in eW case } if(!JetCountCut(true)) return ; //Cut on 3rd Jet Weight // For WH cut on this imediately, whereas for ZZ do later to allow for llbb and llqq //if(!ThirdJetWeightCut(true)) return ; m_passCutThirdJetWeight = ThirdJetWeightCut(true); if (m_VectorBosonOp==eW && !m_passCutThirdJetWeight) return; if(m_VectorBoson2Op==eZnu&&!FirstJetWeightCut(true)) return ; if(!BJetsBackToBackCut(true)) return ; //Fill di-jet histograms MakeJJMass(); if(m_passCutEtmiss&&m_passCutZMass) { m_chronoSvc->chronoStart("FillDijetHistos"); FillDijetHistos(); m_chronoSvc->chronoStop("FillDijetHistos"); } //For these histogram we want the side bands also if(m_passCutEtmiss &&m_jet0 && m_jet1) { m_osHistSvc->fillTH1Flav(m_histDirZ+"MZDijet", m_mZ/GeV, m_event_weight,m_event_flav); if(m_passCutsElectron) m_osHistSvc->fillTH1Flav(m_histDirZ+"MZeeDijet", m_mZ/GeV, m_event_weight,m_event_flav); if(m_passCutsMuon) m_osHistSvc->fillTH1Flav(m_histDirZ+"MZmumuDijet", m_mZ/GeV, m_event_weight,m_event_flav); } //Fill H->ZZ histograms m_chronoSvc->chronoStart("MakeZZHiggsMass"); m_HiggsZZMass=ZZHiggsMass(m_bjet0,m_bjet1,1); m_HiggsZZMassNoScale=ZZHiggsMass(m_bjet0,m_bjet1,1); //Higgs mass based on highest Pt jets for 1 or 0 btags m_HiggsZZMassOPt=ZZHiggsMass(m_jet0,m_jet1,0); m_HiggsZZMassOPtNoScale=ZZHiggsMass(m_jet0,m_jet1,0); m_chronoSvc->chronoStop("MakeZZHiggsMass"); if(m_passCutEtmiss&&m_passCutZMass) { m_chronoSvc->chronoStart("FillFinalHistos"); FillHiggsZZHistos(); FillFinalHistos(""); m_chronoSvc->chronoStop("FillFinalHistos"); } int nrecbjet=NBJets(); if( nrecbjet>1 && m_passCutThirdJetWeight) { if(m_passCutEtmiss&&m_passCutZMass) { m_chronoSvc->chronoStart("FillWCHistos"); m_osHistSvc->fillTH1Flav(m_histDirCuts+"Cut", eCutBTag, m_event_weight,m_event_flav); if(m_zbzlgen) m_osHistSvc->fillTH1Flav(m_histDirCuts+"CutZb", eCutBTag, m_event_weight,m_event_flav); FilllnubbCutHisto(elnubbCutDijetBtag, m_event_weight, m_event_flav); FillDijetHistosAfterBTag(); FillHiggsZZHistosAfterBTag(); m_chronoSvc->chronoStop("FillWCHistos"); } if(m_passCutEtmiss) { m_osHistSvc->fillTH1Flav(m_histDirZ+"MZDijetWC", m_mZ/GeV, m_event_weight); if(m_passCutsElectron) m_osHistSvc->fillTH1Flav(m_histDirZ+"MZeeDijetWC", m_mZ/GeV, m_event_weight,m_event_flav); if(m_passCutsMuon) m_osHistSvc->fillTH1Flav(m_histDirZ+"MZmumuDijetWC", m_mZ/GeV, m_event_weight,m_event_flav); } } if(nrecbjet==1 &&m_passCutEtmiss) { m_osHistSvc->fillTH1Flav(m_histDirZ+"MZDijet1WC", m_mZ/GeV, m_event_weight); } //Fill Z side band histos if(m_VectorBosonOp==eZ&&!m_passCutZMass) { m_chronoSvc->chronoStart("FillZSBHistos"); FillZSBHistos(); m_chronoSvc->chronoStop("FillZSBHistos"); } return ; } bool OSHbb::CutsPhoton() { msg(MSG::DEBUG) << "OSHbb::CutsPhoton()" << endreq; m_selectPhot->select(); if(m_selectPhot->numSelected()>0) { std::cout <<" Photon found " << std::endl; return kFALSE; } return kTRUE; } bool OSHbb::WCutsElectron(bool fillhistos) { msg(MSG::DEBUG) << "OSHbb::WCutsElectron()" << endreq; //Count events with >0,>0 tight, and >1 electrons. if(fillhistos) { if(m_selectElec->numSelected()>0) m_event_counter[eEC1GoodE]+= m_event_weight; if(m_selectElec->numSelectedType(OSSelectorsEnums::ElecTight)+ m_selectElec->numSelectedType(OSSelectorsEnums::ElectronTightIso)>0) m_event_counter[eEC1TightE]+= m_event_weight; if(m_selectElec->numSelected()>1) m_event_counter[eEC2GoodE]+= m_event_weight; } //Cut on number of tight electrons. // if(m_selectElec->numSelectedType(OSSelectorsEnums::ElecTight) // +m_selectElec->numSelectedType(OSSelectorsEnums::ElecTightIso)numSelectedType(OSSelectorsEnums::ElectronTightIso)numSelected()numSelected()numSelected()>m_nelecmax_eW) return 0; // fill Pt, eta for first electron only // Calculate the transverse mass Electron* electron=m_electrons->at(0); //Cut on highest electron pt. float ptlep = m_selectElec->hlvCal(electron).perp(); if(ptlephlvCal(electron).eta(); if(fillhistos) { m_osHistSvc->fillTH1Flav(m_histDirElec+"ElecPt", ptlep/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirElec+"ElecEta", etalep, m_event_weight,m_event_flav); } float pxlep = m_selectElec->hlvCal(electron).px(); float pylep =m_selectElec->hlvCal(electron).py(); float mT = sqrt(pow(m_EtMiss+ptlep,2) - pow(m_EtMiss_x+pxlep,2) - pow(m_EtMiss_y+pylep,2)); m_ptZ=sqrt(pow(m_EtMiss_x+pxlep,2)+pow(m_EtMiss_y+pylep,2)); m_phiZ=atan2(m_EtMiss_y+pylep,m_EtMiss_x+pxlep); m_mZ=mT; //Passed cuts. if(fillhistos) m_event_counter[eECAllElectronCuts]+=m_event_weight; return 1; } bool OSHbb::ZCutsElectron(bool fillhistos) { msg(MSG::DEBUG) << "OSHbb::ZCutsElectron()" << endreq; if(fillhistos&&m_selectElec->numSelected()>0) { m_event_counter[eEC1GoodE]+= m_event_weight; m_osHistSvc->fillTH1Flav(m_histDirCuts+"Cut", eCut1e, m_event_weight,m_event_flav); if(m_zbzlgen) m_osHistSvc->fillTH1Flav(m_histDirCuts+"CutZb", eCut1e, m_event_weight,m_event_flav); // std::cout << " numloose " << m_selectElec->numSelectedType(OSSelectorsEnums::ElecLoose) << " nummedium " << m_selectElec->numSelectedType(OSSelectorsEnums::ElecMedium) << " numtight " << m_selectElec->numSelectedType(OSSelectorsEnums::ElecTight) << " numall " << m_selectElec->numSelected() <size(); ++i) { OSSelectorsEnums::EmType etype = m_selectElec->emType(m_electrons->at(i)); if(etype >= OSSelectorsEnums::ElecMedium) { ++numTight; } } if(fillhistos&&numTight>0) { m_event_counter[eEC1TightE]+= m_event_weight; m_osHistSvc->fillTH1Flav(m_histDirCuts+"Cut", eCut1Tighte, m_event_weight,m_event_flav); if(m_zbzlgen) m_osHistSvc->fillTH1Flav(m_histDirCuts+"CutZb", eCut1Tighte, m_event_weight,m_event_flav); } // Cut on number of electrons if(m_selectElec->numSelected()<2) return 0; if(fillhistos){ m_osHistSvc->fillTH1Flav(m_histDirCuts+"Cut", eCut2e, m_event_weight,m_event_flav); if(m_zbzlgen) m_osHistSvc->fillTH1Flav(m_histDirCuts+"CutZb", eCut2e, m_event_weight,m_event_flav); m_event_counter[eEC2GoodE]+= m_event_weight; } //Cut on number of medium+tight electrons. //CBG: kills LLnoM method //if(numTightnumSelected()>m_nelecmax_eZ) return 0; if(fillhistos){ m_osHistSvc->fillTH1Flav(m_histDirCuts+"Cut", eCutOnly2e, m_event_weight,m_event_flav); if(m_zbzlgen) m_osHistSvc->fillTH1Flav(m_histDirCuts+"CutZb", eCutOnly2e, m_event_weight,m_event_flav); } //Cut on highest Pt electron Pt. if(m_electrons->at(0)->et()at(1)->et()at(0); Electron* elec1 = m_electrons->at(1); if(fillhistos){ m_osHistSvc->fillTH1Flav(m_histDirCuts+"Cut", eCutPte, m_event_weight,m_event_flav); if(m_zbzlgen) m_osHistSvc->fillTH1Flav(m_histDirCuts+"CutZb", eCutPte, m_event_weight,m_event_flav); // fill pt, eta for first electron only m_osHistSvc->fillTH1Flav(m_histDirElec+"ElecPt", (m_selectElec->hlvCal(elec0)).perp()/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirElec+"ElecEta", (m_selectElec->hlvCal(elec0)).eta(), m_event_weight,m_event_flav); for (int i=0; i<2; i++){ float dZ0ZVertex = fabs(m_electrons->at(i)->trackParticle()->measuredPerigee()->parameters()[Trk::z0] - m_primaryVertex->position().z()); m_osHistSvc->fillTH1Flav(m_histDirElec+"ElecDZ0ZVertex",dZ0ZVertex,m_event_weight,m_event_flav); } } // work out the mass //float pxlep = m_electrons->at(0)->px()+m_electrons->at(1)->px(); //float pylep = m_electrons->at(0)->py()+m_electrons->at(1)->py(); // m_mZ = (m_electrons->at(0)->hlv()+m_electrons->at(1)->hlv()).m(); // OSElectronUserDataContainer* udElecs = m_selectElec->selectedElecsUD(); // m_mZ = ((*udElecs)[elec0]->fourMom()+(*udElecs)[elec1]->fourMom()).m(); HepLorentzVector Z=m_selectElec->hlvCal(elec0)+m_selectElec->hlvCal(elec1); m_mZ = Z.m(); m_mZOrig= (elec0->hlv()+elec1->hlv()).m(); m_ptZ = Z.perp(); m_phiZ = Z.phi(); //Passed cuts. if(fillhistos) m_event_counter[eECAllElectronCuts]+=m_event_weight; return 1; } bool OSHbb::ZNuCutsElectron(bool fillhistos) { msg(MSG::DEBUG) << "OSHbb::ZNuCutsElectron()" << endreq; //Count events by lepton numbers and type. if(fillhistos) { if(m_selectElec->numSelected()>0) m_event_counter[eEC1GoodE]+= m_event_weight; if(m_selectElec->numSelected()>1) m_event_counter[eEC2GoodE]+= m_event_weight; if(m_selectElec->numSelectedType(OSSelectorsEnums::ElecTight)>0) m_event_counter[eEC1TightE]+= m_event_weight; } //Cut on max number of electrons. if(m_selectElec->numSelected()>m_nelecmax_eZnu) return 0; //Passes cut. if(fillhistos) m_event_counter[eECAllElectronCuts]+=m_event_weight; return 1; } bool OSHbb::WCutsMuon(bool fillhistos) { msg(MSG::DEBUG) << " OSHbb::WCutsMuon()" << endreq; //Fill event counters. if(fillhistos) { if(m_selectMuon->numSelected()>0) m_event_counter[eEC1GoodMu]+= m_event_weight; if(m_selectMuon->numSelected()>0) m_event_counter[eEC1TightMu]+= m_event_weight; if(m_selectMuon->numSelected()>1) m_event_counter[eEC2GoodMu]+= m_event_weight; } //Cut on minumum number of selected muons if(m_selectMuon->numSelected()numSelected()>m_nmumax_eW) return 0; //Cut on highest muon pt. Muon* muon=m_muons->at(0); if(m_selectMuon->hlvCal(muon).perp()IsIsolated(m_muons->at(0))) return 0; // work out the transverse mass float ptlep = m_selectMuon->hlvCal(muon).perp(); float etalep = m_selectMuon->hlvCal(muon).eta(); // fill pt, eta for first muon only if(fillhistos) { m_osHistSvc->fillTH1Flav(m_histDirMuon+"MuonPt", ptlep/GeV); m_osHistSvc->fillTH1Flav(m_histDirMuon+"MuonEta", etalep); m_osHistSvc->fillTH1Flav(m_histDirMuon+"MuonChi2", m_muons->at(0)->matchChi2OverDoF()); } float pxlep = m_selectMuon->hlvCal(muon).px(); float pylep = m_selectMuon->hlvCal(muon).py(); float mT = sqrt(pow(m_EtMiss+ptlep,2) - pow(m_EtMiss_x+pxlep,2) - pow(m_EtMiss_y+pylep,2)); m_ptZ=sqrt(pow(m_EtMiss_x+pxlep,2)+pow(m_EtMiss_y+pylep,2)); m_phiZ=atan2(m_EtMiss_y+pylep,m_EtMiss_x+pxlep); m_mZ=mT; //Passes muon cuts. if(fillhistos) m_event_counter[eECAllMuonCuts]+=m_event_weight; return 1; } bool OSHbb::ZCutsMuon( bool fillhistos) { msg(MSG::DEBUG) << " OSHbb::ZCutsMuon()" << endreq; m_muons=m_selectMuon->selectedMuons(); if(fillhistos&&m_selectMuon->numSelected()>0) { m_event_counter[eEC1GoodMu]+= m_event_weight; m_osHistSvc->fillTH1Flav(m_histDirCuts+"Cut", eCut1mu, m_event_weight,m_event_flav); if(m_zbzlgen) m_osHistSvc->fillTH1Flav(m_histDirCuts+"CutZb", eCut1mu, m_event_weight,m_event_flav); } if( fillhistos && m_selectMuon->numSelected()>0) { m_event_counter[eEC1TightMu]+= m_event_weight; m_osHistSvc->fillTH1Flav(m_histDirCuts+"Cut", eCut1Tightmu, m_event_weight,m_event_flav); if(m_zbzlgen) m_osHistSvc->fillTH1Flav(m_histDirCuts+"CutZb", eCut1Tightmu, m_event_weight,m_event_flav); } // Cut on number of selected muons if(m_selectMuon->numSelected()<2) return 0; // Cut on number of selected combined muons (allows 1 comb + 1 low Pt) //if(m_selectMuon->numSelectedCombined()<1) return 0; if(fillhistos) { m_osHistSvc->fillTH1Flav(m_histDirCuts+"Cut", eCut2mu, m_event_weight,m_event_flav); if(m_zbzlgen) m_osHistSvc->fillTH1Flav(m_histDirCuts+"CutZb", eCut2mu, m_event_weight,m_event_flav); m_event_counter[eEC2GoodMu]+= m_event_weight; } //Cut on maximum number of selected muons if(m_selectMuon->numSelected()>m_nmumax_eZ) return 0; if(fillhistos) { m_osHistSvc->fillTH1Flav(m_histDirCuts+"Cut", eCutOnly2mu, m_event_weight,m_event_flav); if(m_zbzlgen) m_osHistSvc->fillTH1Flav(m_histDirCuts+"CutZb", eCutOnly2mu, m_event_weight,m_event_flav); } Muon* muon0=m_muons->at(0); Muon* muon1=m_muons->at(1); //Cut on highest muon pt if(m_selectMuon->hlvCal(muon0).perp()hlvCal(muon1).perp()fillTH1Flav(m_histDirCuts+"Cut", eCutPtmu, m_event_weight,m_event_flav); if(m_zbzlgen) m_osHistSvc->fillTH1Flav(m_histDirCuts+"CutZb", eCutPtmu, m_event_weight,m_event_flav); // fill pt, eta for first muon only m_osHistSvc->fillTH1Flav(m_histDirMuon+"MuonPt", m_selectMuon->hlvCal(muon0).perp()/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirMuon+"MuonEta", m_selectMuon->hlvCal(muon0).eta(), m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirMuon+"MuonChi2", m_muons->at(0)->matchChi2OverDoF(), m_event_weight,m_event_flav); } // Fill muon qualtiy for 1st and 2nd muons if(fillhistos) { int muonQuality1(0), muonQuality2(0); if (m_muons->at(0)->isLoose()) muonQuality1 = 1; if (m_muons->at(0)->isMedium()) muonQuality1 = 2; if (m_muons->at(0)->isTight()) muonQuality1 = 3; if (m_muons->at(1)->isLoose()) muonQuality2 = 1; if (m_muons->at(1)->isMedium()) muonQuality2 = 2; if (m_muons->at(1)->isTight()) muonQuality2 = 3; m_osHistSvc->fillTH1Flav(m_histDirMuon+"MuonQuality", muonQuality1, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirMuon+"MuonQuality", muonQuality2, m_event_weight,m_event_flav); for (int i=0; i<2; i++){ m_osHistSvc->fillTH1Flav(m_histDirMuon+"MuonD0",m_selectMuon->d0(m_muons->at(i)),m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirMuon+"MuonDZ0ZVertex",m_selectMuon->dZ0ZVertex(m_muons->at(i)),m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirMuon+"MuonIPSig",m_selectMuon->IPSig(m_muons->at(i)),m_event_weight,m_event_flav); } } // work out the mass etc. // Require muons to be opposite charged if (muon0->charge()*muon1->charge() != -1) return false; HepLorentzVector Z= m_selectMuon->hlvCal(muon0)+m_selectMuon->hlvCal(muon1); m_mZ = Z.m(); m_mZOrig= (muon0->hlv()+muon1->hlv()).m(); m_ptZ = Z.perp(); m_phiZ = Z.phi(); //Passes muon cuts. if(fillhistos) m_event_counter[eECAllMuonCuts]+=m_event_weight; return 1; } bool OSHbb::ZNuCutsMuon(bool fillhistos ) { msg(MSG::DEBUG) << " OSHbb::ZnuCutsMuon()" << endreq; if(fillhistos) { //Fill event counters if(m_selectMuon->numSelected()>0) m_event_counter[eEC1GoodMu]+= m_event_weight; if(m_selectMuon->numSelected()>1) m_event_counter[eEC2GoodMu]+= m_event_weight; if(m_selectMuon->numSelected()>0) m_event_counter[eEC1TightMu]+= m_event_weight; } //Cut on max number of muons. //We only are interested in events with 2 or fewer muons if(m_selectMuon->numSelected()>=3) return 0; // Work out the mass etc. m_mZ = 0; m_ptZ = m_EtMiss; m_phiZ = atan2(-m_EtMiss_y,-m_EtMiss_x); //Passes cut. if(fillhistos) m_event_counter[eECAllMuonCuts]+=m_event_weight; return 1; } bool OSHbb::METCuts(bool fillhistos) { msg(MSG::DEBUG) << "METCuts" << endreq; if(fillhistos) { m_osHistSvc->fillTH1Flav(m_histDirMET+"METOrig",m_EtMissOrig/GeV, m_event_weight,m_event_flav); // Fill individual terms const MissingET * pMissing = 0; if (!evtStore()->retrieve(pMissing, "MET_RefEle").isFailure()) { m_osHistSvc->fillTH1Flav(m_histDirMET+"METEle", pMissing->et()/GeV, m_event_weight,m_event_flav); } if (!evtStore()->retrieve(pMissing, "MET_RefGamma").isFailure()) { m_osHistSvc->fillTH1Flav(m_histDirMET+"METGamma", pMissing->et()/GeV, m_event_weight,m_event_flav); } if (!evtStore()->retrieve(pMissing, "MET_RefTau").isFailure()){ m_osHistSvc->fillTH1Flav(m_histDirMET+"METTau", pMissing->et()/GeV, m_event_weight,m_event_flav); } if (!evtStore()->retrieve(pMissing, "MET_RefJet").isFailure()){ m_osHistSvc->fillTH1Flav(m_histDirMET+"METJet", pMissing->et()/GeV, m_event_weight,m_event_flav); } if (!evtStore()->retrieve(pMissing, "MET_RefMuon").isFailure()){ m_osHistSvc->fillTH1Flav(m_histDirMET+"METRefMuon", pMissing->et()/GeV, m_event_weight,m_event_flav); } if (!evtStore()->retrieve(pMissing, "MET_CellOut").isFailure()){ m_osHistSvc->fillTH1Flav(m_histDirMET+"METCellOut", pMissing->et()/GeV, m_event_weight,m_event_flav); } if (!evtStore()->retrieve(pMissing, "MET_Cryo").isFailure()){ m_osHistSvc->fillTH1Flav(m_histDirMET+"METCryo", pMissing->et()/GeV, m_event_weight,m_event_flav); } if (!evtStore()->retrieve(pMissing, "MET_MuonBoy").isFailure()){ m_osHistSvc->fillTH1Flav(m_histDirMET+"METMuonBoy", pMissing->et()/GeV, m_event_weight,m_event_flav); } // Fill truth const MissingEtTruth * pMissingTruth = 0; if (evtStore()->contains("MET_Truth") && !evtStore()->retrieve(pMissingTruth, "MET_Truth").isFailure() && pMissingTruth){ double exMissTrueNonInt = pMissingTruth->exTruth(MissingEtTruth::NonInt); double eyMissTrueNonInt = pMissingTruth->eyTruth(MissingEtTruth::NonInt); double etMissTrueNonInt = std::sqrt(exMissTrueNonInt * exMissTrueNonInt + eyMissTrueNonInt * eyMissTrueNonInt); double exMissTrueIntCent = pMissingTruth->exTruth(MissingEtTruth::IntCentral); double eyMissTrueIntCent = pMissingTruth->eyTruth(MissingEtTruth::IntCentral); double etMissTrueIntCent = std::sqrt(exMissTrueIntCent * exMissTrueIntCent + eyMissTrueIntCent * eyMissTrueIntCent); double exMissTrueMuon = pMissingTruth->exTruth(MissingEtTruth::Muons); double eyMissTrueMuon = pMissingTruth->eyTruth(MissingEtTruth::Muons); double etMissTrueMuon = std::sqrt(exMissTrueMuon * exMissTrueMuon + eyMissTrueMuon * eyMissTrueMuon); double etMissTrueSum = std::sqrt(std::pow(exMissTrueNonInt + exMissTrueIntCent + exMissTrueMuon, 2) + std::pow(eyMissTrueNonInt + eyMissTrueIntCent + eyMissTrueMuon, 2)); m_osHistSvc->fillTH1Flav(m_histDirMET+"METTruthNonInt", etMissTrueNonInt/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirMET+"METTruthIntCent", etMissTrueIntCent/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirMET+"METTruthMuon", etMissTrueMuon/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirMET+"METTruthSum", etMissTrueSum/GeV, m_event_weight,m_event_flav); } m_osHistSvc->fillTH1Flav(m_histDirMET+"HT", m_sumEtMiss/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirMET+"MET",m_EtMiss/GeV, m_event_weight,m_event_flav); if (!m_semiLeptonicBDecay && !m_semiLeptonicCDecay) { m_osHistSvc->fillTH1Flav(m_histDirMET+"METNonSemiLep",m_EtMiss/GeV, m_event_weight,m_event_flav); } else { m_osHistSvc->fillTH1Flav(m_histDirMET+"METSemiLep",m_EtMiss/GeV, m_event_weight,m_event_flav); } m_osHistSvc->fillTH2Flav(m_histDirMET+"MET_vs_HT", m_sumEtMiss,m_EtMiss/GeV, m_event_weight,m_event_flav); // Jet-MET angle float dPhiJetMetClosest(999); if (m_EtMiss > 50*GeV) { float metPhi = std::atan2(m_EtMiss_y,m_EtMiss_x); JetCollection::const_iterator jetItr (m_jets->begin()), jetItrE (m_jets->end()); for (; jetItr != jetItrE; ++jetItr) { Jet* jet = *jetItr; float dPhiJetMet = acos(cos(jet->phi() - metPhi))*180/3.14; if (dPhiJetMet < dPhiJetMetClosest) dPhiJetMetClosest = dPhiJetMet; m_osHistSvc->fillTH1Flav(m_histDirMET+"JetMETPhi", dPhiJetMet, m_event_weight); m_osHistSvc->fillTH2Flav(m_histDirMET+"JetMETPhi2D", dPhiJetMet, m_EtMiss, m_event_weight); } } if (dPhiJetMetClosest < 999) m_osHistSvc->fillTH2Flav(m_histDirMET+"JetMETPhi2DClosest", dPhiJetMetClosest, m_EtMiss, m_event_weight); } //Cut on Et miss. if(m_VectorBoson2Op!=eZnu) { if(m_VectorBosonOp==eW&&m_EtMissm_etmisscut_eZ) return 0; else if(m_VectorBosonOp==eZnu&&m_EtMisssize()); Electron* elec0 = (nelec>0 ? m_electrons->at(0) : static_cast(0)); Electron* elec1 = (nelec>1 ? m_electrons->at(1) : static_cast(0)); if (elec0) { HepLorentzVector hlve0 = m_selectElec->hlvCal(elec0); m_osHistSvc->fillTH1Flav(m_histDirElec+"ElecPtAfterTrig", hlve0.perp()/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirElec+"ElecEtaAfterTrig", hlve0.eta(), m_event_weight,m_event_flav); } if (elec0 && elec1) { HepLorentzVector hlve0 = m_selectElec->hlvCal(elec0); HepLorentzVector hlve1 = m_selectElec->hlvCal(elec1); m_osHistSvc->fillTH1Flav(m_histDirElec+"ElecPtBothAfterTrig", hlve0.perp()/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirElec+"ElecPtBothAfterTrig", hlve1.perp()/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirElec+"ElecEtaBothAfterTrig", hlve0.eta(), m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirElec+"ElecEtaBothAfterTrig", hlve1.eta(), m_event_weight,m_event_flav); } unsigned int nmuon(m_muons->size()); Muon* muon0 = (nmuon>0 ? m_muons->at(0) : static_cast(0)); Muon* muon1 = (nmuon>1 ? m_muons->at(1) : static_cast(0)); if (muon0) { HepLorentzVector hlvm0 = m_selectMuon->hlvCal(muon0); m_osHistSvc->fillTH1Flav(m_histDirMuon+"MuonPtAfterTrig", hlvm0.perp()/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirMuon+"MuonEtaAfterTrig", hlvm0.eta(), m_event_weight,m_event_flav); } if (muon0 && muon1) { HepLorentzVector hlvm0 = m_selectMuon->hlvCal(muon0); HepLorentzVector hlvm1 = m_selectMuon->hlvCal(muon1); m_osHistSvc->fillTH1Flav(m_histDirMuon+"MuonPtBothAfterTrig", hlvm0.perp()/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirMuon+"MuonPtBothAfterTrig", hlvm1.perp()/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirMuon+"MuonEtaBothAfterTrig", hlvm0.eta(), m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirMuon+"MuonEtaBothAfterTrig", hlvm1.eta(), m_event_weight,m_event_flav); } return; } void OSHbb::FillZSBHistos() { msg(MSG::DEBUG) << "FillZSBHistos()" << endreq; m_osHistSvc->fillTH1Flav(m_histDirZSB+"METZSB", m_EtMiss/GeV, m_event_weight); //We have difficulty here since jets can be different flavours. //If that is the case set to zero int flavuse = m_flav0==m_flav1 ? m_flav0 : 0; if( m_passCutEtmiss) { m_osHistSvc->fillTH1Flav(m_histDirZSB+"DijetMassZSBLowMET", m_dibmassopt/GeV,m_event_weight,flavuse); m_osHistSvc->fillTH1Flav(m_histDirZSB+"DijetZMassZSBLowMET", m_HiggsZZMassOPt/GeV,m_event_weight,flavuse); } else { m_osHistSvc->fillTH1Flav(m_histDirZSB+"DijetMassZSBHighMET", m_dibmassopt/GeV,m_event_weight,flavuse); m_osHistSvc->fillTH1Flav(m_histDirZSB+"DijetZMassZSBHighMET", m_HiggsZZMassOPt/GeV,m_event_weight,flavuse); } //Need to fill WC histos if(NBJets()>1 && m_passCutThirdJetWeight) { m_osHistSvc->fillTH1Flav(m_histDirZSB+"METZSBWC", m_EtMiss/GeV, m_event_weight); if( m_passCutEtmiss) { m_osHistSvc->fillTH1Flav(m_histDirZSB+"DijetMassZSBLowMETWC", m_dibmassopt/GeV,m_event_weight,flavuse); m_osHistSvc->fillTH1Flav(m_histDirZSB+"DijetZMassZSBLowMETWC", m_HiggsZZMassOPt/GeV,m_event_weight,flavuse); } else { m_osHistSvc->fillTH1Flav(m_histDirZSB+"DijetMassZSBHighMETWC", m_dibmassopt/GeV,m_event_weight,flavuse); m_osHistSvc->fillTH1Flav(m_histDirZSB+"DijetZMassZSBHighMETWC", m_HiggsZZMassOPt/GeV,m_event_weight,flavuse); } } } void OSHbb::FillJetHistos() { msg(MSG::DEBUG) << "FillJetHistos()" << endreq; m_JJDeltaPhi=0; msg(MSG::DEBUG) << "Get num Good jets" << endreq; int numjets= m_selectJet->numSelected(); if(numjets>=2) m_osHistSvc->fillTH1Flav(m_histDirCuts+"Cut", eCutDijet, m_event_weight,m_event_flav); if(m_zbzlgen&&numjets>=2) m_osHistSvc->fillTH1Flav(m_histDirCuts+"CutZb", eCutDijet, m_event_weight,m_event_flav); // Fill the mass histograms with weight 1 to be consistent with non-cut versions // cut against other wrong lepton ID here but these events make it in the jets static std::string histDirE(fTHistStreamName+"Electron/"); msg(MSG::DEBUG) << "Get Good jets" << endreq; m_nlightjet=0; m_ncjet=0; m_nbjet=0; // store some values for the 1st, 2nd and 3rd highest et jets double etj0=0; double etj1=0; double etj2=0; double rweta0up=1.0; double rweta0do=1.0; double rwet0up=1.0; double rwet0do=1.0; bool highestetjetmatchedtohad=0; msg(MSG::DEBUG) << "jet loop" << endreq; JetCollection::const_iterator jetItr (m_jets->begin()), jetItrE (m_jets->end()); // loop over the AOD jet container for (; jetItr != jetItrE; ++jetItr) { Jet* jet = *jetItr; double etj = m_selectJet->hlvCal(jet).perp(); double etaj = jet->eta(); int flav=GetFlavourHad(jet); // number of jets if (flav == 5) m_nbjet++; else if (flav == 4) m_ncjet++; else m_nlightjet++; m_osHistSvc->fillTH1Flav(m_histDirJets+"JetEt", etj/GeV, m_event_weight, flav); m_osHistSvc->fillTH1Flav(m_histDirJets+"JetEta", etaj, m_event_weight, flav); // weight for combination of IP3D and SV1 double jetweight = (m_btaggerName.compare("") != 0) ? jet-> getFlavourTagWeight(m_btaggerName) : jet-> getFlavourTagWeight(); // Jet Fitter double w_jetweight_standard = jet->getFlavourTagWeight(); m_osHistSvc->fillTH1Flav(m_histDirJets+"JetWeight", w_jetweight_standard, m_event_weight, flav); double w_jetfit = jet->getFlavourTagWeight("JetFitterTag"); double w_jetfitcomb = jet->getFlavourTagWeight("JetFitterCOMB"); double w_jetfitnn = jet->getFlavourTagWeight("JetFitterTagNN"); double w_jetfitcombnn = jet->getFlavourTagWeight("JetFitterCOMBNN"); m_osHistSvc->fillTH1Flav(m_histDirJets+"JetFitterWeight", w_jetfit , m_event_weight, flav); m_osHistSvc->fillTH1Flav(m_histDirJets+"JetFitterCombWeight", w_jetfitcomb , m_event_weight, flav); m_osHistSvc->fillTH1Flav(m_histDirJets+"JetFitterTagNNWeight", w_jetfitnn , m_event_weight, flav); m_osHistSvc->fillTH1Flav(m_histDirJets+"JetFitterCOMBNNWeight", w_jetfitcombnn , m_event_weight, flav); if (numjets==1) { m_osHistSvc->fillTH1Flav(m_histDirJets+"JetWeightNJet1", w_jetweight_standard, m_event_weight, flav); m_osHistSvc->fillTH1Flav(m_histDirJets+"JetFitterWeightNJet1", w_jetfit , m_event_weight, flav); } if (numjets==2) { m_osHistSvc->fillTH1Flav(m_histDirJets+"JetWeightNJet2", w_jetweight_standard, m_event_weight, flav); m_osHistSvc->fillTH1Flav(m_histDirJets+"JetFitterWeightNJet2", w_jetfit , m_event_weight, flav); } //systematics //Get nearest hadronic jet and find its Et and eta double hadjetet=0; double hadjeteta=-99; double rechadjetdist=99; double rwetaup=1.0; double rwetado=1.0; double rwetup=1.0; double rwetdo=1.0; // flags to see how many had jets are matched to this rec jet int nmatchedjets=0; JetCollection* hadjetslowpt = m_selectJetlowpt->selectedHadJets(); if (hadjetslowpt) { JetCollection::const_iterator jetItrH (hadjetslowpt->begin()), jetItrHE (hadjetslowpt->end()); for (; jetItrH != jetItrHE; ++jetItrH) { Jet* hadjet = *jetItrH; hadjetet = hadjet->pt(); hadjeteta = hadjet->eta(); // match to had jets with Et > 8 GeV if(hadjetet<8*GeV) continue; rechadjetdist = m_analysisTools->deltaR(jet,hadjet); if(rechadjetdist<0.2) { nmatchedjets++; rwetaup= pow(fabs(hadjeteta)+m_RwEtaBCen,m_RwEtaBError); rwetado= pow(fabs(hadjeteta)+m_RwEtaBCen,-m_RwEtaBError); rwetup= pow((hadjetet/m_RwEtBCen),m_RwEtBError); rwetdo= pow((hadjetet/m_RwEtBCen),-m_RwEtBError); // Jet Correction histograms m_osHistSvc->fillTH1Flav(m_histDirJetCorr+"JetEtMinusHad", (etj-hadjetet)/GeV, m_event_weight, flav); m_osHistSvc->fillTH2Flav(m_histDirJetCorr+"JetEt_vs_Had", hadjetet/GeV, etj/GeV, m_event_weight,flav); m_osHistSvc->fillTH2Flav(m_histDirJetCorr+"JetEtMinusHad_vs_Had", hadjetet/GeV, (etj-hadjetet)/GeV, m_event_weight,flav); } } //Reweighted histograms - only fill if there is a matched jet if (nmatchedjets > 1 ) std::cout << " Matching more than 1 jet " << nmatchedjets << std::endl; if (nmatchedjets > 0 ) { m_osHistSvc->fillTH1Flav(m_histDirJets+"JetEtSysRwEtUp", etj/GeV, m_event_weight*rwetup, flav); m_osHistSvc->fillTH1Flav(m_histDirJets+"JetEtSysRwEtDo", etj/GeV, m_event_weight*rwetdo, flav); m_osHistSvc->fillTH1Flav(m_histDirJets+"JetEtaSysRwEtaUp", etaj, m_event_weight*rwetaup, flav); m_osHistSvc->fillTH1Flav(m_histDirJets+"JetEtaSysRwEtaDo", etaj, m_event_weight*rwetado, flav); //MC reweight systematics m_osHistSvc->fillTH1Flav(m_histDirJets+"JetWeightSysRwEtUp", jetweight, m_event_weight*rwetup, flav); m_osHistSvc->fillTH1Flav(m_histDirJets+"JetWeightSysRwEtDo", jetweight, m_event_weight*rwetdo, flav); m_osHistSvc->fillTH1Flav(m_histDirJets+"JetWeightSysRwEtaUp", jetweight, m_event_weight*rwetaup, flav); m_osHistSvc->fillTH1Flav(m_histDirJets+"JetWeightSysRwEtaDo", jetweight, m_event_weight*rwetado, flav); } // highest and second highest et jets if(etj>etj0) { etj2=etj1; etj1=etj0; etj0=etj; rwet0up=rwetup; rwet0do=rwetdo; rweta0up=rwetaup; rweta0do=rwetado; highestetjetmatchedtohad=nmatchedjets; } else if(etj>etj1) { etj2=etj1; etj1=etj; } else if(etj>etj2) { etj2=etj; } } // look at weights of individual algorithms double w_ip3d = jet->getFlavourTagWeight("IP3D"); double w_ip2d = jet->getFlavourTagWeight("IP2D"); m_osHistSvc->fillTH1Flav(m_histDirJets+"JetWeightIP3D", w_ip3d, m_event_weight, flav); m_osHistSvc->fillTH1Flav(m_histDirJets+"JetWeightIP2D", w_ip2d, m_event_weight, flav); // Look at basic track distributions for jet tracks and fill SV0 and SV1 as these depend on the info FillJetTrackInfo(jet, flav); // Look at input distributions for the taggers FillJetTagInfoHistos(jet, flav); } //High Jet probability events if (m_jetWeight0>4 && highestetjetmatchedtohad ) { FillHighProbHistos(m_jet0, m_jet1, m_flav0, m_flav1, rwet0up, rwet0do, rweta0up, rweta0do); } // number of jets in the event int ntotjet=m_nlightjet+m_ncjet+m_nbjet; // For higher multiplicities flag as the following // if at least one b jet flag as b // else if at least one c jet flag as c // else light jet int eventflav=0; if (ntotjet >= 1) { eventflav = (m_nbjet > 0) ? 5 : (m_ncjet > 0) ? 4 : 0 ; } m_osHistSvc->fillTH1Flav(m_histDirJets+"NJet", ntotjet, m_event_weight, eventflav); if(m_passCutsElectron) m_osHistSvc->fillTH1Flav(m_histDirJets+"ElecNJet", ntotjet, m_event_weight, eventflav); if(m_passCutsMuon) m_osHistSvc->fillTH1Flav(m_histDirJets+"MuonNJet", ntotjet, m_event_weight, eventflav); int nrecbjet=NBJets(); m_osHistSvc->fillTH1Flav(m_histDirJets+"NBJet", nrecbjet, m_event_weight, eventflav); if(m_passCutsElectron) m_osHistSvc->fillTH1Flav(m_histDirJets+"ElecNBJet",nrecbjet , m_event_weight, eventflav); if(m_passCutsMuon) m_osHistSvc->fillTH1Flav(m_histDirJets+"MuonNBJet", nrecbjet, m_event_weight, eventflav); if(ntotjet>1) { m_osHistSvc->fillTH1Flav(m_histDirJets+"NBJetDijet", nrecbjet, m_event_weight, eventflav); if(m_passCutsElectron) m_osHistSvc->fillTH1Flav(m_histDirJets+"ElecNBJetDijet",nrecbjet , m_event_weight, eventflav); if(m_passCutsMuon) m_osHistSvc->fillTH1Flav(m_histDirJets+"MuonNBJetDijet", nrecbjet, m_event_weight, eventflav); } if(ntotjet==1) { m_osHistSvc->fillTH1Flav(m_histDirJets+"NBJetNJet1", nrecbjet, m_event_weight, eventflav); if(m_passCutsElectron) m_osHistSvc->fillTH1Flav(m_histDirJets+"ElecNBJetNJet1",nrecbjet , m_event_weight, eventflav); if(m_passCutsMuon) m_osHistSvc->fillTH1Flav(m_histDirJets+"MuonNBJetNJet1", nrecbjet, m_event_weight, eventflav); } if(ntotjet==2) { m_osHistSvc->fillTH1Flav(m_histDirJets+"NBJetNJet2", nrecbjet, m_event_weight, eventflav); if(m_passCutsElectron) m_osHistSvc->fillTH1Flav(m_histDirJets+"ElecNBJetNJet2",nrecbjet , m_event_weight, eventflav); if(m_passCutsMuon) m_osHistSvc->fillTH1Flav(m_histDirJets+"MuonNBJetNJet2", nrecbjet, m_event_weight, eventflav); } if(ntotjet==3) { m_osHistSvc->fillTH1Flav(m_histDirJets+"NBJetNJet3", nrecbjet, m_event_weight, eventflav); if(m_passCutsElectron) m_osHistSvc->fillTH1Flav(m_histDirJets+"ElecNBJetNJet3",nrecbjet , m_event_weight, eventflav); if(m_passCutsMuon) m_osHistSvc->fillTH1Flav(m_histDirJets+"MuonNBJetNJet3", nrecbjet, m_event_weight, eventflav); } // Fill Ht histograms for various scenarios // filled m_osHistSvc->fillTH1Flav(m_histDirJets+"Ht",m_Ht/GeV, m_event_weight,eventflav); if(ntotjet>1) { m_osHistSvc->fillTH1Flav(m_histDirJets+"HtDijet",m_Ht/GeV, m_event_weight,eventflav); if (nrecbjet==1) m_osHistSvc->fillTH1Flav(m_histDirJets+"HtDijetNBJet1",m_Ht/GeV, m_event_weight,eventflav); if (nrecbjet==2) m_osHistSvc->fillTH1Flav(m_histDirJets+"HtDijetNBJet2",m_Ht/GeV, m_event_weight,eventflav); } if(ntotjet==2) { if (nrecbjet==1) m_osHistSvc->fillTH1Flav(m_histDirJets+"HtNJet2NBJet1",m_Ht/GeV, m_event_weight,eventflav); if (nrecbjet==2) m_osHistSvc->fillTH1Flav(m_histDirJets+"HtNJet2NBJet2",m_Ht/GeV, m_event_weight,eventflav); } if(ntotjet==3) { m_osHistSvc->fillTH1Flav(m_histDirJets+"HtNJet3NBJet1",m_Ht/GeV, m_event_weight,eventflav); if (nrecbjet==2) m_osHistSvc->fillTH1Flav(m_histDirJets+"HtNJet3NBJet2",m_Ht/GeV, m_event_weight,eventflav); } } void OSHbb::FillJetSystematicsHistos() { //Jet Enery Systematics // use the low et cut jets // Jet shift up m_selectJet->setSys(OSSelectorsEnums::eSysEnergyUp); JetCollection* jetsEUp=m_selectJet->selectedJets(); JetCollection::const_iterator jetUpItr (jetsEUp->begin()), jetUpItrE (jetsEUp->end()); // loop over the AOD jet container for (; jetUpItr != jetUpItrE; ++jetUpItr) { Jet* jet = *jetUpItr; float etj=m_selectJet->hlvCal(jet).perp(); // std::cout <<"etj up" <fillTH1Flav(m_histDirJets+"JetEtSysEUp", etj/GeV, m_event_weight, flav); // weight for combination of IP3D and SV1 double jetweight = jet->getFlavourTagWeight(); m_osHistSvc->fillTH1Flav(m_histDirJets+"JetWeightSysEUp", jetweight, m_event_weight, flav); } //Jet shift down m_selectJet->setSys(OSSelectorsEnums::eSysEnergyDo); JetCollection* jetsEDo=m_selectJet->selectedJets(); JetCollection::const_iterator jetDoItr (jetsEDo->begin()), jetDoItrE (jetsEDo->end()); // std::cout <<" nom " << m_jets->size() << " shift up " << jetsEUp->size() <<" shift down " << jetsEDo->size()<hlvCal(jet).perp(); //std::cout <<"etj do" <fillTH1Flav(m_histDirJets+"JetEtSysEDo", etj/GeV, m_event_weight, flav); // weight for combination of IP3D and SV1 double jetweight = jet->getFlavourTagWeight(); m_osHistSvc->fillTH1Flav(m_histDirJets+"JetWeightSysEDo", jetweight, m_event_weight, flav); } //Jet shift normal m_selectJet->setSys(OSSelectorsEnums::eSysNone); } /// void OSHbb::FillMETJetHistos() { msg(MSG::DEBUG) << "FillMETJetHistos()" << endreq; // number of jets in the event int numjets= m_selectJet->numSelected(); msg(MSG::DEBUG) << "Number of jets " << numjets << endreq; // we require at least 1 jet in order to have the chance of a b jet if (numjets <= 0) return; // Evaluate the event flavour for MC // if at least one b jet flag as b // else if at least one c jet flag as c // else light jet int ntotjet=m_nlightjet+m_ncjet+m_nbjet; int eventflav=0; if (ntotjet >= 1) { eventflav = (m_nbjet > 0) ? 5 : (m_ncjet > 0) ? 4 : 0 ; } int nrecbjet=NBJets(); msg(MSG::DEBUG) << "Number of b jets " << nrecbjet << endreq; // MET per number of b jets for all events with a jet if (nrecbjet==0) m_osHistSvc->fillTH1Flav(m_histDirMET+"METNBJet0",m_EtMiss/GeV, m_event_weight,m_event_flav); if (nrecbjet==1) m_osHistSvc->fillTH1Flav(m_histDirMET+"METNBJet1",m_EtMiss/GeV, m_event_weight,m_event_flav); if (nrecbjet==2) m_osHistSvc->fillTH1Flav(m_histDirMET+"METNBJet2",m_EtMiss/GeV, m_event_weight,m_event_flav); if (nrecbjet>=3) m_osHistSvc->fillTH1Flav(m_histDirMET+"METNBJet3Plus",m_EtMiss/GeV, m_event_weight,m_event_flav); // MET for events with 1 jet and 1 bjet if(ntotjet==1) { if (nrecbjet==1) m_osHistSvc->fillTH1Flav(m_histDirMET+"METNJet1NBJet1",m_EtMiss/GeV, m_event_weight,m_event_flav); else m_osHistSvc->fillTH1Flav(m_histDirMET+"METNJet1NBJet0",m_EtMiss/GeV, m_event_weight,m_event_flav); } // MET per number of b jets for all events with at least 2 jets if(ntotjet>1) { if (nrecbjet==0) m_osHistSvc->fillTH1Flav(m_histDirMET+"METDijetNBJet0",m_EtMiss/GeV, m_event_weight,m_event_flav); if (nrecbjet==1) m_osHistSvc->fillTH1Flav(m_histDirMET+"METDijetNBJet1",m_EtMiss/GeV, m_event_weight,m_event_flav); if (nrecbjet==2) m_osHistSvc->fillTH1Flav(m_histDirMET+"METDijetNBJet2",m_EtMiss/GeV, m_event_weight,m_event_flav); if (nrecbjet>=3) m_osHistSvc->fillTH1Flav(m_histDirMET+"METDijetNBJet3Plus",m_EtMiss/GeV, m_event_weight,m_event_flav); } // MET per number of b jets for all events with exactly 2 jets if(ntotjet==2) { if (nrecbjet==0) m_osHistSvc->fillTH1Flav(m_histDirMET+"METNJet2NBJet0",m_EtMiss/GeV, m_event_weight,m_event_flav); if (nrecbjet==1) m_osHistSvc->fillTH1Flav(m_histDirMET+"METNJet2NBJet1",m_EtMiss/GeV, m_event_weight,m_event_flav); if (nrecbjet==2) m_osHistSvc->fillTH1Flav(m_histDirMET+"METNJet2NBJet2",m_EtMiss/GeV, m_event_weight,m_event_flav); } // MET per number of b jets for all events with exactly 3 jets if(ntotjet==3) { if (nrecbjet==0) m_osHistSvc->fillTH1Flav(m_histDirMET+"METNJet3NBJet0",m_EtMiss/GeV, m_event_weight,m_event_flav); if (nrecbjet==1) m_osHistSvc->fillTH1Flav(m_histDirMET+"METNJet3NBJet1",m_EtMiss/GeV, m_event_weight,m_event_flav); if (nrecbjet==2) m_osHistSvc->fillTH1Flav(m_histDirMET+"METNJet3NBJet2",m_EtMiss/GeV, m_event_weight,m_event_flav); if (nrecbjet==3) m_osHistSvc->fillTH1Flav(m_histDirMET+"METNJet3NBJet3",m_EtMiss/GeV, m_event_weight,m_event_flav); } } void OSHbb::MakeJJMass(){ m_dibmass=0; m_dibmassopt=0; //We only look at di jet events from now on if(m_jet0==0 || m_jet1==0) return ; // Form di-b mass from scaled jets uising 2 highest btagged jets m_dibmass = m_dibmassscale*(m_selectJet->hlvCal(m_bjet0)+m_selectJet->hlvCal(m_bjet1)).m(); // Form di-b mass from scaled jets uising 2 highest pt jets m_dibmassopt = (m_selectJet->hlvCal(m_jet0)+m_selectJet->hlvCal(m_jet1)).m(); } void OSHbb::FillDijetHistos(){ //We only look at di jet events from now on if(m_jet0==0 || m_jet1==0) return ; // At least 2 jets FillllqqCutHisto(ellqqCutDiJet, m_event_weight, m_event_flav); FilllnubbCutHisto(elnubbCutDiJet, m_event_weight, m_event_flav); if (m_nbjet==1) FilllnubbCutHisto(elnubbCutDiJet1BJet, m_event_weight, m_event_flav); if (m_nbjet==2) FilllnubbCutHisto(elnubbCutDiJet2BJet, m_event_weight, m_event_flav); int numjets=m_jets->size(); if (numjets==2) { FilllnubbCutHisto(elnubbCut2Jet, m_event_weight, m_event_flav); if (m_nbjet==1) FilllnubbCutHisto(elnubbCut2Jet1BJet, m_event_weight, m_event_flav); if (m_nbjet==2) FilllnubbCutHisto(elnubbCut2Jet2BJet, m_event_weight, m_event_flav); } //We have at least 2 b jets (tagged or untagged) in the acceptance if(m_nbjet>=2) m_osHistSvc->fillTH1Flav(m_histDirCuts+"Cut", eCutBDijet, m_event_weight,m_event_flav); if(m_zbzlgen&&m_nbjet>=2) m_osHistSvc->fillTH1Flav(m_histDirCuts+"CutZb", eCutBDijet, m_event_weight,m_event_flav); //Without weight cuts msg(MSG::DEBUG) << "jet0 " << m_jet0 << endreq; msg(MSG::DEBUG) << "jet1 " << m_jet1 << endreq; m_event_counter[eECNumJetsGT1]+= m_event_weight; msg(MSG::DEBUG) << " Set NN" << endreq; SetNNVariables(); // look at trigger msg(MSG::DEBUG) << "Pass state L2_xe65 = " << m_trigDec->isPassed("L2_xe65") << endreq; msg(MSG::DEBUG) << "Pass state L2_e10_medium = " << m_trigDec->isPassed("L2_e10_medium") << endreq; bool passL2_xe65=m_trigDec->isPassed("L2_xe65"); m_osHistSvc->fillTH1Flav(m_histDirZ+"PtZDijetNoTrig", m_ptZ/GeV, m_event_weight); if (passL2_xe65) m_osHistSvc->fillTH1Flav(m_histDirZ+"PtZDijetTrig", m_ptZ/GeV, m_event_weight); // ScaleDiJetsUsingMETConstraint(); //We have difficulty here since jets can be different flavours. //If that is the case set to zero int flavuse = m_flav0==m_flav1 ? m_flav0 : 0; //Here all 3 jets are not HF but there is at least one other jet that is if(m_nbjet+m_ncjet>0&&m_flav0<3&&m_flav1<3&&m_flav2<3) m_osHistSvc->fillTH1Flav(m_histDirJets+"LostDijetMass", m_dibmassopt/GeV,m_event_weight,m_event_flav); if(m_ptZ>50*GeV) m_osHistSvc->fillTH1Flav(m_histDirJets+"DijetMassPtZ50", m_dibmassopt/GeV,m_event_weight,flavuse); if(m_nVtx3==1) { m_osHistSvc->fillTH1Flav(m_histDirJets+"DijetMassNVtx1", m_dibmassopt/GeV,m_event_weight,flavuse); } else if(m_nVtx3>1) { m_osHistSvc->fillTH1Flav(m_histDirJets+"DijetMassNVtxGt1", m_dibmassopt/GeV,m_event_weight,flavuse); } m_osHistSvc->fillTH1Flav(m_histDirZ+"PtZDijet", m_ptZ/GeV, m_event_weight); return ; } double OSHbb::ZZHiggsMass(Jet* jet0, Jet* jet1, int type, bool scale) { //For H->ZZ search make Higgs mass double HiggsZZMass=0; //return here if we have no jets and Z->bb selected if(m_VectorBoson2Op==eZbb && (jet0==0 || jet1==0)) return HiggsZZMass; //Make the higgs mass for H->ZZ decay static const double M_Z=91.188*GeV; //Scale factor for b quarks to make sure di-b invariant mass is equal to Z mass HepLorentzVector hlvJ0ZC=m_selectJet->hlvCal(jet0); HepLorentzVector hlvJ1ZC=m_selectJet->hlvCal(jet1); if (scale && jet0>0 && jet1>0) { ////CBG ScaleDijetsUsingZConstraint(hlvJ0ZC, hlvJ1ZC); } //Simple scale factor for leptons to make sure di-lepton invariant mass is equal to Z mass // double scalel=M_Z/m_mZ; double scalel=1.0; //Doesn't really help so set it to 1 for now if( m_VectorBosonOp==eZ) { HepLorentzVector lep0, lep1; int leppdg; if(m_electrons->size()>1) { lep0=m_selectElec->hlvCal(m_electrons->at(0))*scalel; lep1=m_selectElec->hlvCal(m_electrons->at(1))*scalel; leppdg = 11; } else if(m_muons->size()>1) { lep0=m_selectMuon->hlvCal(m_muons->at(0))*scalel; lep1=m_selectMuon->hlvCal(m_muons->at(1))*scalel; leppdg = 13; } if(m_VectorBoson2Op==eZnu){ //We have llnunu final state so we must make some assumptions about the Z->nunu. We make a particle with Pt=Etmiss, eta=0 and mass=M_Z Hep3Vector vec3nn(m_EtMiss_x,m_EtMiss_y,0); HepLorentzVector znn(0.0,0.0,0.0,0.0); znn.setVectM(vec3nn,M_Z); //We do same for Z->ll so we end up with the transverse mass HepLorentzVector zll(0.0,0.0,0.0,0.0); Hep3Vector vec3ll(lep0.px()+lep1.px(),lep0.py()+lep1.py(),0); zll.setVectM(vec3ll,M_Z); HiggsZZMass=(zll+znn).m(); } else if (m_VectorBoson2Op==eZbb) { //HiggsZZMass=(hlvJ0ZC+hlvJ1ZC+lep0+lep1).m(); // Only scale if both on-shell and for b jets for now if (!((m_zmasslow < m_mZ && m_mZ < m_zmasshigh) && ((type == 0 && m_zbmasslow < m_dibmassopt && m_dibmassopt < m_zbmasshigh) || (type == 1 && m_zbmasslow < m_dibmass && m_dibmass < m_zbmasshigh)) && (jet0->e() > 0 && jet1->e() > 0))) { HiggsZZMass=(hlvJ0ZC+hlvJ1ZC+lep0+lep1).m(); } else if (type==1 && m_bjets->size()==2) { m_kinFitter->reset(); m_kinFitter->setZll(lep0, lep1, leppdg, false); m_kinFitter->setZjj(hlvJ0ZC, hlvJ1ZC, false); std::vector otherJets; if (m_jets->size() > 2) { JetCollection::const_iterator jetItr (m_bjets->begin()+2), jetEnd (m_bjets->end()); for (; jetItr != jetEnd; ++jetItr) { otherJets.push_back(m_selectJet->hlvCal(*jetItr)); } } m_kinFitter->setMETHZZllqq(otherJets, m_EtMiss_x, m_EtMiss_y); m_kinFitter->fit(); m_kinFitter->print(); HepLorentzVector* lep0Cal = m_kinFitter->getZllLep1(); HepLorentzVector* lep1Cal = m_kinFitter->getZllLep2(); HepLorentzVector* jet0Cal = m_kinFitter->getZjjJet1(); HepLorentzVector* jet1Cal = m_kinFitter->getZjjJet2(); HiggsZZMass=(*lep0Cal+*lep1Cal+*jet0Cal+*jet1Cal).m(); } } else { std::cout <<" Error in VectorBoson2Op. We don't have this option implemented yet" << m_VectorBoson2Op<nunu. We make a particle with Pt=Etmiss, eta=0 and mass=M_Z Hep3Vector vec3nn(m_EtMiss_x,m_EtMiss_y,0); HepLorentzVector znn(0.0,0.0,0.0,0.0); znn.setVectM(vec3nn,M_Z); //We take the scaled jets HiggsZZMass=(hlvJ0ZC+hlvJ1ZC+znn).m(); } return HiggsZZMass; } void OSHbb::FillHiggsZZHistos() { //We have difficulty here since jets can be different flavours. //If that is the case set to zero int flavuse = m_flav0==m_flav1 ? m_flav0 : 0; //For Z events add in the Z also for H->ZZ search. Only make this quantity if m_dibmassopt is consistent with Z mass // fill a histogram and cut on the angular separation of the jets m_osHistSvc->fillTH1Flav(m_histDirJets+"DPhiJJBeforeMCut", m_JJDeltaPhiOPt,m_event_weight,flavuse); // fill a histogram and cut on the angular separation of the leptons m_osHistSvc->fillTH1Flav(m_histDirZ+"DPhiLLBeforeMCut", m_LLDeltaPhi,m_event_weight,flavuse); if(ZBJetMassCutOPt()) { FillllqqCutHisto(ellqqCutMqq, m_event_weight, m_event_flav); if(m_JJDeltaPhiOPthlvCal(m_jet0).perp()>50*GeV&&m_selectJet->hlvCal(m_jet1).perp()>50*GeV) FillllqqCutHisto(ellqqCutPt50, m_event_weight, m_event_flav); } } // fill a histogram and cut on the angular separation of the jets m_osHistSvc->fillTH1Flav(m_histDirJets+"DPhiJJ", m_JJDeltaPhiOPt,m_event_weight,flavuse); // fill a histogram and cut on the angular separation of the leptons m_osHistSvc->fillTH1Flav(m_histDirZ+"DPhiLL", m_LLDeltaPhi,m_event_weight,flavuse); if(m_JJDeltaPhiOPtfillTH1Flav(m_histDirZ+"DPhiLLDPhiJJCut", m_LLDeltaPhi,m_event_weight,flavuse); if(m_ptZ>m_PtZCut) m_osHistSvc->fillTH1Flav(m_histDirJets+"DijetZMassPtZ", m_HiggsZZMassOPt/GeV,m_event_weight,flavuse); } else if (m_dibmassopt>m_zbmasslowsb&&m_dibmassoptfillTH1Flav(m_histDirJets+"DijetZMassSBLow", m_HiggsZZMassOPt/GeV,m_event_weight,flavuse); if(m_JJDeltaPhiOPtfillTH1Flav(m_histDirJets+"DijetZMassSBLowDPhiCuts", m_HiggsZZMassOPtNoScale/GeV,m_event_weight,flavuse); if(m_selectJet->hlvCal(m_jet0).perp()>50*GeV&&m_selectJet->hlvCal(m_jet1).perp()>50*GeV) { m_osHistSvc->fillTH1Flav(m_histDirJets+"DijetZMassSBLowPtJ50", m_HiggsZZMassOPtNoScale/GeV,m_event_weight,flavuse); if(m_JJDeltaPhiOPtfillTH1Flav(m_histDirJets+"DijetZMassSBLowPtJ50DPhiCuts", m_HiggsZZMassOPtNoScale/GeV,m_event_weight,flavuse); } } else if (m_dibmassopt>m_zbmasshigh&&m_dibmassopt< m_zbmasshighsb) { m_osHistSvc->fillTH1Flav(m_histDirJets+"DijetZMassSBHigh", m_HiggsZZMassOPtNoScale/GeV,m_event_weight,flavuse); if(m_JJDeltaPhiOPtfillTH1Flav(m_histDirJets+"DijetZMassSBHighDPhiCuts", m_HiggsZZMassOPtNoScale/GeV,m_event_weight,flavuse); if(m_selectJet->hlvCal(m_jet0).perp()>50*GeV&&m_selectJet->hlvCal(m_jet1).perp()>50*GeV) { m_osHistSvc->fillTH1Flav(m_histDirJets+"DijetZMassSBHighPtJ50", m_HiggsZZMassOPtNoScale/GeV,m_event_weight,flavuse); if(m_JJDeltaPhiOPtfillTH1Flav(m_histDirJets+"DijetZMassSBHighPtJ50DPhiCuts", m_HiggsZZMassOPtNoScale/GeV,m_event_weight,flavuse); } } } void OSHbb::FillDijetHistosAfterBTag() { int flavuse = m_flav0==m_flav1 ? m_flav0 : 0; m_osHistSvc->fillTH1Flav(m_histDirJets+"DijetMassWCNoCal", m_dibmass/GeV/m_dibmassscale ,m_event_weight,flavuse); //Fill light-light, light b, light c histos int flavuseL = TMath::Max(m_flav0,m_flav1); if( m_flav0==0&&m_flav1==0) m_osHistSvc->fillTH1Flav(m_histDirJets+"DijetMassWCL", m_dibmass/GeV,m_event_weight,flavuseL); // number of jets - event level weight int ntotjet=m_nlightjet+m_ncjet+m_nbjet; int eventflav=0; if (ntotjet >= 1) { eventflav = (m_nbjet > 0) ? 5 : (m_ncjet > 0) ? 4 : 0 ; } if(m_ptZ>50*GeV) m_osHistSvc->fillTH1Flav(m_histDirJets+"DijetMassWCPtZ50", m_dibmass/GeV,m_event_weight,flavuse); if(m_ptZ>100*GeV) m_osHistSvc->fillTH1Flav(m_histDirJets+"DijetMassWCPtZ100", m_dibmass/GeV,m_event_weight,flavuse); m_osHistSvc->fillTH1Flav(m_histDirZ+"PtZDijetWC", m_ptZ/GeV, m_event_weight); } void OSHbb::FillHiggsZZHistosAfterBTag() { int flavuse = m_flav0==m_flav1 ? m_flav0 : 0; //Only fill for masses around the Z mass if(m_dibmass>m_zbmasslowsb&&m_dibmassfillTH1Flav(m_histDirCuts+"Cut", eCutBZMass, m_event_weight,m_event_flav); if(m_zbzlgen) m_osHistSvc->fillTH1Flav(m_histDirCuts+"CutZb", eCutBZMass, m_event_weight,m_event_flav); int flavuseL = TMath::Max(m_flav0,m_flav1); if( m_flav0==0&&m_flav1==0) { m_osHistSvc->fillTH1Flav(m_histDirJets+"DijetZMassWCL", m_HiggsZZMass/GeV,m_event_weight,flavuseL); if(m_JJDeltaPhifillTH1Flav(m_histDirJets+"DijetZMassDPhiCutsWCL", m_HiggsZZMass/GeV,m_event_weight,flavuseL); } m_osHistSvc->fillTH1Flav(m_histDirJets+"DPhiJJWC", m_JJDeltaPhi,m_event_weight,flavuse); // fill a histogram and cut on the angular separation of the leptons m_osHistSvc->fillTH1Flav(m_histDirZ+"DPhiLLWC", m_LLDeltaPhi,m_event_weight,flavuse); if(m_JJDeltaPhifillTH1Flav(m_histDirZ+"DPhiLLDPhiJJCutWC", m_LLDeltaPhi,m_event_weight,flavuse); if(m_JJDeltaPhifillTH1Flav(m_histDirCuts+"TriggersElec", i ,m_event_weight,flavuse); else if(m_trigDec->isPassed(m_triggers[i])) m_osHistSvc->fillTH1Flav(m_histDirCuts+"TriggersElec", i ,m_event_weight,flavuse); } if (m_passCutsMuon) { if (i==0) m_osHistSvc->fillTH1Flav(m_histDirCuts+"TriggersMuon", i ,m_event_weight,flavuse); else if(m_trigDec->isPassed(m_triggers[i])) m_osHistSvc->fillTH1Flav(m_histDirCuts+"TriggersMuon", i ,m_event_weight,flavuse); } } } if( m_ptZ>m_PtZCut) m_osHistSvc->fillTH1Flav(m_histDirJets+"DijetZMassPtZWC", m_HiggsZZMass/GeV,m_event_weight,flavuse); } else if(m_dibmass>m_zbmasslowsb&&m_dibmassfillTH1Flav(m_histDirJets+"DijetZMassSBLowWC", m_HiggsZZMassNoScale/GeV,m_event_weight,flavuse); if(m_JJDeltaPhifillTH1Flav(m_histDirJets+"DijetZMassSBLowDPhiCutsWC", m_HiggsZZMassNoScale/GeV,m_event_weight,flavuse); } else if(m_dibmass>m_zbmasshigh&&m_dibmassfillTH1Flav(m_histDirJets+"DijetZMassSBHighWC", m_HiggsZZMassNoScale/GeV,m_event_weight,flavuse); if(m_JJDeltaPhifillTH1Flav(m_histDirJets+"DijetZMassSBHighDPhiCutsWC", m_HiggsZZMassNoScale/GeV,m_event_weight,flavuse); } } return; } float OSHbb::coshel(TLorentzVector particle, TLorentzVector parent) { //http://www.slac.stanford.edu/BFROOT/www/doc/workbook/analysis/analysis.html#helicity //A useful quantity in many analyses is the helicity angle. In the reaction // Y -> X -> a + b, the helicity angle of particle a is the angle measured //in the rest frame of the decaying parent particle, X, between the //direction of the decay daughter a and the direction of the grandparent //particle Y. (Source: BAD 529 (V5) p.5-6) TVector3 boosttoparent = -(parent.BoostVector()); particle.Boost(boosttoparent); parent.Boost(boosttoparent); return cos(parent.Angle(particle.Vect())); } void OSHbb::BookHighProbHistos() { // jet BHist histograms at high jet probability m_histDirJetHP = fTHistStreamName+"JetHighProb/"; m_osHistSvc->bookFlavHist(m_histDirJetHP+"NJetHighProb","N^{jet}",5,0.,10); m_osHistSvc->bookFlavHist(m_histDirJetHP+"Jet0EtHighProb","E_{T}^{jet}",50,0.,250.); m_osHistSvc->bookFlavHist(m_histDirJetHP+"Jet1EtHighProb","E_{T}^{jet}",50,0.,250.); m_osHistSvc->bookFlavHist(m_histDirJetHP+"Jet0EtaHighProb","{#eta}^{jet}",70,-3.5,3.5); m_osHistSvc->bookFlavHist(m_histDirJetHP+"Jet1EtaHighProb","{#eta}^{jet}",70,-3.5,3.5); // systematics m_osHistSvc->bookFlavHist(m_histDirJetHP+"Jet0EtHighProbSysRwEtUp","E_{T}^{jet}",50,0.,250.); m_osHistSvc->bookFlavHist(m_histDirJetHP+"Jet0EtHighProbSysRwEtDo","E_{T}^{jet}",50,0.,250.); m_osHistSvc->bookFlavHist(m_histDirJetHP+"Jet0EtaHighProbSysRwEtaUp","{#eta}^{jet}",70,-3.5,3.5); m_osHistSvc->bookFlavHist(m_histDirJetHP+"Jet0EtaHighProbSysRwEtaDo","{#eta}^{jet}",70,-3.5,3.5); return; } void OSHbb::FillHighProbHistos(Jet* jet0, Jet* jet1, int flav0, int flav1, double rwetup, double rwetdo, double rwetaup, double rwetado) { msg(MSG::DEBUG) << "FillHighProbHistos" << endreq; if (jet0 == 0) return; double etj0 = m_selectJet->hlvCal(jet0).perp(); double etaj0 = m_selectJet->hlvCal(jet0).eta(); // first jet m_osHistSvc->fillTH1Flav(m_histDirJetHP+"Jet0EtHighProb", etj0/GeV, m_event_weight, flav0); m_osHistSvc->fillTH1Flav(m_histDirJetHP+"Jet0EtaHighProb", etaj0, m_event_weight, flav0); // syst m_osHistSvc->fillTH1Flav(m_histDirJetHP+"Jet0EtHighProbSysRwEtUp", etj0/GeV, m_event_weight*rwetup, flav0); m_osHistSvc->fillTH1Flav(m_histDirJetHP+"Jet0EtHighProbSysRwEtDo", etj0/GeV, m_event_weight*rwetdo, flav0); m_osHistSvc->fillTH1Flav(m_histDirJetHP+"Jet0EtaHighProbSysRwEtaUp", etaj0, m_event_weight*rwetaup, flav0); m_osHistSvc->fillTH1Flav(m_histDirJetHP+"Jet0EtaHighProbSysRwEtaDo", etaj0, m_event_weight*rwetado, flav0); // second jet if (jet1==0 ) return; double etj1 = m_selectJet->hlvCal(jet1).perp(); double etaj1 = m_selectJet->hlvCal(jet1).eta(); m_osHistSvc->fillTH1Flav(m_histDirJetHP+"Jet1EtHighProb", etj1/GeV, m_event_weight, flav1); m_osHistSvc->fillTH1Flav(m_histDirJetHP+"Jet1EtaHighProb", etaj1, m_event_weight, flav1); return; } void OSHbb::FillGenHistos() { msg(MSG::DEBUG) << "FillGenHistos" << endreq; int nlightjet(0); int ncjet(0); int nbjet(0); // Generator hadron level jets JetCollection* hadjets = m_selectJet->selectedHadJets(); if(hadjets) { // When running on data hadJets may not exist JetCollection::const_iterator jetItrH (hadjets->begin()), jetItrHE (hadjets->end()); for (; jetItrH != jetItrHE; ++jetItrH) { Jet* hadjet(*jetItrH); double hadjetet = hadjet->et(); double hadjeteta = hadjet->eta(); int flav = GetFlavourHad(hadjet); m_osHistSvc->fillTH1Flav(m_histDirTruth+"McHadBJets", 0.5, m_event_weight, flav); m_osHistSvc->fillTH1Flav(m_histDirTruth+"JetEtGen", hadjetet/GeV, m_event_weight, flav); m_osHistSvc->fillTH1Flav(m_histDirTruth+"JetEtaGen", hadjeteta, m_event_weight, flav); // number of jets if (flav == 5) nbjet++; else if (flav == 4) ncjet++; else nlightjet++; //systematics double rwetaup=1.0; double rwetado=1.0; double rwetup=1.0; double rwetdo=1.0; rwetaup= pow(fabs(hadjeteta)+m_RwEtaBCen,m_RwEtaBError); rwetado= pow(fabs(hadjeteta)+m_RwEtaBCen,-m_RwEtaBError); rwetup= pow((hadjetet/m_RwEtBCen),m_RwEtBError); rwetdo= pow((hadjetet/m_RwEtBCen),-m_RwEtBError); //Reweighted histograms m_osHistSvc->fillTH1Flav(m_histDirTruth+"McHadBJetsSysRwEtUp", 0.5, m_event_weight*rwetup, flav); m_osHistSvc->fillTH1Flav(m_histDirTruth+"McHadBJetsSysRwEtDo", 0.5, m_event_weight*rwetdo, flav); m_osHistSvc->fillTH1Flav(m_histDirTruth+"McHadBJetsSysRwEtaUp", 0.5, m_event_weight*rwetaup, flav); m_osHistSvc->fillTH1Flav(m_histDirTruth+"McHadBJetsSysRwEtaDo", 0.5, m_event_weight*rwetado, flav); } } // number of jets in the event int ntotjet=nlightjet+ncjet+nbjet; // For higher multiplicities flag as the following // if at least one b jet flag as b // else if at least one c jet flag as c // else light jet int eventflav=0; if (ntotjet >= 1) eventflav = (nbjet > 0) ? 5 : (ncjet > 0) ? 4 : 0 ; m_osHistSvc->fillTH1Flav(m_histDirTruth+"NHadJet", ntotjet, m_event_weight, eventflav); // For 1 jet fill with flavour if (ntotjet==1) m_osHistSvc->fillTH1Flav(m_histDirTruth+"Flav1HadJet", ntotjet, m_event_weight, eventflav); // For 2 jets fill as ll lc cc lb bb cb // 0 1 2 3 4 5 int dijetflag=-1; if (ntotjet==2) { if (nlightjet==2) dijetflag=0; else if (nlightjet==1 && ncjet==1) dijetflag=1; else if (ncjet==2) dijetflag=2; else if (nlightjet==1 && nbjet==1) dijetflag=3; else if (nbjet==2) dijetflag=4; else if (ncjet==1 && nbjet==1) dijetflag=5; m_osHistSvc->fillTH1Flav(m_histDirTruth+"Flav2HadJet", dijetflag, m_event_weight,0); } // Generator parton level jets if (m_doPartonJets) { int nlightpartjet(0), ncpartjet(0), nbpartjet(0); JetCollection* partjets = m_selectJet->selectedPartJets(); if(partjets) { // When running on data partJets may not exist JetCollection::const_iterator partJetItr(partjets->begin()), partJetEnd(partjets->end()); for (; partJetItr != partJetEnd; ++partJetItr) { Jet* partjet(*partJetItr); double partjetet = partjet->et(); double partjeteta = partjet->eta(); int flav = GetFlavourPart(partjet); m_osHistSvc->fillTH1Flav(m_histDirTruth+"McPartBJets", 0.5, m_event_weight, flav); m_osHistSvc->fillTH1Flav(m_histDirTruth+"JetEtGenPart", partjetet/GeV, m_event_weight, flav); m_osHistSvc->fillTH1Flav(m_histDirTruth+"JetEtaGenPart", partjeteta, m_event_weight, flav); // number of jets if (flav == 5) nbpartjet++; else if (flav == 4) ncpartjet++; else nlightpartjet++; } } // number of jets in the event int ntotpartjet=nlightpartjet+ncpartjet+nbpartjet; // For higher multiplicities flag as the following // if at least one b jet flag as b // else if at least one c jet flag as c // else light jet int parteventflav=0; if (ntotjet >= 1) parteventflav = (nbpartjet > 0) ? 5 : (ncpartjet > 0) ? 4 : 0 ; m_osHistSvc->fillTH1Flav(m_histDirTruth+"NPartJet", ntotpartjet, m_event_weight, parteventflav); // For 1 jet fill with flavour if (ntotpartjet==1) m_osHistSvc->fillTH1Flav(m_histDirTruth+"Flav1PartJet", ntotpartjet, m_event_weight, parteventflav); // For 2 jets fill as ll lc cc lb bb cb // 0 1 2 3 4 5 int dipartjetflag(-1); if (ntotpartjet==2) { if (nlightpartjet==2) dipartjetflag=0; else if (nlightpartjet==1 && ncpartjet==1) dipartjetflag=1; else if (ncpartjet==2) dipartjetflag=2; else if (nlightpartjet==1 && nbpartjet==1) dipartjetflag=3; else if (nbpartjet==2) dipartjetflag=4; else if (ncpartjet==1 && nbpartjet==1) dipartjetflag=5; m_osHistSvc->fillTH1Flav(m_histDirTruth+"Flav2PartJet", dipartjetflag, m_event_weight,0); } } return; } void OSHbb::FillJetTrackInfo(Jet* jet, int flav) { // general info on tracks associated to the jet msg(MSG::DEBUG) << "FillJetTrackInfo" << endreq; flav=-1; // not used const Analysis::TrackAssociation *tpc = jet->getAssociation("Tracks"); if (!tpc) return; // int ntrk_j = tpc->nTracks(); // loop over tracks in a jet std::vector* trackVector = tpc->tracks(); std::vector::const_iterator trkItr(trackVector->begin()), trkEnd(trackVector->end()); for(; trkItr != trkEnd; ++trkItr) { const Rec::TrackParticle* aTemp = *trkItr; if (!aTemp) continue; const Trk::TrackSummary* summary = aTemp->trackSummary(); if (!summary) continue; // float nBLayerHits = summary->get(Trk::numberOfBLayerHits); // float nPixelHits = summary->get(Trk::numberOfPixelHits); // float nSCTHits =summary->get(Trk::numberOfSCTHits); // std::cout << " hits " << nBLayerHits << " " << nPixelHits << " " << nSCTHits << std::endl; } return; } void OSHbb::BookJetTagInfoHistos() { // jet BHist histograms for tagger input quantities m_histDirJetTag = fTHistStreamName+"JetTagInfo/"; // d0(r-phi) significance - IP3D m_osHistSvc->bookFlavHist(m_histDirJetTag+"JetTracksd0Sig","",100,-10,10); // z0 significance - IP3D m_osHistSvc->bookFlavHist(m_histDirJetTag+"JetTracksz0Sig","",100,-10,10); // Fraction of jet energy in SV - SV1 m_osHistSvc->bookFlavHist(m_histDirJetTag+"JetEFracSV1","",50,0.,1); // Mass of the SV - SV1 m_osHistSvc->bookFlavHist(m_histDirJetTag+"JetMassSV1","",50,0.,10); // Number of 2 track vertices in the event - SV1 m_osHistSvc->bookFlavHist(m_histDirJetTag+"JetN2TrackSV1","",20,0.,20); // Same for SV0 m_osHistSvc->bookFlavHist(m_histDirJetTag+"JetEFracSV0","",50,0.,1); m_osHistSvc->bookFlavHist(m_histDirJetTag+"JetMassSV0","",50,0.,10); m_osHistSvc->bookFlavHist(m_histDirJetTag+"JetMassSV0Cut","",50,0.,10); m_osHistSvc->bookFlavHist(m_histDirJetTag+"JetN2TrackSV0","",20,0.,20); // SV0 m_osHistSvc->bookFlavHist(m_histDirJetTag+"JetMassSV0CutNJet1","",50,0.,10); m_osHistSvc->bookFlavHist(m_histDirJetTag+"JetMassSV0CutNJet2","",50,0.,10); return; } void OSHbb::FillJetTagInfoHistos(Jet* jet, int flav) { // information on the input to the jet tagger msg(MSG::DEBUG) << "FillJetTagInfoHistos" << endreq; double w_sv1 = jet->getFlavourTagWeight("SV1"); double w_sv0 = jet->getFlavourTagWeight("SV0"); int numjets= m_selectJet->numSelected(); // --- loop on tags & get detailed information for each tag std::vector infoVector = jet->jetTagInfoVector(); for(uint iInfo = 0; iInfo < infoVector.size(); iInfo++) { bool is3D = (infoVector[iInfo]->infoType()=="IP3D"); bool isIPlus = (infoVector[iInfo]->infoType()=="IPInfoPlus"); bool isSV1Plus = (infoVector[iInfo]->infoType()=="SVInfoPlus"); bool isSV0Plus = (infoVector[iInfo]->infoType()=="SV0InfoPlus"); // bool is2D = (infoVector[iInfo]->infoType()=="IP2D"); // bool isS1 = (infoVector[iInfo]->infoType()=="SV1"); // bool isS2 = (infoVector[iInfo]->infoType()=="SV2"); // bool isSoftElectron = (infoVector[iInfo]->infoType()=="SoftElectronTag"); // basic IP3D information: if (is3D) { //const Analysis::IPInfoBase* infob = dynamic_cast(infoVector[iInfo]); // if(infob) { // int ntrk = infob->nbTracks(); // number of tracks used for tagging in the jet // double pb = infob->tagLikelihood().at(0); // b likelihood // double pu = infob->tagLikelihood()[1]; // u likelihood // << " Pb= " << pb << " Pu= " << pu << std::endl; // } } // if(isIPlus) { const Analysis::IPInfoPlus* infop = dynamic_cast(infoVector[iInfo]); if (infop) { int ntrk = infop->numTrackInfo(); for(int itinf = 0; itinf < ntrk; itinf++) { Analysis::IPTrackInfo trackInfo = infop->getTrackInfo(itinf); // impact parameters w.r.t. primary vertex: // double d0val = trackInfo.d0Value(); double d0sig = trackInfo.d0Significance(); //double z0val = trackInfo.z0Value(); double z0sig = trackInfo.z0Significance(); m_osHistSvc->fillTH1Flav(m_histDirJetTag+"JetTracksd0Sig", d0sig, m_event_weight, flav); m_osHistSvc->fillTH1Flav(m_histDirJetTag+"JetTracksz0Sig", z0sig, m_event_weight, flav); // individual contribution of this track to the b-tagging weights: // double wi2D = trackInfo.trackWeight2D(); // IP2D // double wi3D = trackInfo.trackWeight3D(); // IP3D // double piJP = trackInfo.trackProbJP(); // JetProb // track quality information: // Analysis::TrackGrade grd = trackInfo.trackGrade(); // Good, Shared, // bool vzero = trackInfo.isFromV0(); // track from V0, interaction in identified material, ... // // pointer to actual TrackParticle: // const Rec::TrackParticle* trk = trackInfo.track(); // trk->eta(),rk->phi() } } } if(isSV1Plus) { const Analysis::SVInfoPlus* info = dynamic_cast(infoVector[iInfo]); if(info) { // int ntrk = info->getNGTrackInSvx(); // number of tracks in vertex double mass = info->getMass(); // mass of secondary vertex int n2t = info->getN2T(); // number of two-track vertices bool svok = (mass>0. && n2t>0 ); double efrc = info->getEnergyFraction(); // energy fraction svx/jet if(svok) { m_osHistSvc->fillTH1Flav(m_histDirJets+"JetWeightSV1", w_sv1, m_event_weight, flav); m_osHistSvc->fillTH1Flav(m_histDirJetTag+"JetEFracSV1", efrc, m_event_weight, flav); m_osHistSvc->fillTH1Flav(m_histDirJetTag+"JetMassSV1", mass, m_event_weight, flav); m_osHistSvc->fillTH1Flav(m_histDirJetTag+"JetN2TrackSV1", n2t, m_event_weight, flav); if (numjets==1) m_osHistSvc->fillTH1Flav(m_histDirJets+"JetWeightSV1NJet1", w_sv1, m_event_weight, flav); if (numjets==2) m_osHistSvc->fillTH1Flav(m_histDirJets+"JetWeightSV1NJet2", w_sv1, m_event_weight, flav); } } } if(isSV0Plus) { float sv0cut= 5.85; const Analysis::SVInfoPlus* info = dynamic_cast(infoVector[iInfo]); if(info) { // int ntrk = info->getNGTrackInSvx(); // number of tracks in vertex double mass = info->getMass(); // mass of secondary vertex int n2t = info->getN2T(); // number of two-track vertices bool svok = (mass>0. && n2t>0 ); double efrc = info->getEnergyFraction(); // energy fraction svx/jet if(svok) { m_osHistSvc->fillTH1Flav(m_histDirJets+"JetWeightSV0", w_sv0, m_event_weight, flav); m_osHistSvc->fillTH1Flav(m_histDirJetTag+"JetEFracSV0", efrc, m_event_weight, flav); m_osHistSvc->fillTH1Flav(m_histDirJetTag+"JetMassSV0", mass, m_event_weight, flav); if (numjets==1) m_osHistSvc->fillTH1Flav(m_histDirJets+"JetWeightSV0NJet1", w_sv0, m_event_weight, flav); if (numjets==2) m_osHistSvc->fillTH1Flav(m_histDirJets+"JetWeightSV0NJet2", w_sv0, m_event_weight, flav); if(w_sv0>sv0cut) { m_osHistSvc->fillTH1Flav(m_histDirJetTag+"JetMassSV0Cut", mass, m_event_weight, flav); if (numjets==1) m_osHistSvc->fillTH1Flav(m_histDirJetTag+"JetMassSV0CutNJet1", mass, m_event_weight, flav); if (numjets==2) m_osHistSvc->fillTH1Flav(m_histDirJetTag+"JetMassSV0CutNJet2", mass, m_event_weight, flav); } m_osHistSvc->fillTH1Flav(m_histDirJetTag+"JetN2TrackSV0", n2t, m_event_weight, flav); } } } } return; } void OSHbb::FillIsolHistos() { // Isolation hists: track and calo isol for muons/elecs and distances to jets. // Only do for events passing the leptons cuts containing at least 2 jets. // In the case of the Z analysis also do after the Z mass cut if (!m_doIsol || m_jets->size()<2) return; Muon* muon0 = static_cast(m_mu0); Muon* muon1 = static_cast(m_mu1); Electron* elec0 = static_cast(m_e0); Electron* elec1 = static_cast(m_e1); if (m_passCutsMuon && muon0) { float muon0CaloIsol = m_selectMuon->caloIsol(muon0, MuonParameters::etcone20, false); m_osHistSvc->fillTH1Flav(m_histDirIsol+"MuonCaloIsol",muon0CaloIsol/1000., m_event_weight, m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirIsol+"MuonCaloIsolRel",muon0CaloIsol/muon0->pt(), m_event_weight, m_event_flav); if (m_passCutZMass) { m_osHistSvc->fillTH1Flav(m_histDirIsol+"MuonCaloIsolMZ",muon0CaloIsol/1000., m_event_weight, m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirIsol+"MuonCaloIsolRelMZ",muon0CaloIsol/muon0->pt(), m_event_weight, m_event_flav); } if (muon0->hasInDetTrackParticle()) { float muon0TrackIsol = m_selectMuon->trackIsol(muon0, 0.2, false); m_osHistSvc->fillTH1Flav(m_histDirIsol+"MuonTrackIsol",muon0TrackIsol/1000., m_event_weight, m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirIsol+"MuonTrackIsolRel",muon0TrackIsol/muon0->inDetTrackParticle()->pt(), m_event_weight, m_event_flav); if (m_passCutZMass) { m_osHistSvc->fillTH1Flav(m_histDirIsol+"MuonTrackIsolMZ",muon0TrackIsol/1000., m_event_weight, m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirIsol+"MuonTrackIsolRelMZ",muon0TrackIsol/muon0->inDetTrackParticle()->pt(), m_event_weight, m_event_flav); } } if (muon1) { float muon1CaloIsol = m_selectMuon->caloIsol(muon1, MuonParameters::etcone20, false); m_osHistSvc->fillTH1Flav(m_histDirIsol+"MuonCaloIsol",muon1CaloIsol/1000., m_event_weight, m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirIsol+"MuonCaloIsolRel",muon1CaloIsol/muon1->pt(), m_event_weight, m_event_flav); if (m_passCutZMass) { m_osHistSvc->fillTH1Flav(m_histDirIsol+"MuonCaloIsolMZ",muon1CaloIsol/1000., m_event_weight, m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirIsol+"MuonCaloIsolRelMZ",muon1CaloIsol/muon1->pt(), m_event_weight, m_event_flav); } if (muon1->hasInDetTrackParticle()) { float muon1TrackIsol = m_selectMuon->trackIsol(muon1, 0.2, false); m_osHistSvc->fillTH1Flav(m_histDirIsol+"MuonTrackIsol",muon1TrackIsol/1000., m_event_weight, m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirIsol+"MuonTrackIsolRel",muon1TrackIsol/muon1->inDetTrackParticle()->pt(), m_event_weight, m_event_flav); if (m_passCutZMass) { m_osHistSvc->fillTH1Flav(m_histDirIsol+"MuonTrackIsolMZ",muon1TrackIsol/1000., m_event_weight, m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirIsol+"MuonTrackIsolRelMZ",muon1TrackIsol/muon1->inDetTrackParticle()->pt(), m_event_weight, m_event_flav); } } } } if (m_passCutsElectron && elec0) { float elec0CaloIsol = m_selectElec->caloIsol(elec0, egammaParameters::etcone20, false); m_osHistSvc->fillTH1Flav(m_histDirIsol+"ElecCaloIsol",elec0CaloIsol/1000., m_event_weight, m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirIsol+"ElecCaloIsolRel",elec0CaloIsol/elec0->pt(), m_event_weight, m_event_flav); if (m_passCutZMass) { m_osHistSvc->fillTH1Flav(m_histDirIsol+"ElecCaloIsolMZ",elec0CaloIsol/1000., m_event_weight, m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirIsol+"ElecCaloIsolRelMZ",elec0CaloIsol/elec0->pt(), m_event_weight, m_event_flav); } if (elec0->trackParticle()) { float elec0TrackIsol = m_selectElec->trackIsol(elec0, 0.2, false); m_osHistSvc->fillTH1Flav(m_histDirIsol+"ElecTrackIsol",elec0TrackIsol/1000., m_event_weight, m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirIsol+"ElecTrackIsolRel",elec0TrackIsol/elec0->trackParticle()->pt(), m_event_weight, m_event_flav); if (m_passCutZMass) { m_osHistSvc->fillTH1Flav(m_histDirIsol+"ElecTrackIsolMZ",elec0TrackIsol/1000., m_event_weight, m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirIsol+"ElecTrackIsolRelMZ",elec0TrackIsol/elec0->trackParticle()->pt(), m_event_weight, m_event_flav); } } if (elec1) { float elec1CaloIsol = m_selectElec->caloIsol(elec1, egammaParameters::etcone20, false); m_osHistSvc->fillTH1Flav(m_histDirIsol+"ElecCaloIsol",elec1CaloIsol/1000., m_event_weight, m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirIsol+"ElecCaloIsolRel",elec1CaloIsol/elec1->pt(), m_event_weight, m_event_flav); if (m_passCutZMass) { m_osHistSvc->fillTH1Flav(m_histDirIsol+"ElecCaloIsolMZ",elec1CaloIsol/1000., m_event_weight, m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirIsol+"ElecCaloIsolRelMZ",elec1CaloIsol/elec1->pt(), m_event_weight, m_event_flav); } if (elec1->trackParticle()) { float elec1TrackIsol = m_selectElec->trackIsol(elec1, 0.2, false); m_osHistSvc->fillTH1Flav(m_histDirIsol+"ElecTrackIsol",elec1TrackIsol/1000., m_event_weight, m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirIsol+"ElecTrackIsolRel",elec1TrackIsol/elec1->trackParticle()->pt(), m_event_weight, m_event_flav); if (m_passCutZMass) { m_osHistSvc->fillTH1Flav(m_histDirIsol+"ElecTrackIsolMZ",elec1TrackIsol/1000., m_event_weight, m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirIsol+"ElecTrackIsolRelMZ",elec1TrackIsol/elec1->trackParticle()->pt(), m_event_weight, m_event_flav); } } } } float dRMuon0(999.), dRMuon1(999.), dRElec0(999.), dRElec1(999.); bool dRMuon0isB(false), dRMuon1isB(false), dRElec0isB(false), dRElec1isB(false); int dRMuon0Flav(0), dRMuon1Flav(0), dRElec0Flav(0), dRElec1Flav(0); JetCollection::const_iterator jetItr (m_jets->begin()), jetItrE (m_jets->end()); for (; jetItr != jetItrE; ++jetItr) { if (!(*jetItr)) continue; int flav=GetFlavourHad(*jetItr); float tagWeight = (!m_btaggerName.empty()) ? (*jetItr)-> getFlavourTagWeight(m_btaggerName) : (*jetItr)-> getFlavourTagWeight(); if (muon0) { float dR = m_analysisTools->deltaR(muon0,*jetItr); if (dR < dRMuon0) { dRMuon0 = dR; dRMuon0isB = (tagWeight > m_jetWeight0); dRMuon0Flav = flav; } } if (muon1) { float dR = m_analysisTools->deltaR(muon1,*jetItr); if (dR < dRMuon1) { dRMuon1 = dR; dRMuon1isB = (tagWeight > m_jetWeight0); dRMuon1Flav = flav; } } if (elec0) { float dR = m_analysisTools->deltaR(elec0,*jetItr); if (dR < dRElec0) { dRElec0 = dR; dRElec0isB = (tagWeight > m_jetWeight0); dRElec0Flav = flav; } } if (elec1) { float dR = m_analysisTools->deltaR(elec1,*jetItr); if (dR < dRElec1) { dRElec1 = dR; dRElec1isB = (tagWeight > m_jetWeight0); dRElec1Flav = flav; } } } if (muon0) m_osHistSvc->fillTH1Flav(m_histDirIsol+"MuonJetDR",dRMuon0, m_event_weight, dRMuon0Flav); if (muon0 && dRMuon0isB) m_osHistSvc->fillTH1Flav(m_histDirIsol+"MuonJetDRWC",dRMuon0, m_event_weight, dRMuon0Flav); if (muon1) m_osHistSvc->fillTH1Flav(m_histDirIsol+"MuonJetDR",dRMuon1, m_event_weight, dRMuon1Flav); if (muon1 && dRMuon1isB) m_osHistSvc->fillTH1Flav(m_histDirIsol+"MuonJetDRWC",dRMuon1, m_event_weight, dRMuon1Flav); if (elec0) m_osHistSvc->fillTH1Flav(m_histDirIsol+"ElecJetDR",dRElec0, m_event_weight, dRElec0Flav); if (elec0 && dRElec0isB) m_osHistSvc->fillTH1Flav(m_histDirIsol+"ElecJetDRWC",dRElec0, m_event_weight, dRElec0Flav); if (elec1) m_osHistSvc->fillTH1Flav(m_histDirIsol+"ElecJetDR",dRElec1, m_event_weight, dRElec1Flav); if (elec1 && dRElec1isB) m_osHistSvc->fillTH1Flav(m_histDirIsol+"ElecJetDRWC",dRElec1, m_event_weight, dRElec1Flav); return; } int OSHbb::GetFlavourHad(Jet* jet) { //Hadron definition of jet flavour if( m_selectTruth->GetDjetNearestBHadron(jet->eta(),jet->phi())<0.4) { return 5; } if( m_selectTruth->GetDjetNearestCHadron(jet->eta(),jet->phi())<0.4) { return 4; } return 0; } int OSHbb::GetFlavourPart(Jet* jet) { //Parton definition of jet flavour if( m_selectTruth->GetDjetNearestBQuark(jet->eta(),jet->phi())<0.4) { return 5; } if( m_selectTruth->GetDjetNearestCQuark(jet->eta(),jet->phi())<0.4) { return 4; } return 0; } int OSHbb::GetFlavour(Jet* jet) { msg(MSG::DEBUG) << "GetFlavour" << endreq; std::string label("N/A"); const Analysis::TruthInfo* mcinfo = jet->tagInfo("TruthInfo"); if (mcinfo) label = mcinfo->jetTruthLabel(); int flav(0); if(label == "B") flav = 5; // b-jet if(label == "C") flav = 4; // c-jet if(label == "T") flav = 15; // tau return flav; } void OSHbb::TrackQuantities(){ //Fill some track quantities with main aim to reduce top background //since dileptonic decay should have a high pt track somewhere even //if it wasn't reconstructed as a lepton m_trPt0=0; m_trPt1=0; m_trPt2=0; // read the AOD track container from persistecy storage const Rec::TrackParticleContainer* trackTES(0); if (evtStore()->retrieve(trackTES, m_trackParticleContainerName).isFailure()) { msg(MSG::ERROR) << "Unable to retrieve " << m_trackParticleContainerName << endreq; return; } Rec::TrackParticleContainer::const_iterator trackItr=trackTES->begin(); Rec::TrackParticleContainer::const_iterator trackItrE=trackTES->end(); for(;trackItr !=trackItrE; ++trackItr){ const Rec::TrackParticle * trackParticle=(*trackItr); float trpt=trackParticle->pt(); float DeltaRMin=9.99; JetCollection::const_iterator jetItr (m_jets->begin()), jetItrE (m_jets->end()); for (; jetItr != jetItrE; ++jetItr) { Jet* jet = *jetItr; float DeltaR= jet->hlv().deltaR(trackParticle->hlv()); if(DeltaRsize()>0){ float DeltaR= m_electrons->at(0)->hlv().deltaR(trackParticle->hlv()); if(DeltaR<0.1) break; } if(m_muons->size()>0){ float DeltaR= m_muons->at(0)->hlv().deltaR(trackParticle->hlv()); if(DeltaR<0.1) break; } if(m_electrons->size()>1){ float DeltaR= m_electrons->at(1)->hlv().deltaR(trackParticle->hlv()); if(DeltaR<0.1) break; } if(m_muons->size()>1){ float DeltaR= m_muons->at(1)->hlv().deltaR(trackParticle->hlv()); if(DeltaR<0.1) break; } if(trpt> m_trPt0) { m_trPt2=m_trPt1; m_trPt1=m_trPt0; m_trPt0=trpt; } else if(trpt> m_trPt1) { m_trPt2=m_trPt1; m_trPt1=trpt; } else if(trpt> m_trPt2) { m_trPt2=trpt; } } //Stop tracks going to very high values. if(m_trPt0>100000.0) m_trPt0=100000.0; if(m_trPt1>100000.0) m_trPt1=100000.0; if(m_trPt2>100000.0) m_trPt2=100000.0; } void OSHbb::FillHt() { // Fill the Ht, transverse sum data members m_Ht=0; //transverse sum of all objects m_SEtjets=0; //transverse sum of all jets m_SEtnobjets=0; //transverse sum of all jets wo b jets //Add jets using low pt container // JetCollection* jetslowpt=m_selectJetlowpt->selectedJets(); // JetCollection::const_iterator jetItr (jetslowpt->begin()), jetItrE (jetslowpt->end()); //option to use standard jets JetCollection::const_iterator jetItr (m_jets->begin()), jetItrE (m_jets->end()); for (; jetItr != jetItrE; ++jetItr) { Jet* jet = *jetItr; m_SEtjets += m_selectJet->hlvCal(jet).perp(); } m_SEtnobjets=m_SEtjets-m_selectJet->hlvCal(m_jet0).perp()-m_selectJet->hlvCal(m_jet1).perp(); //Add et miss to Ht m_Ht =m_SEtjets+m_EtMiss; //Add leptons if(m_electrons->size()>0) m_Ht +=m_selectElec->hlvCal(m_electrons->at(0)).perp(); if(m_muons->size()>0) m_Ht += m_selectMuon->hlvCal(m_muons->at(0)).perp(); if(m_electrons->size()>1) m_Ht +=m_selectElec->hlvCal(m_electrons->at(1)).perp(); if(m_muons->size()>1) m_Ht += m_selectMuon->hlvCal(m_muons->at(1)).perp(); } void OSHbb::FillNNVar(float var){ //, std::string name) { fNNVariables[m_idxNN] = var; //fNNNames[m_idxNN] = name; ++m_idxNN; } void OSHbb::SetNNVariables(){ //Set the variables for use in an NN discriminant if(m_jet0<=0||m_jet1<=0) return; float jetEt0= (m_jet0 !=0) ? m_selectJet->hlvCal(m_jet0).et() :0; float jetEt1= (m_jet1 !=0) ? m_selectJet->hlvCal(m_jet1).et() :0; float jetEt2= (m_jet2 !=0) ? m_selectJet->hlvCal(m_jet2).et() :0; float jjDeltaR=m_selectJet->hlvCal(m_jet0).deltaR(m_selectJet->hlvCal(m_jet1)); float jjDeltaEta=acos(cos(m_selectJet->hlvCal(m_jet0).eta() - m_selectJet->hlvCal(m_jet1).eta())); int numjets=m_jets->size(); HepLorentzVector jj=m_selectJet->hlvCal(m_jet0)+m_selectJet->hlvCal(m_jet1); float phijj=jj.phi(); float mjj=jj.m(); float jjZDeltaPhi=acos(cos(m_phiZ-phijj)); // Cos Theta star TLorentzVector bj0; TLorentzVectorConv(m_selectJet->hlvCal(m_jet0),bj0); TLorentzVector bj1; TLorentzVectorConv(m_selectJet->hlvCal(m_jet1),bj1); TLorentzVector higgs=bj0+bj1; float costhetastar=coshel(bj0, higgs); //Make a vector of neutrino with Pz=0 TLorentzVector nu0=TLorentzVector(m_EtMiss_x,m_EtMiss_y,0,sqrt(m_EtMiss_x*m_EtMiss_x+m_EtMiss_y*m_EtMiss_y)); TLorentzVector nu1=TLorentzVector(m_EtMiss_x,m_EtMiss_y,0,sqrt(m_EtMiss_x*m_EtMiss_x+m_EtMiss_y*m_EtMiss_y)); // Lepton four vectors TLorentzVector lep0, lep1; if(m_electrons->size()>0) TLorentzVectorConv(m_selectElec->hlvCal(m_electrons->at(0)),lep0); else if(m_muons->size()>0) TLorentzVectorConv(m_selectMuon->hlvCal(m_muons->at(0)),lep0); else lep0.SetPtEtaPhiM(0.,0.,0.,0.); if(m_electrons->size()>1) TLorentzVectorConv(m_selectElec->hlvCal(m_electrons->at(1)),lep1); else if(m_muons->size()>1) TLorentzVectorConv(m_selectMuon->hlvCal(m_muons->at(1)),lep1); else lep1.SetPtEtaPhiM(0.,0.,0.,0.); float llDeltaEta=acos(cos(lep0.Eta() - lep1.Eta())); // Top mass bool solnu = NeutrinoWConstraint(lep0, nu0, nu1); float topmass0= (bj0+lep0+nu0).M(); float topmass1= (solnu)? (bj0+lep0+nu1).M() : 0; float topmass2= (bj1+lep0+nu0).M(); float topmass3= (solnu)? (bj1+lep0+nu1).M() : 0; std::cout <<"topmass0 "<< topmass0 << " topmass1 " << topmass1 <<"topmass2 "<< topmass2 << " topmass3 " << topmass3 <hlvCal(m_electrons->at(0)).phi()-m_selectElec->hlvCal(m_electrons->at(1)).phi())); m_idxNN=0; FillNNVar(m_selectJet->hlvCal(m_jet0).perp()); FillNNVar(m_selectJet->hlvCal(m_jet0).eta()); FillNNVar(m_selectJet->hlvCal(m_jet0).phi()); FillNNVar(m_selectJet->hlvCal(m_jet0).m()); FillNNVar((!m_btaggerName.empty()) ? m_jet0->getFlavourTagWeight(m_btaggerName) : m_jet0-> getFlavourTagWeight()); //FillNNVar(m_jetWeight0); FillNNVar(m_selectJet->hlvCal(m_jet1).perp()); FillNNVar(m_selectJet->hlvCal(m_jet1).eta()); FillNNVar(m_selectJet->hlvCal(m_jet1).phi()); FillNNVar(m_selectJet->hlvCal(m_jet1).m()); FillNNVar((!m_btaggerName.empty()) ? m_jet1->getFlavourTagWeight(m_btaggerName) : m_jet1-> getFlavourTagWeight()); //FillNNVar(m_jetWeight1); if (m_jet2) { FillNNVar(m_selectJet->hlvCal(m_jet2).perp()); FillNNVar(m_selectJet->hlvCal(m_jet2).eta()); FillNNVar(m_selectJet->hlvCal(m_jet2).phi()); FillNNVar(m_selectJet->hlvCal(m_jet2).m()); FillNNVar((!m_btaggerName.empty()) ? m_jet2->getFlavourTagWeight(m_btaggerName) : m_jet2-> getFlavourTagWeight()); //FillNNVar(m_jetWeight2); } else { FillNNVar(0.); FillNNVar(0.); FillNNVar(0.); FillNNVar(0.); FillNNVar(0.); } FillNNVar(lep0.Pt()); FillNNVar(lep0.Eta()); FillNNVar(lep0.Phi()); FillNNVar(lep0.M()); FillNNVar(lep1.Pt()); FillNNVar(lep1.Eta()); FillNNVar(lep1.Phi()); FillNNVar(lep1.M()); FillNNVar(m_ptZ); FillNNVar(m_mZ); FillNNVar(jjDeltaEta); FillNNVar(m_JJDeltaPhiOPt); FillNNVar(llDeltaEta); FillNNVar(m_LLDeltaPhi); FillNNVar(jjZDeltaPhi); FillNNVar(m_Ht); FillNNVar(m_EtMiss); FillNNVar(numjets); FillNNVar(mjj); FillNNVar(fabs(costhetastar)); FillNNVar(topmass); FillNNVar(m_SEtjets); FillNNVar(m_SEtnobjets); FillNNVar(m_event_weight); m_osHistSvc->fillTree("/file1/NNVariables"); if(mjj<40000.0) mjj=40000.0; float jetEt0Sc=jetEt0/mjj; float jetEt1Sc=jetEt1/mjj; float HtSc=m_Ht/mjj; //Fill historams of inputs m_histDirNN = fTHistStreamName+"NN/"; m_osHistSvc->fillTH1Flav(m_histDirNN+"NNptZ",m_ptZ/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNjetEt0Sc",jetEt0Sc, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNjetEt1Sc",jetEt1Sc, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNjetEt0",jetEt0/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNjetEt1",jetEt1/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNjjDeltaR",jjDeltaR, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNjjZDeltaPhi",jjZDeltaPhi, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNHtSc",HtSc, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNHt",m_Ht/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNmjj",mjj/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNEtmiss",m_EtMiss/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNnumjets",numjets, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNmZ",m_mZ/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNtrPt0",m_trPt0/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNtrPt1",m_trPt1/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNjetEt2",jetEt2/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNSEtjets",m_SEtjets/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNSEtnobjets",m_SEtnobjets/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNtopmassAll",topmass0/GeV, m_event_weight,m_event_flav); // if(solnu) m_osHistSvc->fillTH1Flav(m_histDirNN+"NNtopmassAll",topmass1/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNtopmassAll",topmass2/GeV, m_event_weight,m_event_flav); // if(solnu) m_osHistSvc->fillTH1Flav(m_histDirNN+"NNtopmassAll",topmass3/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNtopmassNear",topmass/GeV, m_event_weight,m_event_flav); if(NBJets() ==2 && m_passCutThirdJetWeight) { m_osHistSvc->fillTH1Flav(m_histDirNN+"NNptZWC",m_ptZ/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNjetEt0ScWC",jetEt0Sc, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNjetEt1ScWC",jetEt1Sc, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNjetEt0WC",jetEt0/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNjetEt1WC",jetEt1/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNjjDeltaRWC",jjDeltaR, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNjjZDeltaPhiWC",jjZDeltaPhi, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNHtScWC",HtSc, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNHtWC",m_Ht, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNmjjWC",mjj/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNEtmissWC",m_EtMiss/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNnumjetsWC",numjets, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNmZWC",m_mZ/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNtrPt0WC",m_trPt0/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNtrPt1WC",m_trPt1/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNjetEt2WC",jetEt2/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNSEtjetsWC",m_SEtjets/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNSEtnobjetsWC",m_SEtnobjets/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNtopmassAllWC",topmass0/GeV, m_event_weight,m_event_flav); if(solnu) m_osHistSvc->fillTH1Flav(m_histDirNN+"NNtopmassAllWC",topmass1/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNtopmassAllWC",topmass2/GeV, m_event_weight,m_event_flav); if(solnu) m_osHistSvc->fillTH1Flav(m_histDirNN+"NNtopmassAllWC",topmass3/GeV, m_event_weight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirNN+"NNtopmassNearWC",topmass/GeV, m_event_weight,m_event_flav); } //Cut based analysis if(m_VectorBosonOp==eW){ m_cutsloose = mjj>40000.0&&jjDeltaR<3&&jjZDeltaPhi>2.75&&jetEt2<20000&&m_SEtnobjets<40000; m_cutsmedium = mjj>40000.0&&jjDeltaR<3&&jjZDeltaPhi>2.75&&jetEt2<20000&&m_SEtnobjets<40000&&numjets<=3; m_cutstight = mjj>40000.0&&jjDeltaR<3&&jjZDeltaPhi>2.75&&m_SEtnobjets<40000&&numjets==2; } else if (m_VectorBosonOp==eZnu){ m_cutsloose = mjj>40000.0&&m_ptZ>100000.0&&jetEt0>0.6&&jetEt1>0.2&&jjDeltaR<1.5&&jetEt2<20000&&m_SEtnobjets<40000; m_cutsmedium = mjj>40000.0&&m_ptZ>150000.0&&jetEt0>0.6&&jetEt1>0.2&&jjDeltaR<1.5&&jetEt2<20000&&m_SEtnobjets<40000; m_cutstight = mjj>40000.0&&m_ptZ>200000.0&&jetEt0>0.6&&jetEt1>0.2&&jjDeltaR<1.5&&jetEt2<20000&&m_SEtnobjets<40000; } else if (m_VectorBosonOp==eZ){ m_cutsloose = mjj>40000.0&&m_ptZ>100000.0&&jetEt0>0.6&&jetEt1>0.2&&jjDeltaR<1.5&&jjZDeltaPhi>2.75&&jetEt2<20000&&m_SEtnobjets<40000; m_cutsmedium = mjj>40000.0&&m_ptZ>150000.0&&jetEt0>0.6&&jetEt1>0.2&&jjDeltaR<1.5&&jjZDeltaPhi>2.75&&jetEt2<20000&&m_SEtnobjets<40000; m_cutstight = mjj>40000.0&&m_ptZ>200000.0&&jetEt0>0.6&&jetEt1>0.2&&jjDeltaR<1.5&&jjZDeltaPhi>2.75&&jetEt2<20000&&m_SEtnobjets<40000; } } bool OSHbb::ApplyCuts(bool fillhistos) { //Make sure histograms have been booked in case ApplyCuts() is called from child of OSHbb. //(because the cut methods fill histograms) if (!m_bookedHistos) OSHbb::BookHistos(); return ( ApplyMuonCuts(fillhistos) || ApplyElectronCuts(fillhistos) ) && ZLeptonMassCut() && METCuts(fillhistos) && ApplyJetCuts(fillhistos); } bool OSHbb::ApplyMuonCuts(bool fillhistos) { return (m_VectorBosonOp==eW)? WCutsMuon(fillhistos) : (m_VectorBosonOp==eZ) ? ZCutsMuon(fillhistos) : (m_VectorBosonOp==eZnu) ? ZNuCutsMuon(fillhistos) : false; } bool OSHbb::ApplyElectronCuts(bool fillhistos) { return (m_VectorBosonOp==eW)? WCutsElectron(fillhistos) : (m_VectorBosonOp==eZ)? ZCutsElectron(fillhistos) : (m_VectorBosonOp==eZnu)? ZNuCutsElectron(fillhistos) : false; } bool OSHbb::ZLeptonMassCut() { //This now does not cut the events in the Z mass sidebands to allow us to plot backgrounds for the side bands msg(MSG::DEBUG) << "ZLeptonMassCut()" << endreq; m_passCutZMass = true; if(m_VectorBosonOp==eZ&&(m_mZm_zmasshigh)) m_passCutZMass=false; else if(m_VectorBosonOp==eW&&(m_mZm_zmasslowsb&&m_mZm_zbmasslow&&m_dibmassm_zbmasslow&&m_dibmassopt1 && ZBJetMassCut(); } bool OSHbb::JetCountCut(bool fillhistos) { //Cut on number of jets //Return 0 this point if number of jets >3 (or >2 for Znunu analysis) int numjets = m_selectJet->numSelected(); // if(m_VectorBosonOp!=eZ&&numjets>2) return 0; // disabled for W case for now //if(m_VectorBosonOp==eZnu&&numjets>2) return 0; if(fillhistos&&m_passCutZMass&&m_passCutEtmiss) { m_event_counter[eECNumJetsLT4]+= m_event_weight; FillllqqCutHisto(ellqqCutJetVeto, m_event_weight, m_event_flav); FilllnubbCutHisto(elnubbCutJetVeto, m_event_weight, m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirCuts+"Cut", eCutExtraJets, m_event_weight,m_event_flav); if(m_zbzlgen) m_osHistSvc->fillTH1Flav(m_histDirCuts+"CutZb", eCutExtraJets, m_event_weight,m_event_flav); } return 1; } bool OSHbb::JetVertexFractionCut(bool fillhistos) { //Cut on Jet Vertex Fraction to remove pileup jets float jvf_cut=0.75; JetCollection::const_iterator jetItr (m_jets->begin()), jetItrE (m_jets->end()); for (; jetItr != jetItrE; ++jetItr) { Jet* jet = *jetItr; //if(fabs(jet->eta())>2.5) continue; float jvf=jet->getMoment("JVF"); if(fillhistos&&m_passCutZMass&&m_passCutEtmiss) { m_osHistSvc->fillTH1Flav(m_histDirJets+"JetVertexFraction", jvf, m_event_weight,m_event_flav); if(fabs(jet->eta())<2.5) m_osHistSvc->fillTH1Flav(m_histDirJets+"JetVertexFractionBarrel", jvf, m_event_weight,m_event_flav); if(fabs(jet->eta())<2.1) m_osHistSvc->fillTH1Flav(m_histDirJets+"JetVertexFractionCentral", jvf, m_event_weight,m_event_flav); else m_osHistSvc->fillTH1Flav(m_histDirJets+"JetVertexFractionForward", jvf, m_event_weight,m_event_flav); int numjets= m_selectJet->numSelected(); if (numjets>1) { m_osHistSvc->fillTH1Flav(m_histDirJets+"JetVertexFractionDijet", jvf, m_event_weight,m_event_flav); if(fabs(jet->eta())<2.5) m_osHistSvc->fillTH1Flav(m_histDirJets+"JetVertexFractionDijetBarrel", jvf, m_event_weight,m_event_flav); if(fabs(jet->eta())<2.1) m_osHistSvc->fillTH1Flav(m_histDirJets+"JetVertexFractionDijetCentral", jvf, m_event_weight,m_event_flav); else m_osHistSvc->fillTH1Flav(m_histDirJets+"JetVertexFractionDijetForward", jvf, m_event_weight,m_event_flav); } } // Keep jets with no track associated, where JVF=-1 if(fabs(jvf)Remove(jet); } //Same for lowpt jets JetCollection* jetslowpt=m_selectJetlowpt->selectedJets(); JetCollection::const_iterator jetlItr (jetslowpt->begin()), jetlItrE (jetslowpt->end()); for (; jetlItr != jetlItrE; ++jetlItr) { Jet* jet = *jetlItr; //if(fabs(jet->eta())>2.5) continue; float jvf=fabs(jet->getMoment("JVF")); // Keep jets with no track associated, where JVF=-1 if(jvfRemove(jet); } //Do not reject events since jets are already rejected return 1; } bool OSHbb::JetCleaningCut(bool fillhistos) { // fill cut histo for all jet events, regardless of MZ and Etmiss cut if(fillhistos&&m_passCutEtmiss&&m_passCutZMass) m_osHistSvc->fillTH1Flav(m_histDirCuts+"Cut", eCutJetAll, m_event_weight,m_event_flav); int nbadjets = m_selectJet->nBad(); //int nuglyjets = m_selectJet->nUgly(); if(nbadjets != 0 ) return 0; // if(nbadjets != 0 || nuglyjets != 0 ) return 0; //if(m_VectorBosonOp==eZnu&&numjets>2) return 0; if(fillhistos&&m_passCutEtmiss&&m_passCutZMass) { m_osHistSvc->fillTH1Flav(m_histDirCuts+"Cut", eCutJetClean, m_event_weight,m_event_flav); if(m_zbzlgen) m_osHistSvc->fillTH1Flav(m_histDirCuts+"CutZb", eCutJetClean, m_event_weight,m_event_flav); } return 1; } bool OSHbb::FirstJetWeightCut(bool fillhistos){ //Cut on 1st Jet Weight as a top veto //Return 0 if 1st jet might be a b or c if(m_jetWeight0>0) return 0; if(fillhistos&&m_passCutEtmiss&&m_passCutZMass) m_osHistSvc->fillTH1Flav(m_histDirCuts+"Cut", eCut1stJetWC, m_event_weight,m_event_flav); return 1; } bool OSHbb::ThirdJetWeightCut(bool fillhistos ){ //Cut on 3rd Jet Weight //Return 0 if 3rd jet might be a b or c if(m_jetWeight2>m_jetWeight2Cut) return 0; if(fillhistos&&m_passCutEtmiss&&m_passCutZMass) { m_event_counter[eECJetWeight2]+= m_event_weight; FillllqqCutHisto(ellqqCutBJetVeto, m_event_weight, m_event_flav); FilllnubbCutHisto(elnubbCutBJetVeto, m_event_weight, m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirCuts+"Cut", eCut3rdJetWC, m_event_weight,m_event_flav); if(m_zbzlgen) m_osHistSvc->fillTH1Flav(m_histDirCuts+"CutZb", eCut3rdJetWC, m_event_weight,m_event_flav); } return 1; } bool OSHbb::BJetsBackToBackCut(bool fillhistos ){ m_JJDeltaPhi= (m_bjet0>0&&m_bjet1>0) ? (m_bjet0->hlv()).v().deltaPhi((m_bjet1->hlv()).v()) : 0; m_JJDeltaPhiOPt= (m_jet0>0&&m_jet1>0) ? (m_jet0->hlv()).v().deltaPhi((m_jet1->hlv()).v()) : 0; //Sign doesn't matter m_JJDeltaPhi= fabs(m_JJDeltaPhi); m_JJDeltaPhiOPt= fabs(m_JJDeltaPhiOPt); //Fill counter with delta phi cuts. if(fillhistos&& m_JJDeltaPhi2.0944) return 0; if(m_JJDeltaPhi>1.5707) return 0; //or if there is a lot of energy outside the 2 leading jets //if(m_SEtnobjets>50*GeV) return 0; } if(fillhistos&&m_passCutEtmiss&&m_passCutZMass) { m_osHistSvc->fillTH1Flav(m_histDirCuts+"Cut", eCutDPhiJJ, m_event_weight,m_event_flav); if(m_zbzlgen) m_osHistSvc->fillTH1Flav(m_histDirCuts+"CutZb", eCutDPhiJJ, m_event_weight,m_event_flav); } return 1; } int OSHbb::NBJets(){ //Cut on b-jet weights. int nbjets=0; if(m_jetWeight0>m_jetWeight0Cut) nbjets++; if(m_jetWeight1>m_jetWeight1Cut) nbjets++; if(m_jetWeight2>m_jetWeight2Cut) nbjets++; return nbjets; } StatusCode OSHbb::FetchToolsAndServices(){ msg(MSG::INFO) << "OSHbb::FetchToolsAndServices()" << endreq; // event timing if (service("ChronoStatSvc", m_chronoSvc).isFailure()) { msg(MSG::ERROR) << "Unable to retrieve pointer to ChronoStatSvc" << endreq; return StatusCode::FAILURE; } if (m_selectElec.retrieve().isFailure()) { msg(MSG::ERROR) << "Can't get handle on electron selector" << endreq; return StatusCode::FAILURE; } if (m_selectPhot.retrieve().isFailure()) { msg(MSG::ERROR) << "Can't get handle on photon selector" << endreq; return StatusCode::FAILURE; } if (m_selectMuon.retrieve().isFailure()) { msg(MSG::ERROR) << "Can't get handle on muon selector" << endreq; return StatusCode::FAILURE; } if (m_selectMuon.retrieve().isFailure()) { msg(MSG::ERROR) << "Can't get handle on muon selector" << endreq; return StatusCode::FAILURE; } if (m_selectTau.retrieve().isFailure()) { msg(MSG::ERROR) << "Can't get handle on tau selector" << endreq; return StatusCode::FAILURE; } if (m_selectTruth.retrieve().isFailure()) { msg(MSG::ERROR) << "Can't get handle on truth selector" << endreq; return StatusCode::FAILURE; } if (m_selectJet.retrieve().isFailure()) { msg(MSG::ERROR) << "Can't get handle on jet selector" << endreq; return StatusCode::FAILURE; } if (m_selectJetlowpt.retrieve().isFailure()) { msg(MSG::ERROR) << "Can't get handle on jet low pt selector" << endreq; return StatusCode::FAILURE; } if (m_selectJetformet.retrieve().isFailure()) { msg(MSG::ERROR) << "Can't get handle on jet for met selector" << endreq; return StatusCode::FAILURE; } if (m_osWeightSvc.retrieve().isFailure()) { msg(MSG::ERROR) << "Unable to retrieve pointer to WeightSvc" << endreq; return StatusCode::FAILURE; } if (service( "THistSvc", fTHistSvc, true).isFailure()) { msg(MSG::ERROR) << "THistSvc could NOT be loaded." << endreq; return StatusCode::FAILURE; } if (m_trigDec.retrieve().isFailure() ){ msg(MSG::ERROR) << "Can't get handle on TrigDecisionTool" << endreq; return StatusCode::FAILURE; } if (m_beamCondSvc.retrieve().isFailure()) { msg(MSG::WARNING) << "Failed to retrieve beamspot service" << endreq; return StatusCode::FAILURE; } if (m_kinFitter.retrieve().isFailure()) { msg(MSG::WARNING) << "Failed to retrieve kinematic fitter tool" << endreq; return StatusCode::FAILURE; } return StatusCode::SUCCESS; } void OSHbb::FindReconstructedParticles() { //Find Electrons msg(MSG::DEBUG) << "Elecs()" << endreq; m_chronoSvc->chronoStart("SelectElectron"); Check( m_selectElec->select() ); m_electrons = m_selectElec->selectedElecs(); m_chronoSvc->chronoStop("SelectElectron"); //Find Jets msg(MSG::DEBUG) << "Jets()" << endreq; m_chronoSvc->chronoStart("SelectJet"); m_selectJet->setElecsForOverlap(m_electrons); m_selectJetlowpt->setElecsForOverlap(m_electrons); m_selectJetformet->setElecsForOverlap(m_electrons); Check( m_selectJet->select() ); Check( m_selectJetlowpt->select() ); Check( m_selectJetformet->select() ); m_jets =m_selectJet->selectedJets(); m_jetsformet =m_selectJetformet->selectedJets(); // m_jetsUD=*m_selectJet->selectedJetsUD(); m_chronoSvc->chronoStop("SelectJet"); //comment next line if you want to try without this cut out non-isolated muons m_selectMuon->setJetsForOverlap(m_jets);//cut out non-isolated muons //Find Muons msg(MSG::DEBUG) << "Muons()" << endreq; m_chronoSvc->chronoStart("SelectMuon"); Check( m_selectMuon->select() ); m_muons = m_selectMuon->selectedMuons(); // m_muonsUD=*m_selectMuon->selectedMuonsUD(); m_chronoSvc->chronoStop("SelectMuon"); //Find Taus (tau cuts are removed) //m_selectTau->setElecsForOverlap(m_electrons); //m_selectTau->setMuonsForOverlap(m_muons); //m_selectTau->setJetsForOverlap(m_jets); //Check( m_selectTau->select() ); //m_taus = m_selectTau->selectedMuons(); //Set Candidate Lepton+Jet pointers. msg(MSG::DEBUG) << "leptons()" << endreq; m_chronoSvc->chronoStart("FindLeptons"); FindLeptons(); m_chronoSvc->chronoStop("FindLeptons"); //Calculate and Set MET Variables m_chronoSvc->chronoStart("MET"); m_metsyshadshiftdo=false; m_metsyshadshiftup=false; m_metsyshadsmear=false; m_metsysfwhadshiftdo=false; m_metsysfwhadshiftup=false; m_metsysfwhadsmear=false; msg(MSG::DEBUG) << "MET()" << endreq; CalculateMissingEt(); m_chronoSvc->chronoStop("MET"); return; } void OSHbb::FindLeptons() { //Find leptons. //Lepton pointers m_e0 = 0; m_e1 = 0; m_mu0 = 0; m_mu1 = 0; m_l0 = 0; m_l1 = 0; {//Always select the highest PT pair of electrons. // sort(m_electrons->begin(),m_electrons->end(),&OSHbb::SortPt); if(m_electrons->size()>0) m_e0 = m_electrons->at(0); if(m_electrons->size()>1) m_e1 = m_electrons->at(1); } {//Always select the highest PT pair of muons. // sort(m_muons->begin(),m_muons->end(),&OSHbb::SortPt); if(m_muons->size()>0) m_mu0 = m_muons->at(0); if(m_muons->size()>1) m_mu1 = m_muons->at(1); } //Set the selected Higgs Lepton pointers. if(m_mu0 && m_e0) {//If both electrons and muons are found, choose the higher pt lepton flavour. if( m_e0->pt() > m_mu0->pt() ) {//Highest Pt lepton is an electron m_l0 = m_e0; m_l1 = m_e1; } else {//Highest Pt lepton is a muon m_l0 = m_mu0; m_l1 = m_mu1; } } else if (m_e0) {//Only electrons found m_l0 = m_e0; m_l1 = m_e1; } else {//Only muons found m_l0 = m_mu0; m_l1 = m_mu1; } return; } void OSHbb::FindJets() { msg(MSG::DEBUG) << "FindJets()" << endreq; //Set Jet pointers to null. m_bjet0 = 0; m_bjet1 = 0; m_bjet2 = 0; m_jet0 = 0; m_jet1 = 0; m_jet2 = 0; m_jetWeight0 =-99; m_jetWeight1 =-99; m_jetWeight2 =-99; //Set Highest ET Jets //sort(m_jets->begin(),m_jets->end(),OSHbb::SortPt);//sort by PT if(m_jets->size()>0) m_jet0 = m_jets->at(0); if(m_jets->size()>1) m_jet1 = m_jets->at(1); if(m_jets->size()>2) m_jet2 = m_jets->at(2); //Set Highest Bweight Jets //sort(m_jets->begin(),m_jets->end(),OSHbb::SortBWeight);//Sort by B weight m_bjets=new JetCollection(*m_jets); sort(m_bjets->begin(),m_bjets->end(),OSHbb::SortBWeight);//Sort by B weight if(m_jets->size()>0) { m_bjet0 = m_bjets->at(0); if(fabs(m_selectJet->hlvCal(m_bjet0).eta()) getFlavourTagWeight(m_btaggerName) : m_bjet0-> getFlavourTagWeight(); } else m_jetWeight0 =-10; m_flav0 = GetFlavourHad(m_bjet0); } if(m_jets->size()>1) { m_bjet1 = m_bjets->at(1); if(fabs(m_selectJet->hlvCal(m_bjet1).eta()) getFlavourTagWeight(m_btaggerName) : m_bjet1-> getFlavourTagWeight(); } else m_jetWeight1 =-10; m_flav1 = GetFlavourHad(m_bjet1); } if(m_jets->size()>2) { m_bjet2 = m_bjets->at(2); if(fabs(m_selectJet->hlvCal(m_bjet2).eta()) getFlavourTagWeight(m_btaggerName) : m_bjet2-> getFlavourTagWeight(); } else m_jetWeight2 =-10; m_flav2 = GetFlavourHad(m_bjet2); } return; } /* void OSHbb::CalculateMissingEt() { // Get/Calculate MissingEt m_sumEtMiss = 0; m_EtMiss = 0; m_EtMiss_x = 0.; m_EtMiss_y = 0.; bool useRefFinal=1; // use ref final, else for early data Calo(LocHadTopo)+muons by hand const MissingET * pMissing = 0; std::string metcont("MET_RefFinal"); if (!useRefFinal) metcont="MET_LocHadTopo"; if (evtStore()->retrieve(pMissing, metcont).isFailure()) { msg(MSG::ERROR) << "Could not retrieve Missing ET Object" << endreq; return; } m_EtMissOrig = pMissing->et(); if (useRefFinal) { //This is supposedly the definition of EtMiss //MET_Ex(y)MissRefFinal = MET_Ex(y)MissRefEle + MET_Ex(y)MissRefGamma + MET_Ex(y)MissRefTau + MET_Ex(y)MissRefJet + MET_Ex(y)MissRefMuon + MET_Ex(y)MissCellOut + MET_Ex(y)MissMuonBoy + MET_Ex(y)MissCryo m_sumEtMiss = pMissing->sumet(); m_EtMiss = pMissing->et(); m_EtMiss_x = pMissing->etx(); m_EtMiss_y = pMissing->ety(); } else { // Early data retrieve info from Eta Regions const MissingEtRegions* caloReg = pMissing->getRegions(); if ( caloReg != 0 ) { // for the Central Region ( |eta| < 1.5) double ExMissTopoCENTReg = caloReg->exReg(MissingEtRegions::Central); double EyMissTopoCENTReg = caloReg->eyReg(MissingEtRegions::Central); double EtSumTopoCENTReg = caloReg->etSumReg(MissingEtRegions::Central); // for the EndCap region ( 1.5 < |eta| < 3.2) double ExMissTopoECReg = caloReg->exReg(MissingEtRegions::EndCap); double EyMissTopoECReg = caloReg->eyReg(MissingEtRegions::EndCap); double EtSumTopoECReg = caloReg->etSumReg(MissingEtRegions::EndCap); // for the Forward region (default value: |eta| < 5 => set to |eta|<4.5 in 7 TeV data) double ExMissTopoFWReg = caloReg->exReg(MissingEtRegions::Forward); double EyMissTopoFWReg = caloReg->eyReg(MissingEtRegions::Forward); double EtSumTopoFWReg = caloReg->etSumReg(MissingEtRegions::Forward); // m_sumEtMiss = EtSumTopoCENTReg + EtSumTopoECReg + EtSumTopoFWReg; m_EtMiss_x = ExMissTopoCENTReg + ExMissTopoECReg + ExMissTopoFWReg; m_EtMiss_y = EyMissTopoCENTReg + EyMissTopoECReg + EyMissTopoFWReg; m_EtMiss = sqrt(m_EtMiss_x*m_EtMiss_x+m_EtMiss_y*m_EtMiss_y); } } // Fix muon term ... if (useRefFinal) { const MissingET * pMissingMuon = 0; if (evtStore()->retrieve(pMissingMuon, "MET_MuonBoy").isFailure()) { msg(MSG::ERROR) << "Could not retrieve MuonBoy Missing ET Object" << endreq; return; } const MissingET * pMissingMuCalo = 0; if (evtStore()->retrieve(pMissingMuCalo, "MET_RefMuon").isFailure()) { msg(MSG::ERROR) << "Could not retrieve RefMuon Missing ET Object" << endreq; return; } const MuonContainer* muonContainer (0); //Now use our selected muons muonContainer=m_muons; // if (evtStore()->retrieve(muonContainer, "StacoMuonCollection").isFailure() || !muonContainer) { // msg(MSG::WARNING) << "No muon container called found " << endreq; // return ; // } // Subtract off MET_MuonBoy and non-isolated calo term m_EtMiss_x -= (pMissingMuon->etx() + pMissingMuCalo->etx()); m_EtMiss_y -= (pMissingMuon->ety() + pMissingMuCalo->ety()); // Add back in combined muons passing our Pt and Eta selection (use AOD muons so before isolation) MuonContainer::const_iterator muonItr(muonContainer->begin()), muonEnd(muonContainer->end()); for (; muonItr != muonEnd; ++muonItr) { Muon* muon(*muonItr); if (!muon) continue; // if (muon->pt() < 15*GeV || fabs(muon->eta()) > 2.5) continue; // if (muon->isStandAloneMuon()) continue; // Need to be careful with calo term depending on if isol/non-isol muons as they use comb/standalone muons // float muonCaloEt = muon->energyLoss().first*muon->sinTh(); float muonEt = m_selectMuon->hlvCal(muon).perp(); // - muonCaloEt; m_EtMiss_x -= muonEt * std::cos(m_selectMuon->hlvCal(muon).phi()); m_EtMiss_y -= muonEt * std::sin(m_selectMuon->hlvCal(muon).phi()); } // Scale jet term for calibration const MissingET * pMissingJet = 0; if (evtStore()->retrieve(pMissingJet, "MET_RefJet").isFailure()) { msg(MSG::ERROR) << "Could not retrieve MET_RefJet Missing ET Object" << endreq; return; } // Add in extra scaled energy (e.g. 1.05 - 1 = .05 = 5%) in same direction float scale = (m_selectJet->jetScale() - 1.); m_EtMiss_x += pMissingJet->etx() * scale; m_EtMiss_y += pMissingJet->ety() * scale; } else { // early data const MissingET * pMissingMuon = 0; if (evtStore()->retrieve(pMissingMuon, "MET_MuonBoy").isFailure()) { msg(MSG::ERROR) << "Could not retrieve MuonBoy Missing ET Object" << endreq; return; } const MissingET * pMissingMuTrack = 0; if (evtStore()->retrieve(pMissingMuTrack, "MET_RefMuon_Track").isFailure()) { msg(MSG::ERROR) << "Could not retrieve RefMuon_Track Missing ET Object" << endreq; return; } // Add MET_MuonBoy(track+spect) and remove calo contribution of isolated muons m_EtMiss_x += (pMissingMuon->etx() - pMissingMuTrack->etx()); m_EtMiss_y += (pMissingMuon->ety() - pMissingMuTrack->ety()); } m_EtMiss=sqrt(m_EtMiss_x*m_EtMiss_x+m_EtMiss_y*m_EtMiss_y); return; } */ /* void OSHbb::FillFakeHistos() { //Fill histograms to be used for fake rates // for fakes in data we only want JetTauEtmiss/calo stream, MC will have no stream //CBG if ( m_Stream!=eJetTauEtmiss || m_Stream!=eL1Calo || m_Stream!=eNoStream ) return; static const float jetEscale=0.85; //ratio of electron scale to jet scale //This is the standard scale for jets // static const float jetscale=m_selectJet->jetScale(); Check( m_selectJetformet->select() ); int numjets= m_selectJetforfakes->numSelected(); if(numjets<2) return; //Can't do anything with only 1 jet JetCollection* jets =m_selectJetforfakes->selectedJets(); // OSJetUserDataContainer jetsUD=*m_selectJetforfakes->selectedJetsUD(); float zfakeweight=0; // for fake histograms require the J5 trigger // approx 50% efficient at pt=25 GeV, 100% at ~45 GeV std::string J5("EF_j20_jetNoCut"); if (m_periodABCD) J5="L1_J5"; if (m_isMC) J5="EF_j20"; std::string MinBias("EF_L1MinBias_NoAlg"); if (m_periodABCD) MinBias="L1_MBTS_1"; if (m_isMC) MinBias="L1_MBTS_1"; bool passL1_J5=m_trigDec->isPassed(J5); // actual bool passL1_J5_raw=m_trigDec->isPassedBits(J5)&TrigDefs::L1_isPassedBeforePrescale; // raw bool passL1_MBTS=m_trigDec->isPassed(MinBias); // monitor //This is to find Z->ee fakes from jets JetCollection::const_iterator jetItr0 (jets->begin()), jetItrE (jets->end()); // loop over the AOD jet container for (; jetItr0 != jetItrE; ++jetItr0) { float zmass =0; Jet* jet0 = *jetItr0; double ptj = jet0->pt(); //double etaj = jet->eta(); //Plot Et spectrum before and after electron cut m_osHistSvc->fillTH1Flav(m_histDirFake+"PtJet", ptj/GeV, m_event_weight,m_event_flav); if (passL1_J5) m_osHistSvc->fillTH1Flav(m_histDirFake+"PtJetTrig", ptj/GeV, m_event_weight,m_event_flav); // monitor trigger efficiency if (passL1_MBTS) m_osHistSvc->fillTH1Flav(m_histDirFake+"PtJetTrigMon", ptj/GeV, m_event_weight,m_event_flav); // trigger efficiency if (passL1_MBTS && passL1_J5) m_osHistSvc->fillTH1Flav(m_histDirFake+"PtJetTrigActual", ptj/GeV, m_event_weight,m_event_flav); if (passL1_MBTS&&passL1_J5_raw) m_osHistSvc->fillTH1Flav(m_histDirFake+"PtJetTrigRaw", ptj/GeV, m_event_weight,m_event_flav); //Only fill the test if J5. This will limit stats badly // if(!passL1_J5) continue; int nume=m_selectElec->numSelected(); for(int i=0; iat(i)->eta()-jet0->eta(); float dphi=acos(cos(m_electrons->at(i)->phi()-jet0->phi())); float dr=sqrt(deta*deta+dphi*dphi); if(dr<0.4){ m_osHistSvc->fillTH1Flav(m_histDirFake+"PtJetPassE", ptj/GeV, m_event_weight,m_event_flav); float pteptjrat=m_electrons->at(i)->pt()/jet0->pt(); m_osHistSvc->fillTH1Flav(m_histDirFake+"PtEPtJetRat", pteptjrat, m_event_weight,m_event_flav); } } //Same for muons int nummu=m_selectMuon->numSelected(); for(int i=0; iat(i)->eta()-jet0->eta(); float dphi=acos(cos(m_muons->at(i)->phi()-jet0->phi())); float dr=sqrt(deta*deta+dphi*dphi); if(dr<0.4){ m_osHistSvc->fillTH1Flav(m_histDirFake+"PtJetPassMu", ptj/GeV, m_event_weight,m_event_flav); float pteptjrat=m_muons->at(i)->pt()/jet0->pt(); m_osHistSvc->fillTH1Flav(m_histDirFake+"PtMuPtJetRat", pteptjrat, m_event_weight,m_event_flav); } } if(jet0->pt()<20*GeV) continue; // std::cout << " ptj " << std::endl; //start counting at 1 more than first jet to loop over all pairs of jets JetCollection::const_iterator jetItr1 (jetItr0); // loop over the AOD jet container jetItr1++; for (; jetItr1 != jetItrE; ++jetItr1) { Jet* jet1 = *jetItr1; if(jet1->pt()<20*GeV) continue; zfakeweight=GetEleFakeRate(jet0->pt())*GetEleFakeRate(jet1->pt()); //correct for difference between jet scale and zmass = jetEscale*(jet0->hlv()+jet1->hlv()).m(); //std::cout << " zmass " << zmass << std::endl; if(m_passCutEtmiss) { m_osHistSvc->fillTH1Flav(m_histDirZ+"MZFake", zmass/GeV, m_event_weight*zfakeweight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirZ+"MZeeFake", zmass/GeV, m_event_weight*zfakeweight,m_event_flav); } // if(numjets<4) continue; //work out z->bb mass and 4 leptons mass int ntagtight=0; int ntagloose=0; Jet* jetw0=0; Jet* jetw1=0; //count number of b tags by looping over all other jets in the event //These are the default jets at the hadron scale JetCollection::const_iterator jetItr2 (m_jets->begin()), jetItr2E (m_jets->end()); // loop over the AOD jet container //These are ordered by jet weight for (; jetItr2 != jetItr2E; ++jetItr2) { Jet* jet2 = *jetItr2; // float deta=jet2->eta()-jet0->eta(); float dphi=acos(cos(jet2->phi()-jet0->phi())); float dr0=sqrt(deta*deta+dphi*dphi); deta=jet2->eta()-jet1->eta(); dphi=acos(cos(jet2->phi()-jet1->phi())); float dr1=sqrt(deta*deta+dphi*dphi); if(dr0<0.2||dr1<0.2){ // std:: cout <<" same jet found " <ptCal()minPt()) continue; float jetweight = !m_btaggerName.empty() ? jet2-> getFlavourTagWeight(m_btaggerName) : jet2-> getFlavourTagWeight(); if(jetweight>m_jetWeight0Cut) ntagtight++; if(jetweight>m_jetWeight1Cut) ntagloose++; //These are the 2 jets with the highest weights (besides the 2 selected as electrons) if(jetw0 ==0 ) jetw0=jet2; else if(jetw1 ==0 ) jetw1=jet2; } if(jetw0 ==0 || jetw1 ==0 ) { // std::cout <<"OSHbb::FillFakeHistos(): 4 jets but cannot find 2 tag jets. Something went wrong" <hlvCal(jetw0)), hlvjw1(m_selectJet->hlvCal(jetw1)); float dibmass =(hlvjw0+hlvjw1).m(); static const double M_Z=91.188*GeV; float M_bb = dibmass; float scale0 = 1 + (M_Z*M_Z - M_bb*M_bb) * hlvjw1.e() / (2*M_bb*M_bb*(hlvjw0.e()+hlvjw1.e())) ; float scale1 = 1 + (M_Z*M_Z - M_bb*M_bb) * hlvjw0.e() / (2*M_bb*M_bb*(hlvjw0.e()+hlvjw1.e())) ; hlvjw0 *= scale0; hlvjw1 *= scale1; float HiggsZZMass=(jetEscale*jet0->hlv()+jetEscale*jet1->hlv()+hlvjw0+hlvjw1).m(); float JJDeltaPhi=(jetw0->hlv()).v().deltaPhi((jetw1->hlv()).v()); if(m_passCutEtmiss) { m_osHistSvc->fillTH1Flav(m_histDirZ+"MZDijetFake", zmass/GeV, m_event_weight*zfakeweight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirZ+"MZeeDijetFake", zmass/GeV, m_event_weight*zfakeweight,m_event_flav); } if(m_passCutEtmiss&&m_passCutZMass) { if(m_event_weight*zfakeweight>0.5) std::cout <<"Something wrong with weights" << " m_event_weight " <fillTH1Flav(m_histDirJets+"DijetMassFake", dibmass/GeV,m_event_weight*zfakeweight,flavuse); if (dibmass>m_zbmasslow&&dibmassfillTH1Flav(m_histDirJets+"DijetZMassFake", HiggsZZMass/GeV,m_event_weight*zfakeweight,flavuse); if(fabs(JJDeltaPhi)fillTH1Flav(m_histDirJets+"DijetZMassDPhiCutsFake", HiggsZZMass/GeV,m_event_weight*zfakeweight,flavuse); } if(dibmass>m_zbmasslowsb&&dibmassfillTH1Flav(m_histDirJets+"DijetZMassSBLowFake", HiggsZZMass/GeV,m_event_weight*zfakeweight,flavuse); if(fabs(JJDeltaPhi)fillTH1Flav(m_histDirJets+"DijetZMassSBLowDPhiCutsFake", HiggsZZMass/GeV,m_event_weight*zfakeweight,flavuse); } else if(dibmass>m_zbmasshigh&&dibmassfillTH1Flav(m_histDirJets+"DijetZMassSBHighFake", HiggsZZMass/GeV,m_event_weight*zfakeweight,flavuse); if(fabs(JJDeltaPhi)fillTH1Flav(m_histDirJets+"DijetZMassSBHighDPhiCutsFake", HiggsZZMass/GeV,m_event_weight*zfakeweight,flavuse); } } if((zmass>m_zmasslowsb &&zmassm_zmasshigh &&zmassfillTH1Flav(m_histDirZSB+"DijetMassZSBLowMETFake", dibmass/GeV,m_event_weight*zfakeweight,flavuse); if (dibmass>m_zbmasslow&&dibmassfillTH1Flav(m_histDirZSB+"DijetZMassZSBLowMETFake", HiggsZZMass/GeV,m_event_weight*zfakeweight,flavuse); } else { m_osHistSvc->fillTH1Flav(m_histDirZSB+"DijetMassZSBHighMETFake", dibmass/GeV,m_event_weight*zfakeweight,flavuse); if (dibmass>m_zbmasslow&&dibmassfillTH1Flav(m_histDirZSB+"DijetZMassZSBHighMETFake", HiggsZZMass/GeV,m_event_weight*zfakeweight,flavuse); } } if(ntagtight>0&&ntagloose>1) { if(m_passCutEtmiss) { m_osHistSvc->fillTH1Flav(m_histDirZ+"MZDijetWCFake", zmass/GeV, m_event_weight*zfakeweight,m_event_flav); m_osHistSvc->fillTH1Flav(m_histDirZ+"MZeeDijetWCFake", zmass/GeV, m_event_weight*zfakeweight,m_event_flav); } if(m_passCutEtmiss&&m_passCutZMass) { m_osHistSvc->fillTH1Flav(m_histDirJets+"DijetMassWCFake", dibmass/GeV,m_event_weight*zfakeweight,flavuse); if(dibmass>m_zbmasslow&&dibmassfillTH1Flav(m_histDirJets+"DijetZMassWCFake", HiggsZZMass/GeV,m_event_weight*zfakeweight,flavuse); if(fabs(JJDeltaPhi)fillTH1Flav(m_histDirJets+"DijetZMassDPhiCutsWCFake", HiggsZZMass/GeV,m_event_weight*zfakeweight,flavuse); } if(dibmass>m_zbmasslowsb&&dibmassfillTH1Flav(m_histDirJets+"DijetZMassSBLowWCFake", HiggsZZMass/GeV,m_event_weight*zfakeweight,flavuse); if(fabs(JJDeltaPhi)fillTH1Flav(m_histDirJets+"DijetZMassSBLowDPhiCutsWCFake", HiggsZZMass/GeV,m_event_weight*zfakeweight,flavuse); } else if(dibmass>m_zbmasshigh&&dibmassfillTH1Flav(m_histDirJets+"DijetZMassSBHighWCFake", HiggsZZMass/GeV,m_event_weight*zfakeweight,flavuse); if(fabs(JJDeltaPhi)fillTH1Flav(m_histDirJets+"DijetZMassSBHighDPhiCutsFake", HiggsZZMass/GeV,m_event_weight*zfakeweight,flavuse); } } if((zmass>m_zmasslowsb &&zmassm_zmasshigh &&zmassfillTH1Flav(m_histDirZSB+"DijetMassZSBLowMETWCFake", dibmass/GeV,m_event_weight*zfakeweight,flavuse); if (dibmass>m_zbmasslow&&dibmassfillTH1Flav(m_histDirZSB+"DijetZMassZSBLowMETWCFake", HiggsZZMass/GeV,m_event_weight*zfakeweight,flavuse); } else { m_osHistSvc->fillTH1Flav(m_histDirZSB+"DijetMassZSBHighMETWCFake", dibmass/GeV,m_event_weight*zfakeweight,flavuse); if (dibmass>m_zbmasslow&&dibmassfillTH1Flav(m_histDirZSB+"DijetZMassZSBHighMETWCFake", HiggsZZMass/GeV,m_event_weight*zfakeweight,flavuse); } } } } } } float OSHbb::GetEleFakeRate(float ptj) { static const float p0 = 2.64727e-03; static const float p1 = 2.62354e+04; static const float p2 = 1.97190e+03; // //From data // static const float p0 =1.72217e-03 ; // static const float p1 =2.78632e+04 ; // static const float p2 =3.06408e+03 ; return p0/(1+exp(-(ptj-p1)/p2)); } // float OSHbb::GetMuFakeRate(float ptj) { // return 0.0007; // // } */ void OSHbb::ScaleDijetsUsingMETConstraint(HepLorentzVector& hlvJ0,HepLorentzVector& hlvJ1) { // Using the MET constraint to scale the two jets if (m_selectJet->numSelected() != 2) return; float scale0 = 1 + (m_EtMiss_y*std::cos(hlvJ1.phi()) - m_EtMiss_x*std::sin(hlvJ1.phi())) / (hlvJ0.et()*std::sin(hlvJ0.phi() - hlvJ1.phi())); float scale1 = 1 + (m_EtMiss_y*std::cos(hlvJ0.phi()) - m_EtMiss_x*std::sin(hlvJ0.phi())) / (hlvJ1.et()*std::sin(hlvJ1.phi() - hlvJ0.phi())); hlvJ0 *= scale0; hlvJ1 *= scale1; } void OSHbb::ScaleDijetsUsingZConstraint(HepLorentzVector& hlvJ0ZC,HepLorentzVector& hlvJ1ZC) { bool simplescale=true; static const double M_Z=91.188*GeV; float M_bb = (hlvJ0ZC + hlvJ1ZC).m(); float scale0,scale1; if(simplescale){ scale0=scale1=M_Z/M_bb; } else { scale0 = 1 + (M_Z*M_Z - M_bb*M_bb) * hlvJ1ZC.e() / (2*M_bb*M_bb*(hlvJ0ZC.e()+hlvJ1ZC.e())) ; scale1 = 1 + (M_Z*M_Z - M_bb*M_bb) * hlvJ0ZC.e() / (2*M_bb*M_bb*(hlvJ0ZC.e()+hlvJ1ZC.e())) ; } hlvJ0ZC *= scale0; hlvJ1ZC *= scale1; } void OSHbb::FindTruthParticles() { //TODO: modify for Znunu and WH analysis. //Clear pointers from last event mc_muPlus = 0; mc_muMinus = 0; mc_ePlus = 0; mc_eMinus = 0; mc_b = 0; mc_bbar = 0; mc_Zll = 0; mc_Zbb = 0; mc_H = 0; m_bCentralEvent=0; m_cCentralEvent=0; m_semiLeptonicBDecay=0; m_semiLeptonicCDecay=0; m_chronoSvc->chronoStart("truthselect"); Check( m_selectTruth->select() );//cuts out geant/documentary particles m_chronoSvc->chronoStop("truthselect"); // For pythia, status = 1 for stable particles, 2 for unstable particles //Get electrons - results include particle & anti-particle TruthParticleContainer* truE = m_selectTruth->selectedTruthParticlesPDG(PDG::e_minus,1); //Get muons - results include particle & anti-particle TruthParticleContainer* truMu = m_selectTruth->selectedTruthParticlesPDG(PDG::mu_minus,1); //Get jets TruthParticleContainer* truB = m_selectTruth->selectedTruthParticlesPDG(PDG::b,2); TruthParticleContainer* truC = m_selectTruth->selectedTruthParticlesPDG(PDG::c,2); //Get Z's TruthParticleContainer* truZ = m_selectTruth->selectedTruthParticlesPDG(PDG::Z0,2); //Get W's TruthParticleContainer* truW = m_selectTruth->selectedTruthParticlesPDG(PDG::W_plus,2); //Get Nu's TruthParticleContainer* truNu = m_selectTruth->selectedTruthParticlesPDG(PDG::nu_e,1); TruthParticleContainer* truNuMu = m_selectTruth->selectedTruthParticlesPDG(PDG::nu_mu,1); std::copy(truNuMu->begin(), truNuMu->end(), std::back_inserter(*truNu)); delete truNuMu; //set up a container of truth leptons m_emufromwztruth= new TruthParticleContainer(SG::VIEW_ELEMENTS); m_emufromwztruth->reserve(10); TruthParticleContainer::const_iterator it; //Assign truth electons for(it = truE->begin(); it < truE->end(); it++) { TruthParticle* mcpart = *it; TruthParticle* myZ = 0; if( m_selectTruth->IsFrom(mcpart,PDG::Z0,myZ,0) && myZ) { if( (*it)->pdgId() > 0 ) mc_eMinus = mcpart; else mc_ePlus = mcpart; m_emufromwztruth->push_back(mcpart); m_osHistSvc->fillTH1Flav(m_histDirTruth+"PtLepfromZGen", mcpart->pt()/GeV, m_event_weight); } //Add e,mus from W to array TruthParticle* myW = 0; if( m_selectTruth->IsFrom(mcpart,PDG::W_plus,myW,0) && myW ) m_emufromwztruth->push_back(mcpart); } //Assign truth muons for(it = truMu->begin(); it < truMu->end(); it++) { TruthParticle* mcpart = *it; TruthParticle* myZ = 0; if( m_selectTruth->IsFrom(mcpart,PDG::Z0,myZ,0) && myZ ) { if( (mcpart)->pdgId() > 0 ) mc_muMinus = mcpart; else mc_muPlus = mcpart; m_emufromwztruth->push_back(mcpart); m_osHistSvc->fillTH1Flav(m_histDirTruth+"PtLepfromZGen", mcpart->pt()/GeV, m_event_weight); } //Add e,mus from W to array TruthParticle* myW = 0; if( m_selectTruth->IsFrom(mcpart,PDG::W_plus,myW,0) && myW ) m_emufromwztruth->push_back(mcpart); } //Assign truth b for(it = truB->begin(); it < truB->end(); it++) { TruthParticle* mcpart = *it; if(fabs(mcpart->eta())<3&&mcpart->pt()>10) m_bCentralEvent++; TruthParticle* myZ = 0; //B's can also be direct from the H or Z depending on signal (ZH,H->ZZ). if( (m_selectTruth->IsFrom(mcpart,PDG::Z0,myZ,0) || m_selectTruth->IsFrom(mcpart,PDG::Higgs0,myZ,0)) && myZ ) { if( (mcpart)->pdgId() > 0 ) mc_b = mcpart; else mc_bbar = mcpart; } } //Assign truth c for(it = truC->begin(); it < truC->end(); it++) { TruthParticle* mcpart = *it; if(fabs(mcpart->eta())<3&&mcpart->pt()>10) m_cCentralEvent++; } if(m_bCentralEvent>1) m_event_flav = 5; // b-jet else if( m_cCentralEvent>1) m_event_flav = 4; // c-jet // for(it = truNu->begin(); it < truNu->end(); it++) { TruthParticle* mcpart = *it; TruthParticle* myNu = 0; if (m_selectTruth->IsFrom(mcpart,PDG::b,myNu,true) && myNu) ++m_semiLeptonicBDecay; if (m_selectTruth->IsFrom(mcpart,PDG::c,myNu,true) && myNu) ++m_semiLeptonicCDecay; } //Assign truth Z for(it = truZ->begin(); it < truZ->end(); it++) { TruthParticle* mcpart = *it; if( mcpart->hasChild(PDG::b) ) { mc_Zbb = mcpart; } else if( mcpart->hasChild(PDG::e_minus)) { mc_Zll = mcpart; mc_Zee = mcpart; } else if ( mcpart->hasChild(PDG::mu_minus)){ mc_Zll = mcpart; mc_Zmumu = mcpart; } m_osHistSvc->fillTH1Flav(m_histDirTruth+"PtZGen", mcpart->pt()/GeV, m_event_weight, m_event_flav); } //Assign truth W for(it = truW->begin(); it < truW->end(); it++) { TruthParticle* mcpart = *it; m_osHistSvc->fillTH1Flav(m_histDirTruth+"PtZGen", mcpart->pt()/GeV, m_event_weight, m_event_flav); } m_zbzlgen=mc_Zbb>0&&mc_Zll>0; //Get Higgs TruthParticleContainer* truH = m_selectTruth->selectedTruthParticlesPDG(PDG::Higgs0,2); //Assign truth Higgs for(it = truH->begin(); it < truH->end(); it++) mc_H = *it; if(mc_H) m_osHistSvc->fillTH1Flav(m_histDirTruth+"TruthHiggsMass", mc_H->m()/GeV, m_event_weight, m_event_flav); // generated jets at the hadron level m_chronoSvc->chronoStart("hadronjetselect"); if (m_removeHadOverlapWZDecays) m_selectJet->setTruthParticlesForOverlap(m_emufromwztruth); if (m_removeHadOverlapWZDecays) m_selectJetlowpt->setTruthParticlesForOverlap(m_emufromwztruth); Check( m_selectJet->selectHad()); Check( m_selectJetlowpt->selectHad()); m_chronoSvc->chronoStop("hadronjetselect"); if (m_doPartonJets) { m_chronoSvc->chronoStart("partonjetselect"); Check( m_selectJet->selectPart()); m_chronoSvc->chronoStop("partonjetselect"); } return; } TruthParticle* OSHbb::FindTruthMatch(IParticle* recpart, PDG::pidType inpdg, int status, double deltaRcone) { TruthParticleContainer* tru = m_selectTruth->selectedTruthParticlesPDG(inpdg,status); TruthParticle* p = 0; for(TruthParticleContainer::iterator it = tru->begin(); it < tru->end(); ++it) { if(P4Helpers::deltaR(recpart,*it)pt()-p->pt()) < fabs(recpart->pt()-(*it)->pt()) ) ? p : *it; } else { p=*it; } } } return p; } bool OSHbb::SortBWeight(const Jet* a,const Jet* b) { return (m_btaggerName.compare("") != 0) ? a->getFlavourTagWeight(m_btaggerName) > b->getFlavourTagWeight(m_btaggerName) : a->getFlavourTagWeight() > b->getFlavourTagWeight(); } void OSHbb::Check(const StatusCode& sc) { if(!sc.isSuccess()) { msg(MSG::ERROR) << "StatusCode check failed! Printing stack trace." << endreq; Athena::DebugAids::stacktrace(); } return; } void OSHbb::Check(const StatusCode& sc,const std::string& s) { if(!sc.isSuccess()) { msg(MSG::ERROR) << "StatusCode check failed :" << s << endreq; } return; } float OSHbb::AngleBetweenLeptons() { msg(MSG::DEBUG) << "AngleBetweenLeptons()" << endreq; if(m_selectElec->numSelected()>1) { return m_electrons->at(0)->hlv().angle(m_electrons->at(1)->hlv()); } else if(m_selectMuon->numSelected()>1) { return m_muons->at(0)->hlv().angle(m_muons->at(1)->hlv()); } else return -1.0; } void OSHbb::FillSystematics(){ // shifts up, down and resol msg(MSG::DEBUG) << "Systematics()" << endreq; // Only apply systematics if they have been requested by all selectors // if (!m_selectElec->getUseSys() || !m_selectMuon->getUseSys() || !m_selectJet->getUseSys() ) return; for(int isys=0; isyssetSys(OSSelectorsEnums::eSysNone); m_selectElec->setSys(OSSelectorsEnums::eSysNone); m_selectMuon->setSys(OSSelectorsEnums::eSysNone); if(isys==0) { m_selectJet->setSys(OSSelectorsEnums::eSysEnergyUp); } else if(isys==1) { m_selectJet->setSys(OSSelectorsEnums::eSysEnergyDo); } else if(isys==2) { m_selectJet->setSys(OSSelectorsEnums::eSysEnergyResol); } else if(isys==3) { m_selectElec->setSys(OSSelectorsEnums::eSysEnergyUp); } else if(isys==4) { m_selectElec->setSys(OSSelectorsEnums::eSysEnergyDo); } else if(isys==5) { m_selectElec->setSys(OSSelectorsEnums::eSysEnergyResol); } else if(isys==6) { m_selectMuon->setSys(OSSelectorsEnums::eSysEnergyUp); } else if(isys==7) { m_selectMuon->setSys(OSSelectorsEnums::eSysEnergyDo); } else if(isys==8) { m_selectMuon->setSys(OSSelectorsEnums::eSysEnergyResolIDUp); } else if(isys==9) { m_selectMuon->setSys(OSSelectorsEnums::eSysEnergyResolIDDo); } else if(isys==10) { m_selectMuon->setSys(OSSelectorsEnums::eSysEnergyResolMSUp); } else if(isys==11) { m_selectMuon->setSys(OSSelectorsEnums::eSysEnergyResolMSDo); } else if(isys==12) { m_selectMuon->setSys(OSSelectorsEnums::eSysEnergyResolNoSmear); } else if(isys==13) { m_metsyshadshiftdo=false; m_metsyshadshiftup=true; m_metsyshadsmear=false; } else if(isys==14) { m_metsyshadshiftdo=true; m_metsyshadshiftup=false; m_metsyshadsmear=false; } else if(isys==15) { m_metsyshadshiftdo=false; m_metsyshadshiftup=false; m_metsyshadsmear=true; } m_electrons = m_selectElec->selectedElecs(); m_muons = m_selectMuon->selectedMuons(); m_jets = m_selectJet->selectedJets(); msg(MSG::DEBUG) << "sys findleptons()" << endreq; FindLeptons(); bool passCutsElectron = ApplyElectronCuts(false); bool passCutsMuon = ApplyMuonCuts(false); if(!passCutsElectron&&!passCutsMuon) continue; if(!ZLeptonMassCut() ) continue; if(!m_passCutZMass) continue; //For W reject events if they have 2 leptons of different types if(m_VectorBosonOp==eW&&m_muons->size()>0&&m_electrons->size()>0) continue; //For Z reject events if they have 2 lepton pairs of different types if(m_VectorBosonOp==eZ&&m_passCutsMuon&&m_passCutsElectron) continue; MakeDPhiLL(); CalculateMissingEt(); if(!METCuts(false)) continue; //Cut on number of jets //if(!JetCleaningCut(false)) continue; if(!JetVertexFractionCut(false)) continue; FindJets(); if(!JetCountCut(false)) continue; //Cut on 3rd Jet Weight m_passCutThirdJetWeight = ThirdJetWeightCut(false); if (m_VectorBosonOp==eW && !m_passCutThirdJetWeight) return; if(m_VectorBoson2Op==eZnu&&!FirstJetWeightCut(false)) continue; if(!BJetsBackToBackCut(false)) continue; MakeJJMass(); m_HiggsZZMass=ZZHiggsMass(m_bjet0,m_bjet1,1); m_HiggsZZMassNoScale=ZZHiggsMass(m_bjet0,m_bjet1,1); //Higgs mass based on highest Pt jets for 1 or 0 btags m_HiggsZZMassOPt=ZZHiggsMass(m_jet0,m_jet1,0); m_HiggsZZMassOPtNoScale=ZZHiggsMass(m_jet0,m_jet1,0); FillFinalHistos( m_sysNames[isys]); } //Set all systematics to standard m_selectJet->setSys(OSSelectorsEnums::eSysNone); m_selectElec->setSys(OSSelectorsEnums::eSysNone); m_selectMuon->setSys(OSSelectorsEnums::eSysNone); m_electrons = m_selectElec->selectedElecs(); m_muons = m_selectMuon->selectedMuons(); m_jets = m_selectJet->selectedJets(); m_metsyshadshiftdo=false; m_metsyshadshiftup=false; m_metsyshadsmear=false; m_metsysfwhadshiftdo=false; m_metsysfwhadshiftup=false; m_metsysfwhadsmear=false; } void OSHbb::FillFinalHistos(std::string tag){ if(!m_passCutsElectron&&!m_passCutsMuon) return; for (int i=0; i<2; i++) { std::string pre= ""; if(i==1) pre= (m_passCutsElectron)? "Elec" : "Muon"; m_osHistSvc->fillTH1Flav(m_histDirZ+pre+"ZMass"+tag, m_mZ/GeV,m_event_weight,m_event_flav); if(m_passCutZMass) m_osHistSvc->fillTH1Flav(m_histDirMET+pre+"MissingET"+tag,m_EtMiss/GeV, m_event_weight,m_event_flav); //VH Histos and Mjj FillFinalVHHistos(pre,tag); //ZZ histos FillFinalHZZHistos(pre,tag); FillFinalHZZHistos0or1WC(pre,tag); } } void OSHbb::FillFinalVHHistos(std::string& pre, std::string& tag){ if(m_jet0 ==0 || m_jet1 ==0) return; int flavuse = m_flav0==m_flav1 ? m_flav0 : 0; int nrecbjet=NBJets(); m_osHistSvc->fillTH1Flav(m_histDirJets+pre+"DijetMass"+tag, m_dibmassopt/GeV,m_event_weight,flavuse); if(nrecbjet==0 && m_passCutThirdJetWeight) m_osHistSvc->fillTH1Flav(m_histDirJets+pre+"DijetMass0WC"+tag, m_dibmassopt/GeV,m_event_weight,flavuse); else if (nrecbjet==1 && m_passCutThirdJetWeight) m_osHistSvc->fillTH1Flav(m_histDirJets+pre+"DijetMass1WC"+tag, m_dibmassopt/GeV,m_event_weight,flavuse); // Nbtag jets after dijet cut m_osHistSvc->fillTH1Flav(m_histDirJets+"NJetBTag", nrecbjet, m_event_weight, flavuse); //Cut based analysis if(m_cutsloose) m_osHistSvc->fillTH1Flav(m_histDirJets+pre+"DijetMassCutsLoose"+tag, m_dibmassopt/GeV,m_event_weight,flavuse); if(m_cutsmedium) m_osHistSvc->fillTH1Flav(m_histDirJets+pre+"DijetMassCutsMedium"+tag, m_dibmassopt/GeV,m_event_weight,flavuse); if(m_cutstight) m_osHistSvc->fillTH1Flav(m_histDirJets+pre+"DijetMassCutsTight"+tag, m_dibmassopt/GeV,m_event_weight,flavuse); if( nrecbjet>1 && m_passCutThirdJetWeight) { m_osHistSvc->fillTH1Flav(m_histDirJets+pre+"DijetMassWC"+tag, m_dibmass/GeV,m_event_weight,flavuse); //Cut based analysis if(m_cutsloose) m_osHistSvc->fillTH1Flav(m_histDirJets+pre+"DijetMassCutsLooseWC"+tag, m_dibmass/GeV,m_event_weight,flavuse); if(m_cutsmedium) m_osHistSvc->fillTH1Flav(m_histDirJets+pre+"DijetMassCutsMediumWC"+tag, m_dibmass/GeV,m_event_weight,flavuse); if(m_cutstight) m_osHistSvc->fillTH1Flav(m_histDirJets+pre+"DijetMassCutsTightWC"+tag, m_dibmass/GeV,m_event_weight,flavuse); } } void OSHbb::FillFinalHZZHistos(std::string& pre, std::string& tag){ if(m_jet0 ==0 || m_jet1 ==0) return; //if(!ZBJetMassCut()) return; // WRONG for the llqq analysis where should cut in the highest Pt jets not the 2 highest b weights int flavuse = m_flav0==m_flav1 ? m_flav0 : 0; if (ZBJetMassCutOPt()) { m_osHistSvc->fillTH1Flav(m_histDirJets+pre+"DijetZMass"+tag, m_HiggsZZMassOPt/GeV,m_event_weight,flavuse); if(m_jets->size()>=4) m_osHistSvc->fillTH1Flav(m_histDirJets+pre+"DijetZMassNJ4"+tag, m_HiggsZZMassOPt/GeV,m_event_weight,flavuse); if(m_JJDeltaPhiOPtfillTH1Flav(m_histDirJets+pre+"DijetZMassDPhiCuts"+tag, m_HiggsZZMassOPt/GeV,m_event_weight,flavuse); if(m_jets->size()>=4) m_osHistSvc->fillTH1Flav(m_histDirJets+pre+"DijetZMassDPhiJJNJ4"+tag, m_HiggsZZMassOPt/GeV,m_event_weight,flavuse); } if(m_selectJet->hlvCal(m_jet0).perp()>50*GeV&&m_selectJet->hlvCal(m_jet1).perp()>50*GeV) { m_osHistSvc->fillTH1Flav(m_histDirJets+pre+"DijetZMassPtJ50"+tag, m_HiggsZZMassOPt/GeV,m_event_weight,flavuse); if(m_JJDeltaPhiOPtfillTH1Flav(m_histDirJets+pre+"DijetZMassDPhiJJPtJ50"+tag, m_HiggsZZMassOPt/GeV,m_event_weight,flavuse); // Final events for scanning if (tag.empty()) std::cout << "Final high Mh events: " << m_eventInfo->event_ID()->run_number() << " " << m_eventInfo->event_ID()->event_number() << " " << m_HiggsZZMassOPt/GeV << std::endl; } } if( ZBJetMassCut() && NBJets()>1 && m_passCutThirdJetWeight) { m_osHistSvc->fillTH1Flav(m_histDirJets+pre+"DijetZMassWC"+tag, m_HiggsZZMass/GeV,m_event_weight,flavuse); if(m_jets->size()>=4) m_osHistSvc->fillTH1Flav(m_histDirJets+pre+"DijetZMassNJ4WC"+tag, m_HiggsZZMass/GeV,m_event_weight,flavuse); if(m_JJDeltaPhifillTH1Flav(m_histDirJets+pre+"DijetZMassDPhiCutsWC"+tag, m_HiggsZZMass/GeV,m_event_weight,flavuse); if(m_jets->size()>=4) m_osHistSvc->fillTH1Flav(m_histDirJets+pre+"DijetZMassDPhiJJNJ4WC"+tag, m_HiggsZZMass/GeV,m_event_weight,flavuse); } if(m_selectJet->hlvCal(m_bjet0).perp()>50*GeV&&m_selectJet->hlvCal(m_bjet1).perp()>50*GeV) { m_osHistSvc->fillTH1Flav(m_histDirJets+pre+"DijetZMassPtJ50WC"+tag, m_HiggsZZMass/GeV,m_event_weight,flavuse); if(m_JJDeltaPhifillTH1Flav(m_histDirJets+pre+"DijetZMassDPhiJJPtJ50WC"+tag, m_HiggsZZMass/GeV,m_event_weight,flavuse); } } } void OSHbb::FillFinalHZZHistos0or1WC(std::string& pre, std::string& tag){ if(m_jet0 ==0 || m_jet1 ==0) return; int flavuse = m_flav0==m_flav1 ? m_flav0 : 0; //ZZ histos if(!ZBJetMassCutOPt()) return; if(NBJets()==0 && m_passCutThirdJetWeight) { m_osHistSvc->fillTH1Flav(m_histDirJets+pre+"DijetZMass0WC"+tag, m_HiggsZZMassOPt/GeV,m_event_weight,flavuse); if(m_selectJet->hlvCal(m_jet0).perp()>50*GeV&&m_selectJet->hlvCal(m_jet1).perp()>50*GeV) { m_osHistSvc->fillTH1Flav(m_histDirJets+pre+"DijetZMassPtJ500WC"+tag, m_HiggsZZMassOPt/GeV,m_event_weight,flavuse); if(m_JJDeltaPhiOPtfillTH1Flav(m_histDirJets+pre+"DijetZMassDPhiJJPtJ500WC"+tag, m_HiggsZZMassOPt/GeV,m_event_weight,flavuse); } } else if(NBJets()==1 && m_passCutThirdJetWeight) { m_osHistSvc->fillTH1Flav(m_histDirJets+pre+"DijetZMass1WC"+tag, m_HiggsZZMass/GeV,m_event_weight,flavuse); if(m_selectJet->hlvCal(m_jet0).perp()>50*GeV&&m_selectJet->hlvCal(m_jet1).perp()>50*GeV) { m_osHistSvc->fillTH1Flav(m_histDirJets+pre+"DijetZMassPtJ501WC"+tag, m_HiggsZZMassOPt/GeV,m_event_weight,flavuse); if(m_JJDeltaPhiOPtfillTH1Flav(m_histDirJets+pre+"DijetZMassDPhiJJPtJ501WC"+tag, m_HiggsZZMassOPt/GeV,m_event_weight,flavuse); } } } bool OSHbb::NeutrinoWConstraint(TLorentzVector& lep, TLorentzVector& nu0, TLorentzVector& nu1){ //Work out PZ of neutrino based on W mass constraint float alp=lep.Px()*nu0.Px()+lep.Py()*nu0.Py(); float ptn2=pow(nu0.Pt(),2); float ctl=cos(lep.Theta()); static const float Mw=80.399*GeV; float bet=(Mw*Mw+2*alp)/(2*lep.P()); //These follow the terms of the quadratic eqn ax**2+bx+c=0; float a=(1-ctl*ctl); float b=-2*bet*ctl; float c=ptn2-bet*bet; if(a<0) std::cout <<" OSHbb::NeutrinoRec has problems in calculation a<0" <4*a*c) { //We have a solution float pznu0=(-b+sqrt(b*b-4*a*c))/(2*a); float pznu1=(-b-sqrt(b*b-4*a*c))/(2*a); //order solutions according to lowest pz nu if(fabs(pznu1)size()>1) m_LLDeltaPhi=acos(cos(m_selectElec->hlvCal(m_electrons->at(0)).phi()-m_selectElec->hlvCal(m_electrons->at(1)).phi())); else if(m_muons->size()>1) m_LLDeltaPhi=acos(cos(m_selectMuon->hlvCal(m_muons->at(0)).phi()-m_selectMuon->hlvCal(m_muons->at(1)).phi())); } }