NonMCBackground/src/select_fakes.cc

//
// System headers
//
#include <vector>                   // STL vector class
#include <iostream>                 // standard I/O
#include <iomanip>                  // functions to format standard I/O
#include <bitset>  
#include <map>  
using namespace std;

//
// root headers
//
#include <TFile.h>                  
#include <TTree.h>                  
#include <TChain.h>                 
#include <TBranch.h>                
#include <TClonesArray.h>           

//
// TMVA
//
#include "TMVA/Reader.h"
#include "TMVA/Tools.h"
#include "TMVA/Config.h"
#include "TMVA/MethodBDT.h"

//
// ntuple format headers
//
#include "EventHeader.h"
#include "PFMet.h"
#include "Electron.h"
#include "Muon.h"
#include "Vertex.h"
#include "PFCandidate.h"
#include "PFCandidateCol.h"
#include "PileupInfoCol.h"
#include "PileupEnergyDensity.h"
#include "MCParticle.h"
#include "MCEventInfo.h"
#include "TriggerMask.h"
#include "TriggerTable.h"
#include "Names.h"
#include "BaseMod.h"

//
// our headers
//
#include "ParseArgs.h"
#include "MuonSelection.h"
#include "ElectronSelection.h"
#include "IsolationSelection.h"
#include "SelectionFuncs.h"
#include "CommonDefs.h"
#include "Angles.h"
#include "AngleTuple.h"
#include "FOTuple.h"
#include "TriggerUtils.h"
#include "PassHLT.h"
#include "RunLumiRangeMap.h"
#include "fake_defs.h"
#include "Various.h"
#include "CommonDefs.h"

#ifndef CMSSW_BASE
#define CMSSW_BASE "../"
#endif

using namespace mithep;
//
// globals
//----------------------------------------------------------------------------
TH1D *hpu_2011, *hpu_2012; 
mithep::RunLumiRangeMap rlrm;
vector<SimpleLepton> failingLeptons ; // for fake estimation
vector<SimpleLepton> passingLeptons;      // for fake estimation
vector<unsigned> cutvec;
vector<vector<unsigned> > zcutvec;
vector<vector<unsigned> > zzcutvec;
map<unsigned,float>       evtrhoMap;
vector<string> cutstrs;
bool passes_HLT_MC;
vector<bool>   PFnoPUflag;;

//
// prototypes
//----------------------------------------------------------------------------
SelectionStatus analysisSelection(ControlFlags &ctrl, const mithep::Muon *mu, const mithep::Vertex *vtx, mithep::Array<mithep::PFCandidate> *pfArr,
                                  mithep::Array<mithep::PileupEnergyDensity>  *puDArr);
SelectionStatus analysisSelection(ControlFlags &ctrl, const mithep::Electron *ele, const mithep::Vertex *vtx, mithep::Array<mithep::PFCandidate> *pfArr,
                                  mithep::Array<mithep::PileupEnergyDensity>  *puDArr);
void initRunLumiRangeMap(ControlFlags &ctrl); 
float getMass(const mithep::Muon *zmu1, const mithep::Muon *zmu2, const mithep::Electron *zele1, const mithep::Electron *zele2);
mithep::MuonTools::EMuonEffectiveAreaTarget eraMu;
mithep::ElectronTools::EElectronEffectiveAreaTarget eraEle;
vector<const mithep::PFCandidate*> photonsToVeto;
vector<SimpleLepton> findFakes(ControlFlags &ctrl,
                               mithep::Array<mithep::Muon> *muonArr,
                               mithep::Array<mithep::Electron> *electronArr,
                               mithep::Array<mithep::PFJet> *jetArr,
                               mithep::Array<mithep::PFCandidate> *pfArr,
                               mithep::Array<mithep::PileupInfo>  *puArr,
                               mithep::Array<mithep::PileupEnergyDensity>  *puDArr,
                               mithep::Array<mithep::Vertex>               *vtxArr,
                               const mithep::Vertex         *vtx,
                               const mithep::Muon *zmu1,
                               const mithep::Muon *zmu2,
                               const mithep::Electron *zele1,
                               const mithep::Electron *zele2);
