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 "ChargedParticle.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 "SampleWeight.h"
#include "CommonDefs.h"

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

using namespace mithep;
//
// globals
//----------------------------------------------------------------------------
TH1D *hpu_2011, *hpu_2012; 
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
//----------------------------------------------------------------------------
bool passResKill(const Muon *zmu1, const Muon *zmu2, const Electron *zele1, const Electron *zele2,
                 const Muon *mu, const Electron *ele);
bool passGhost(vector<const Muon*> goodMuons, const Muon *mu);
SelectionStatus denominatorSelection(ControlFlags &ctrl, const Muon *mu, const Electron *ele, const Vertex *vtx, Array<PFCandidate> *pfArr);
SelectionStatus analysisSelection(ControlFlags &ctrl, const Muon *mu, const Vertex *vtx, Array<PFCandidate> *pfArr,
                                  Array<PileupEnergyDensity>  *puDArr);
SelectionStatus analysisSelection(ControlFlags &ctrl, const Electron *ele, const Vertex *vtx, Array<PFCandidate> *pfArr,
                                  Array<PileupEnergyDensity>  *puDArr);
void initRunLumiRangeMap(ControlFlags &ctrl); 
float getMass(const Muon *zmu1, const Muon *zmu2, const Electron *zele1, const Electron *zele2);
MuonTools::EMuonEffectiveAreaTarget eraMu;
ElectronTools::EElectronEffectiveAreaTarget eraEle;
vector<const PFCandidate*> photonsToVeto;
vector<SimpleLepton> findFakes(ControlFlags &ctrl,
                               Array<Muon> *muonArr,
                               Array<Electron> *electronArr,
                               Array<PFJet> *jetArr,
                               Array<PFCandidate> *pfArr,
                               Array<PileupInfo>  *puArr,
                               Array<PileupEnergyDensity>  *puDArr,
                               Array<Vertex>               *vtxArr,
                               const Vertex         *vtx,
                               const Muon *zmu1,
                               const Muon *zmu2,
                               const Electron *zele1,
                               const Electron *zele2);