SimpleLepton fillSimpleLepton(const Muon *mu, const Electron *ele, SelectionStatus *ctrl);
//=== MAIN =================================================================================================
int main(int argc, char** argv) 
{

  string cmsswpath(CMSSW_BASE);
  cmsswpath.append("/src");
  std::bitset<1024> triggerBits;
  vector<vector<string> > inputFiles;
  inputFiles.push_back(vector<string>());
  map<string,unsigned> entrymap; // number of unskimmed entries in each file

  //
  // args
  //--------------------------------------------------------------------------------------------------------------
  ControlFlags ctrl;
  parse_args( argc, argv, ctrl );
  ctrl.dump();

  // // hack: store low and high mass values in mcfmfname
  // float massLo,massHi;
  // TString massWindow(ctrl.mcfmfname);
  // massWindow.ReplaceAll("..."," ");
  // string line(massWindow);
  // stringstream ss(line);
  // ss >> massLo >> massHi;
  // assert(massLo>=0 && massLo <= 8000);
  // assert(massHi>=0 && massHi <= 8000);

  //
  // File I/O
  //--------------------------------------------------------------------------------------------------------------
  TChain * chain = new TChain("Events");
  TChain * hltchain = new TChain("HLT");

  string fname;
  unsigned total_unskimmed=0;
  if( !ctrl.inputfiles.empty() ) { 
    ifstream f(ctrl.inputfiles.c_str());
    while (f >> fname) { 
      if( !(strncmp( fname.c_str(), "#", 1 ) ) ) continue; // skip commented lines
      cout << "adding inputfile : " << fname.c_str() << endl;
      entrymap[string(fname.c_str())] = unskimmedEntries(fname.c_str());
      total_unskimmed += entrymap[string(fname.c_str())];
      chain->AddFile(fname.c_str());
      hltchain->AddFile(fname.c_str());
    }
  } else { 
    cout << "adding inputfile : " << ctrl.inputfile.c_str() << endl;
    unsigned unsk_ents = unskimmedEntries(ctrl.inputfile.c_str());
    entrymap[string(ctrl.inputfile.c_str())] = unsk_ents;
    total_unskimmed += unsk_ents;
    chain->AddFile(ctrl.inputfile.c_str());
    hltchain->AddFile(ctrl.inputfile.c_str());
  }
  // write the total number of unskimmed events that went into making this output file to a text file
  writeEntries(ctrl,total_unskimmed);

  UInt_t nEvents=0,nMuDenom=0,nMuNumer=0,nEleDenom=0,nEleNumer=0;

  // move this to compute_fakes
  // double metMax = 25;

  //
  // Setup
  //--------------------------------------------------------------------------------------------------------------
  const char * ofname;
  if( strcmp( ctrl.outputfile.c_str(), "") ) {
    ofname = ctrl.outputfile.c_str();
  } else {
    ofname = "fo.root";
  }

  // TFile *outFile = new TFile(ofname, "RECREATE");

  Angles angles;
  KinematicsStruct kinematics;
  InfoStruct evtinfo;
  // GenInfoStruct geninfo;

  AngleTuple nt( (const char*)ofname, (const char*)"zznt");
  nt.makeAngleBranch(angles);
  nt.makeKinematicsBranch(kinematics);
  nt.makeInfoBranch(evtinfo);

  vector<SimpleLepton> failingLeptons,passingLeptons;
  FOTuple foTree( nt.getFile(), (const char*)"FOtree");
  foTree.makeFailedLeptonBranch(failingLeptons);
  foTree.makePassedLeptonBranch(passingLeptons);

  NtuplerBase foFlatTree( nt.getFile(), (const char*)"FO");
  const char * brname = "FO";
  const char * varlist = VAR_LIST_FOINFO;
  foinfo foi;
  foFlatTree.makeBranch( brname, (void*) &foi, varlist); 

  //  initMuonIDMVA();
  // getEATargets(ctrl,eraMu,eraEle);
  initElectronIDMVAV1();
  initMuonIsoMVA();
  initElectronIDMVA();
  initElectronIsoMVA();
  initTrigger();
  if(!ctrl.mc)
    initRunLumiRangeMap(ctrl);

  //
  // Access samples and fill histograms
  TFile *inputFile=0;
  TTree *eventTree=0;  
  string currentFile("");
   
  mithep::EventHeader *info    = new mithep::EventHeader();
  mithep::Array<mithep::PFMet>                *metArr        = new mithep::Array<mithep::PFMet>();
  mithep::Array<mithep::Electron>             *electronArr   = new mithep::Array<mithep::Electron>();
  mithep::Array<mithep::Muon>                 *muonArr       = new mithep::Array<mithep::Muon>();
  mithep::Array<mithep::Vertex>               *vtxArr        = new mithep::Array<mithep::Vertex>();
  mithep::Array<mithep::PFCandidate>          *pfArr         = new mithep::Array<mithep::PFCandidate>();
  mithep::Array<mithep::PFJet>                *jetArr        = new mithep::Array<mithep::PFJet>();
  mithep::Array<mithep::PileupInfo>           *puArr         = new mithep::Array<mithep::PileupInfo>();
  mithep::Array<mithep::PileupEnergyDensity>  *puDArr        = new mithep::Array<mithep::PileupEnergyDensity>();
  mithep::Array<mithep::Track>                *trkArr        = new mithep::Array<mithep::Track>();
  mithep::Array<mithep::MCParticle>           *mcArr         = new mithep::Array<mithep::MCParticle>();
  mithep::MCEventInfo                         *mcEvtInfo     = new mithep::MCEventInfo();
  mithep::TriggerMask                         *trigMask      = new mithep::TriggerMask();
  mithep::TriggerTable                        *hltTable      = new mithep::TriggerTable();
  vector<string>                              *hltTableStrings  = new vector<string>();

  TString fElectronName(Names::gkElectronBrn);
  TString fMuonName(Names::gkMuonBrn);
  TString fInfoName(Names::gkEvtHeaderBrn);
  TString fPFMetName("PFMet");
  TString fPrimVtxName(Names::gkPVBrn);
  TString fPFCandidateName(Names::gkPFCandidatesBrn);
  TString fPFJetName(Names::gkPFJetBrn);
  TString fPileupInfoName(Names::gkPileupInfoBrn);
  TString fPileupEnergyDensityName(Names::gkPileupEnergyDensityBrn);
  TString fTrackName(Names::gkTrackBrn);
  TString fMCParticleName(Names::gkMCPartBrn);
  TString fMCEvtInfoName(Names::gkMCEvtInfoBrn);
  TString fTriggerMaskName(Names::gkHltBitBrn);
  TString fTriggerTableName(Names::gkHltTableBrn);
  
  chain->SetBranchAddress(fInfoName,        &info);
  chain->SetBranchAddress(fPFMetName,       &metArr);
  chain->SetBranchAddress(fElectronName,    &electronArr);
  chain->SetBranchAddress(fMuonName,        &muonArr);
  chain->SetBranchAddress(fPrimVtxName,     &vtxArr);
  chain->SetBranchAddress(fPFCandidateName, &pfArr);
  chain->SetBranchAddress(fPFJetName,       &jetArr);
  if( ctrl.mc ) { 
    chain->SetBranchAddress(fPileupInfoName,  &puArr);
    chain->SetBranchAddress(fMCParticleName,  &mcArr);
    chain->SetBranchAddress(fMCEvtInfoName,  &mcEvtInfo);
  }
  chain->SetBranchAddress(fPileupEnergyDensityName, &puDArr);
  chain->SetBranchAddress(fTrackName, &trkArr);
  chain->SetBranchAddress(fTriggerMaskName, &trigMask);
  cout << "hlttable: " << fTriggerTableName << endl;
  hltchain->SetBranchAddress(fTriggerTableName, &hltTableStrings);

  const mithep::Vertex              *vtx;          // best primary vertex in the event

  //  ginfo = NULL;
  // if( ctrl.mc ) { chain->SetBranchAddress("Gen",        &ginfo);}

  // only 1 HLT table / file ???
  hltchain->GetEntry(0);

  int imax = chain->GetEntries();
  cout << "nEntries: " << imax << endl;
  for(UInt_t ientry=0; ientry<imax; ientry++) {
    chain->GetEntry(ientry);
    if(!(ientry%10000)) cout << "entry: " << ientry << "\tfrac: " << setw(15) << float(ientry)/chain->GetEntries() << endl;
    
    string fname = string(chain->GetFile()->GetEndpointUrl()->GetFile());
    if( ctrl.era == 0  )
      setEra( fname, ctrl );
    getEATargets(ctrl,eraMu,eraEle);

    // pfNoPU 
    PFnoPUflag.clear();
    int pfnopu_size = makePFnoPUArray( pfArr, PFnoPUflag, vtxArr );
    assert(pfnopu_size == pfArr->GetEntries());

    if(!(ctrl.noJSON) ) {
      RunLumiRangeMap::RunLumiPairType rl(info->RunNum(), info->LumiSec());      
      if( !(rlrm.HasRunLumi(rl)) ) continue;
    }
    
    nEvents++;

    bool goodVertex = setPV( ctrl, vtxArr, vtx );
    if(!goodVertex) continue;

    // move this to compute_fakes
    // TLorentzVector pfmet; pfmet.SetPxPyPzE(metArr->At(0)->Mex(),metArr->At(0)->Mey(),0,0);
    // if(pfmet.Pt() > metMax) continue;

    // good lepton selection
    failingLeptons.clear();
    passingLeptons.clear();
    vector<const mithep::Muon*> goodMuons;
    vector<const mithep::Electron*> goodElectrons;
    vector<SelectionStatus> muselv,eleselv;
    for(int imu=0; imu<muonArr->GetEntries(); imu++) {
      const mithep::Muon *mu = (*muonArr)[imu];
      SelectionStatus status = analysisSelection(ctrl,mu,vtx,pfArr,puDArr);
      if(status.loose()) {
        goodMuons.push_back(mu);
        muselv.push_back(status);
        passingLeptons.push_back(fillSimpleLepton(mu,&status));
      } else {
        failingLeptons.push_back(fillSimpleLepton(mu,&status));
      }
    }
    for(Int_t i=0; i<electronArr->GetEntries(); i++) {
      const mithep::Electron *ele = (*electronArr)[i];
      bool muonMatch=false;
      for(unsigned imu=0; imu<muonArr->GetEntries(); imu++) {
        const mithep::Muon *mu = (*muonArr)[imu];
        if(dr(mu,ele) < 0.05) muonMatch = true;
      }
      if(muonMatch) continue;
      SelectionStatus status = analysisSelection(ctrl,ele,vtx,pfArr,puDArr);
      if(status.loose()) {
        goodElectrons.push_back(ele);
        eleselv.push_back(status);
        passingLeptons.push_back(fillSimpleLepton(ele,&status));
      } else {
        failingLeptons.push_back(fillSimpleLepton(ele,&status));
      }
    }

    // look for z pairs
    int NZCandidates=0;
    float bestMass=-999;
    const mithep::Muon *zmu1=0,*zmu2=0;
    const mithep::Electron *zele1=0,*zele2=0;
    SelectionStatus *status1,*status2;
    for(unsigned imu=0; imu<goodMuons.size(); imu++) {
      for(unsigned jmu=imu+1; jmu<goodMuons.size(); jmu++) {
        if(goodMuons[imu]->Charge() == goodMuons[jmu]->Charge()) continue;
        float mass = getMass(goodMuons[imu],goodMuons[jmu],0,0);
        if( fabs(mass - Z_MASS) < fabs(bestMass - Z_MASS)) {
          bestMass = mass;
          zmu1 = goodMuons[imu];   status1 = &muselv[imu];
          zmu2 = goodMuons[jmu];   status2 = &muselv[jmu];
        }
      }
    }
    for(unsigned iele=0; iele<goodElectrons.size(); iele++) {
      for(unsigned jele=iele+1; jele<goodElectrons.size(); jele++) {
        if(goodElectrons[iele]->Charge() == goodElectrons[jele]->Charge()) continue;
        float mass = getMass(0,0,goodElectrons[iele],goodElectrons[jele]);
        if( fabs(mass - Z_MASS) < fabs(bestMass - Z_MASS)) {
          bestMass = mass;
          zele1 = goodElectrons[iele];   status1 = &eleselv[iele];
          zele2 = goodElectrons[jele];   status2 = &eleselv[jele];
          zmu1 = zmu2 = 0;
        }
      }
    }
    if(bestMass < 0) continue;
    assert( (zmu1 && zmu2) || (zele1 && zele2));
    assert( !(zmu1!=0 && zele1!=0) );
    assert( !(zmu2!=0 && zele2!=0) );

    EventData data;
    data.Z1leptons.push_back(fillSimpleLepton(zmu1,zele1,status1));
    data.Z1leptons.push_back(fillSimpleLepton(zmu2,zele2,status2));
    data.Z2leptons.push_back(fillSimpleLepton(zmu1,zele1,status1));
    data.Z2leptons.push_back(fillSimpleLepton(zmu2,zele2,status2));

    // loop over fakes
    vector<SimpleLepton> fakes = findFakes(ctrl,muonArr,electronArr,jetArr,pfArr,puArr,puDArr,vtxArr,vtx,zmu1,zmu2,zele1,zele2);
    //????????????????????????????????????????????????????????????????????????????????????????
    // if(fakes.size() > 1) continue;
    //????????????????????????????????????????????????????????????????????????????????????????
    for(unsigned ifake=0; ifake<fakes.size(); ifake++) {
      SimpleLepton fake = fakes[ifake];
      if(zmu1) {
        if(dr(fake,zmu1)  < 0.3) continue;
        if(dr(fake,zmu2)  < 0.3) continue;
      }      
      if(zele1) {
        if(dr(fake,zele1) < 0.3) continue;
        if(dr(fake,zele2) < 0.3) continue;
      }

      if(abs(fake.type) == 13) {
        nMuDenom++;
        if(fake.isLoose) nMuNumer++;
      } else {
        nEleDenom++;
        if(fake.isLoose) nEleNumer++;
      }

      // fill output TTree
      fillKinematics(data,kinematics);
      TLorentzVector pfmet; pfmet.SetPxPyPzE(metArr->At(0)->Mex(),metArr->At(0)->Mey(),0,0);
      fillEventInfo( info, pfmet, evtinfo, ctrl.mc ? getNPU(puArr) : 0, vtxArr->GetEntries());
      foTree.Fill();
      nt.Fill();

      foi.run    = info->RunNum();
      foi.evt    = info->EvtNum();
      foi.lumi   = info->LumiSec();
      foi.type   = fake.type;
      foi.npu    = ctrl.mc ? getNPU(puArr) : 0;
      foi.npv    = vtxArr->GetEntries(); // note: should maybe apply some selection here
      foi.fopt   = fake.vec.Pt();
      foi.foeta  = fake.vec.Eta();
      foi.fophi  = fake.vec.Phi();
      foi.fopass = fake.isLoose ? 1 : 0;
      foi.mass   = getMass(zmu1,zmu2,zele1,zele2);
      foi.l1pt   = zmu1 ? zmu1->Pt() : zele1->Pt();
      foi.l2pt   = zmu2 ? zmu2->Pt() : zele2->Pt();
      foi.met    = pfmet.Pt();
      foFlatTree.Fill();
    }
  }
  
  // outFile->Write();
  // outFile->Close();
  foTree.getFile()->cd();
  foTree.getTree()->Write();
  foFlatTree.getFile()->cd();
  foFlatTree.getTree()->Write();
  nt.WriteClose();
  
  //--------------------------------------------------------------------------------------------------------------
  // Summary print out
  //==============================================================================================================
  cout << endl;
  cout << "*" << endl;
  cout << "* SUMMARY" << endl;
  cout << "*--------------------------------------------------" << endl;
  cout << endl;

  cout << " >>>    No. of events processed: " << nEvents << endl;
  cout << " >>> No. of denominator muons: "   << nMuDenom << endl;
  cout << " >>>   No. of numerator muons: "   << nMuNumer << endl;     
  cout << " >>>                      ratio: " << float(nMuNumer)/nMuDenom << endl;     
  cout << " >>> No. of denominator eles:  "   << nEleDenom << endl;
  cout << " >>>   No. of numerator eles:  "   << nEleNumer << endl;     
  cout << " >>>                      ratio: " << float(nEleNumer)/nEleDenom << endl;     
  cout << endl;

  TString outputDir(ofname);
  outputDir = outputDir(0,outputDir.Last('/'));
  ofstream txtfile;
  char txtfname[100];
  TString txtname(ctrl.outputfile);
  txtname.ReplaceAll(".root",".txt");
  cout << "writing text to: " << txtname << endl;
  txtfile.open(txtname);
  txtfile << "*" << endl;
  txtfile << "* SUMMARY" << endl;
  txtfile << "*--------------------------------------------------" << endl;
  txtfile << endl;

  txtfile << " >>>    No. of events processed: " << nEvents << endl;
  txtfile << " >>> No. of denominator muons: "   << nMuDenom << endl;
  txtfile << " >>>   No. of numerator muons: "   << nMuNumer << endl;     
  txtfile << " >>>                      ratio: " << float(nMuNumer)/nMuDenom << endl;     
  txtfile << " >>> No. of denominator eles:  "   << nEleDenom << endl;
  txtfile << " >>>   No. of numerator eles:  "   << nEleNumer << endl;     
  txtfile << " >>>                      ratio: " << float(nEleNumer)/nEleDenom << endl;     
  txtfile << endl;  
  txtfile.close();
  
  cout << " <> Output saved in " << outputDir << "/" << endl;      
  cout << endl;
  
  cout << "done!" << endl;
}
//----------------------------------------------------------------------------
void initRunLumiRangeMap(ControlFlags &ctrl) 
//----------------------------------------------------------------------------
{
  /*
    rlrm.AddJSONFile(std::string("./data/Cert_136033-149442_7TeV_Apr21ReReco_Collisions10_JSON.txt"));
    //  rlrm.AddJSONFile(std::string("./data/Cert_160404-173244_7TeV_PromptReco_Collisions11_JSON_v2.txt"));
    rlrm.AddJSONFile(std::string("./data/Cert_160404-178078_7TeV_PromptReco_Collisions11_JSON.txt"));
    rlrm.AddJSONFile(std::string("./data/Cert_160404-163869_7TeV_May10ReReco_Collisions11_JSON_v3.txt"));  
    rlrm.AddJSONFile(std::string("./data/Cert_170249-172619_7TeV_ReReco5Aug_Collisions11_JSON.txt"));  
    // r11b
    rlrm.AddJSONFile(std::string("./data/Cert_160404-180252_7TeV_PromptReco_Collisions11_JSON.txt"));  
  */

  // 2012 only for now ...
  rlrm.AddJSONFile(string(ctrl.jsonFile));

};
//----------------------------------------------------------------------------------------
float getMass(const mithep::Muon *zmu1, const mithep::Muon *zmu2, const mithep::Electron *zele1, const mithep::Electron *zele2)
{
  TLorentzVector lep1,lep2;
  if(zmu1 && zmu2) {
    lep1.SetPtEtaPhiM(zmu1->Pt(),zmu1->Eta(),zmu1->Phi(),MUON_MASS);
    lep2.SetPtEtaPhiM(zmu2->Pt(),zmu2->Eta(),zmu2->Phi(),MUON_MASS);
  } else if(zele1 && zele2) {
    lep1.SetPtEtaPhiM(zele1->Pt(),zele1->Eta(),zele1->Phi(),ELECTRON_MASS);
    lep2.SetPtEtaPhiM(zele2->Pt(),zele2->Eta(),zele2->Phi(),ELECTRON_MASS);
  } else assert(0);

  return (lep1 + lep2).M();
}
//----------------------------------------------------------------------------------------
SelectionStatus analysisSelection(ControlFlags &ctrl, const mithep::Muon *mu, const mithep::Vertex *vtx, mithep::Array<mithep::PFCandidate> *pfArr,
                                  mithep::Array<mithep::PileupEnergyDensity>  *puDArr)
{
  SelectionStatus status;
  status |= muonReferencePreSelection(ctrl,mu,vtx,pfArr);
  status |= muonIDPFSelection(ctrl,mu,vtx,pfArr);
  status |= muonReferenceIsoSelection(ctrl,mu,vtx,pfArr,puDArr,eraMu,photonsToVeto);
  return status;
}
//----------------------------------------------------------------------------------------
SelectionStatus analysisSelection(ControlFlags &ctrl, const mithep::Electron *ele, const mithep::Vertex *vtx, mithep::Array<mithep::PFCandidate> *pfArr,
                                  mithep::Array<mithep::PileupEnergyDensity>  *puDArr)
{
  SelectionStatus status;
  status |= electronReferencePreSelection(ctrl,ele,vtx);
  status |= electronReferenceIDMVASelectionV1(ctrl,ele,vtx);
  status |= electronReferenceIsoSelection(ctrl,ele,vtx,pfArr,puDArr,eraEle,photonsToVeto);
  return status;
}
//----------------------------------------------------------------------------------------
vector<SimpleLepton> findFakes(ControlFlags &ctrl,
                               mithep::Array<mithep::Muon> *muonArr,
                               mithep::Array<mithep::Electron> *electronArr,
                               mithep::Array<mithep::PFJet> *jetArr,
                               mithep::Array<mithep::PFCandidate> *pfArr,
                               mithep::Array<mithep::PileupInfo>  *puArr,
                               mithep::Array<mithep::PileupEnergyDensity>  *puDArr,
                               mithep::Array<mithep::Vertex>               *vtxArr,
                               const mithep::Vertex         *vtx,
                               const mithep::Muon *zmu1,
                               const mithep::Muon *zmu2,
                               const mithep::Electron *zele1,
                               const mithep::Electron *zele2)
{
  vector<SimpleLepton> fakes;

  for(Int_t imu=0; imu<muonArr->GetEntries(); imu++) {  
    const mithep::Muon* mu = (*muonArr)[imu];

    if(mu==zmu1 || mu==zmu2)     continue;
      
    SelectionStatus denomsel;
    denomsel |= muonPreSelectionNoD0DzIP(ctrl,mu,vtx,pfArr);;
    if( !denomsel.passPre() ) continue;

    SelectionStatus status = analysisSelection(ctrl,mu,vtx,pfArr,puDArr);

    SimpleLepton fake;
    fake.vec.SetPtEtaPhiM(mu->Pt(),mu->Eta(),mu->Phi(),MUON_MASS);
    fake.type = 13;
    fake.isLoose = status.loose();
    
    fakes.push_back(fake);
  }
  for(Int_t iele=0; iele<electronArr->GetEntries(); iele++) {   
    const mithep::Electron* ele = (*electronArr)[iele];

    bool muonMatch=false;
    for(unsigned imu=0; imu<muonArr->GetEntries(); imu++) {
      const mithep::Muon *mu = (*muonArr)[imu];
      if(dr(mu,ele) < 0.05) muonMatch = true;
    }
    if(muonMatch) continue;

    if(ele==zele1 || ele==zele2)     continue;
      
    SelectionStatus denomsel;
    denomsel |= electronPreSelectionNoD0DzIP(ctrl,ele,vtx);
    if( !denomsel.passPre() ) continue;

    SelectionStatus status = analysisSelection(ctrl,ele,vtx,pfArr,puDArr);

    SimpleLepton fake;
    fake.vec.SetPtEtaPhiM(ele->Pt(),ele->Eta(),ele->Phi(),ELECTRON_MASS);
    fake.type = 11;
    fake.isLoose = status.loose();
    
    fakes.push_back(fake);
  }

  return fakes;
}
//----------------------------------------------------------------------------------------
SimpleLepton fillSimpleLepton(const Muon *mu, const Electron *ele, SelectionStatus *ctrl)
{
  if(mu)       return fillSimpleLepton(mu, ctrl);
  else if(ele) return fillSimpleLepton(ele,ctrl);
  else assert(0);
  return SimpleLepton();
}
Revision:	1.10
Committed:	Thu Jul 12 15:10:53 2012 UTC (12 years, 10 months ago) by dkralph
Content type:	text/plain
Branch:	MAIN
Changes since 1.9:	+11 -7 lines
Log Message:	Added a PlotHeader that my plot exe's use heavily, and cleaned up the plotters. Cleaned up my fake rate code so it almost all runs off of ntuples made by applyZPlusX in one of either ZPlusF or ZPlusFF mode, rather than different selectors. Added a plotter for their fakes, NonMCBackground/src/plotTheirFakes.cc Made some adjustments to fake selector to be in sync tih ZPlusX and vice versa. Removed all the json_spirit and RunLumiRangeMap stuff from Util -- it was only there because we were compiling without MitCommon, now that we pick up MitCommon, it just causes problems to have this in two different places.