SimpleLepton fillSimpleLepton(const Muon *mu, const Electron *ele, SelectionStatus *ctrl, const Vertex *vtx);
//=== 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();

  //
  // 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
      if(fname.find("/store/") != string::npos) fname = "root://eoscms//"+fname;
      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 { 
    if(ctrl.inputfile.find("/store/") != string::npos) ctrl.inputfile = "root://eoscms//"+ctrl.inputfile;
    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("");
   
  EventHeader *info    = new EventHeader();
  Array<PFMet>                *metArr        = new Array<PFMet>();
  Array<Electron>             *electronArr   = new Array<Electron>();
  Array<Muon>                 *muonArr       = new Array<Muon>();
  Array<Vertex>               *vtxArr        = new Array<Vertex>();
  Array<PFCandidate>          *pfArr         = new Array<PFCandidate>();
  Array<PFJet>                *jetArr        = new Array<PFJet>();
  Array<PileupInfo>           *puArr         = new Array<PileupInfo>();
  Array<PileupEnergyDensity>  *puDArr        = new Array<PileupEnergyDensity>();
  Array<Track>                *trkArr        = new Array<Track>();
  Array<MCParticle>           *mcArr         = new Array<MCParticle>();
  MCEventInfo                         *mcEvtInfo     = new MCEventInfo();
  TriggerMask                         *trigMask      = new TriggerMask();
  TriggerTable                        *hltTable      = new 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 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;

    // good lepton selection
    failingLeptons.clear();
    passingLeptons.clear();
    vector<const Muon*> goodMuons;
    vector<const Electron*> goodElectrons;
    vector<SelectionStatus> muselv,eleselv;
    for(int imu=0; imu<muonArr->GetEntries(); imu++) {
      const Muon *mu = (*muonArr)[imu];
      SelectionStatus denomsel;
      denomsel |= denominatorSelection(ctrl,mu,NULL,vtx,pfArr);
      if( !denomsel.passPre() ) continue;
      if(!passGhost(goodMuons,mu)) continue;
      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,vtx));
      } else {
        failingLeptons.push_back(fillSimpleLepton(mu,&status,vtx));
      }
    }
    for(Int_t i=0; i<electronArr->GetEntries(); i++) {
      const Electron *ele = (*electronArr)[i];
      bool muonMatch=false;
      for(unsigned imu=0; imu<muonArr->GetEntries(); imu++) {
        const Muon *mu = (*muonArr)[imu];
        if(dr(mu,ele) < 0.05) muonMatch = true;
      }
      if(muonMatch) continue;
      SelectionStatus denomsel;
      denomsel |= denominatorSelection(ctrl,NULL,ele,vtx,NULL);
      if( !denomsel.passPre() ) 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,vtx));
      } else {
        failingLeptons.push_back(fillSimpleLepton(ele,&status,vtx));
      }
    }

    // look for z pairs
    int NZCandidates=0;
    float bestMass=-999;
    const Muon *zmu1=0,*zmu2=0;
    const 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],NULL,NULL);
        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(NULL,NULL,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,vtx));
    data.Z1leptons.push_back(fillSimpleLepton(zmu2,zele2,status2,vtx));
    data.Z2leptons.push_back(fillSimpleLepton(zmu1,zele1,status1,vtx));
    data.Z2leptons.push_back(fillSimpleLepton(zmu2,zele2,status2,vtx));

    // 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(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) 
{
  if(ctrl.jsonFile != "")
    rlrm.AddJSONFile(string(ctrl.jsonFile));
}
//----------------------------------------------------------------------------------------
float getMass(const Muon *zmu1, const Muon *zmu2, const Electron *zele1, const 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 if(zmu1 && zele1) {
    lep1.SetPtEtaPhiM(zmu1->Pt(),zmu1->Eta(),zmu1->Phi(),MUON_MASS);
    lep2.SetPtEtaPhiM(zele1->Pt(),zele1->Eta(),zele1->Phi(),ELECTRON_MASS);
  } else assert(0);

  return (lep1 + lep2).M();
}
//----------------------------------------------------------------------------------------
SelectionStatus denominatorSelection(ControlFlags &ctrl, const Muon *mu, const Electron *ele, const Vertex *vtx, Array<PFCandidate> *pfArr)
{
  assert(mu || ele);
  assert(!(mu && ele));
  if(mu) {
    assert(pfArr);
    if(ctrl.muSele=="noD0Ip")
      return muonPreSelectionNoD0IP(ctrl,mu,vtx,pfArr);
    else if(ctrl.muSele=="d0Ip")
      return muonReferencePreSelection(ctrl,mu,vtx,pfArr);
    else assert(0);
  } else if(ele) {
    if(ctrl.eleSele=="noD0Ip")
      return electronPreSelectionNoD0IP(ctrl,ele,vtx);
    else if(ctrl.eleSele=="d0Ip")
      return electronReferencePreSelection(ctrl,ele,vtx);
    else assert(0);
  } else assert(0);
}  
//----------------------------------------------------------------------------------------
SelectionStatus analysisSelection(ControlFlags &ctrl, const Muon *mu, const Vertex *vtx, Array<PFCandidate> *pfArr,
                                  Array<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 Electron *ele, const Vertex *vtx, Array<PFCandidate> *pfArr,
                                  Array<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,
                               Array<Muon> *muonArr,
                               Array<Electron> *electronArr,
                               Array<PFJet> *jetArr,
                               Array<PFCandidate> *pfArr,
                               Array<PileupInfo>  *puArr,
                               Array<PileupEnergyDensity>  *puDArr,
                               Array<Vertex>               *vtxArr,
                               const Vertex         *vtx,
                               const Muon *zmu1,
                               const Muon *zmu2,
                               const Electron *zele1,
                               const Electron *zele2)
{
  vector<SimpleLepton> fakes;

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

    if(mu==zmu1 || mu==zmu2)     continue;
    if(!passResKill(zmu1,zmu2,zele1,zele2,mu,NULL)) continue;
      
    SelectionStatus denomsel;
    denomsel |= denominatorSelection(ctrl,mu,NULL,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 Electron* ele = (*electronArr)[iele];

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

    if(ele==zele1 || ele==zele2)     continue;
    if(!passResKill(zmu1,zmu2,zele1,zele2,NULL,ele)) continue;
      
    SelectionStatus denomsel;
    denomsel |= denominatorSelection(ctrl,NULL,ele,vtx,NULL);
    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, const Vertex *vtx)
{
  if(mu)       return fillSimpleLepton(mu, ctrl, vtx);
  else if(ele) return fillSimpleLepton(ele,ctrl, vtx);
  else assert(0);
  return SimpleLepton();
}
//----------------------------------------------------------------------------------------
bool passGhost(vector<const Muon*> goodMuons, const Muon *mu)
{
  for(unsigned imu=0; imu<goodMuons.size(); imu++) {
    const Muon *goodMu = goodMuons[imu];
    if(dr(goodMu,mu) < 0.02) return false;
  }
  return true;
}
//----------------------------------------------------------------------------------------
bool passResKill(const Muon *zmu1, const Muon *zmu2, const Electron *zele1, const Electron *zele2,
                 const Muon *mu, const Electron *ele)
// require mll > 4 between a lepton (either mu or ele) and the OS lepton from the Z
{
  float lepCharge = mu ? mu->Charge() : ele->Charge();
  float mass(0);
  if(zmu1) {
    if(zmu1->Charge()*lepCharge < 0)      mass = getMass(zmu1, mu ? mu : NULL, mu ? NULL : ele, NULL);
    else if(zmu2->Charge()*lepCharge < 0) mass = getMass(zmu2, mu ? mu : NULL, mu ? NULL : ele, NULL);
    else assert(0);
  } else if(zele1) {
    if(zele1->Charge()*lepCharge < 0)      mass = getMass(mu ? mu : NULL, NULL, zele1, mu ? NULL : ele);
    else if(zele2->Charge()*lepCharge < 0) mass = getMass(mu ? mu : NULL, NULL, zele2, mu ? NULL : ele);
    else assert(0);
  } else assert(0);

  return mass > 4;
}
Revision:	1.16
Committed:	Mon Dec 17 17:25:55 2012 UTC (12 years, 5 months ago) by dkralph
Content type:	text/plain
Branch:	MAIN
CVS Tags:	compiled
Changes since 1.15:	+11 -2 lines
Log Message:	* empty log message *