Selection/src/ReferenceSelection.cc

// system headers
#include <map>
#include <utility>

// mit headers
#include "Vertex.h"

// 4L stuff
#include "SelectionStatus.h" 
#include "EventData.h" 
#include "SimpleLepton.h" 
#include "EfficiencyWeightsInterface.h"
#include "ElectronSelection.h"
#include "MuonSelection.h"
#include "IsolationSelection.h"
#include "ReferenceSelection.h"
#include "Selection.h"
#include "CommonDefs.h"
#include "SelectionDefs.h"
#include "FSR.h"
#include "SelectionFuncs.h"
#include "ElectronMomentumCorrection.h"
#include "MuonMomentumCorrection.h"


extern vector<SimpleLepton> failingLeptons;
extern vector<SimpleLepton> passingLeptons;

extern map<unsigned,float>       evtrhoMap;
extern bool passes_HLT;
extern map<TString, map<TString,int>* > counts;
extern ElectronMomentumCorrection electron_momentum_correction;
extern MuCorr *muCorr;

//----------------------------------------------------------------------------------------
void updateSimpleLepton(SimpleLepton &tmplep); 
//----------------------------------------------------------------------------------------
void increment(TString cutName, int z1type=0, int z2type=0)
//----------------------------------------------------------------------------------------
// increment counters that keep track up event counts at each cut
{
  assert(counts.find(cutName) != counts.end());
  (*(counts[cutName]))["all"] ++;
  if(z1type!=0 && z2type!=0) {
    assert(counts.find(cutName) != counts.end());
    (*(counts[cutName]))[getChannel(z1type,z2type)] ++;
  }
}
//--------------------------------------------------------------------------------------------------
EventData apply_HZZ4L_reference_selection(ControlFlags &ctrl,           // input control 
                                           const EventHeader *info,     // input event info
                                           const Array<Vertex>       * vtxArr ,
                                           const Array<PFCandidate>  *pfCandidates,
                                           const Array<PileupEnergyDensity>  *puEnergyDensity,
                                           const Array<Electron> *electronArr,    // input electrons
                                           SelectionStatus (*ElectronPreSelector)( ControlFlags &, 
                                                                                   const Electron*,
                                                                                   const Vertex *),
                                           SelectionStatus (*ElectronIDSelector)( ControlFlags &, 
                                                                                  const Electron*,
                                                                                  const Vertex *), 
                                           SelectionStatus (*ElectronIsoSelector)( ControlFlags &, 
                                                                                   const Electron*, 
                                                                                   const Vertex *, 
                                                                                   const Array<PFCandidate> *, 
                                                                                   const Array<PileupEnergyDensity> *, 
                                                                                   ElectronTools::EElectronEffectiveAreaTarget,
                                                                                   vector<const PFCandidate*>),
                                           const Array<Muon> *muonArr,    // input muons
                                           SelectionStatus (*MuonPreSelector)( ControlFlags &, 
                                                                               const Muon*,
                                                                               const Vertex *,
                                                                               const Array<PFCandidate> *),
                                           SelectionStatus (*MuonIDSelector)( ControlFlags &, 
                                                                              const Muon*, 
                                                                              // const Vertex &),
                                                                              const Vertex *,
                                                                              const Array<PFCandidate> *),
                                           SelectionStatus (*MuonIsoSelector)( ControlFlags &, 
                                                                               const Muon*, 
                                                                               const Vertex *, 
                                                                               const Array<PFCandidate> *,
                                                                               const Array<PileupEnergyDensity> *, 
                                                                               MuonTools::EMuonEffectiveAreaTarget,
                                                                               vector<const PFCandidate*>)
                                          )
//--------------------------------------------------------------------------------------------------
{
  increment("start");
  EventData ret;
  failingLeptons.clear();
  passingLeptons.clear();

  MuonTools::EMuonEffectiveAreaTarget eraMu;
  ElectronTools::EElectronEffectiveAreaTarget eraEle;
  getEATargets(ctrl,eraMu,eraEle);

  if( ctrl.debug ) { 
    cout << "presel nlep: " << muonArr->GetEntries() + electronArr->GetEntries()
         << "\tnmuon: "    << muonArr->GetEntries()
         << "\tnelectron: " << electronArr->GetEntries()
         << endl;
  }

  // correct muon momentum
  if(ctrl.correct_muon_momentum) {
    assert(muCorr->corr2011 && muCorr->corr2012);
    for(int i=0; i<muonArr->GetEntries(); i++) {
      const Muon *const_mu = (Muon*)((*muonArr)[i]);
      Muon *mu = const_cast<Muon*>(const_mu);
      correct_muon_momentum(ctrl,muCorr,mu,info->RunNum());
    }
  }

  const Vertex * vtx;
  bool goodVertex = setPV( ctrl, vtxArr, vtx );
  if(!goodVertex) {
    if(ctrl.debug) cout << "found bad vertex" << endl;
    ret.status.setStatus(SelectionStatus::FAIL);
    return ret;
  }
  if(ctrl.debug)
    cout << "vtx :: ntrks: " << vtx->NTracksFit() << endl;

  //***********************************************************
  // Trigger Selection  
  //***********************************************************
  if( passes_HLT ) { 
    increment("trigger");
  } else {
    if(ctrl.debug) cout << "fails trigger" << endl;
    ret.status.setStatus(SelectionStatus::FAIL);
    return ret;
  }
  
  //***********************************************************
  // Lepton Selection 
  //***********************************************************
  vector<SimpleLepton> lepvec;
  vector<const PFCandidate*> photonsToVeto;

  if( ctrl.debug ) cout << "\tnMuons: " << muonArr->GetEntries() << endl;
  for(int i=0; i<muonArr->GetEntries(); i++) 
    {
      const Muon *mu = (Muon*)((*muonArr)[i]);      
      
      SelectionStatus musel;
      if(ctrl.debug) cout << "musel.status  before anything: " << musel.getStatus() << endl;

      musel |= (*MuonPreSelector)(ctrl,mu,vtx,pfCandidates);
      if(ctrl.debug) cout << "musel.status  after presel: " << musel.getStatus() << endl;
      if( !(musel.getStatus() & SelectionStatus::PRESELECTION) ) continue;
      if( ctrl.debug ) cout << endl;

      musel |= (*MuonIDSelector)(ctrl,mu,vtx,pfCandidates );
      if(ctrl.debug) cout << "musel.status  after ID: " << musel.getStatus() << endl;
      if( ctrl.debug ) cout << endl;

      // if we *do* FSR, isolation is calculated only after FSR recovery
      if( !(ctrl.doFSR) ) { 
        musel |=  (*MuonIsoSelector)(ctrl,mu,vtx,pfCandidates,puEnergyDensity,eraMu,photonsToVeto);
        if(ctrl.debug) cout << "musel.status  after iso: " << musel.getStatus() << endl;
      }

      if( ctrl.debug )
        cout << "muon:: pt: " << mu->Pt() << "\teta: " << mu->Eta() << "\tstatus: " << hex << musel.getStatus() << dec << endl;

      SimpleLepton tmplep;
      tmplep.vec.SetPtEtaPhiM(mu->Pt(), mu->Eta(), mu->Phi(), MUON_MASS);
      tmplep.type    = 13;
      tmplep.index   = i;
      tmplep.charge  = mu->Charge();
      tmplep.ip3dSig = mu->Ip3dPVSignificance();
      tmplep.is4l    = false;
      tmplep.isEB    = (fabs(mu->Eta()) < 1.479 ? 1 : 0 ); 
      tmplep.isLoose = musel.loose();
      tmplep.isTight = mu->IsPFMuon() || mu->IsGlobalMuon();
      tmplep.status  = musel;     
      tmplep.fsrRecoveryAttempted = false; // NOTE: this is *used* inside FSR.cc
      lepvec.push_back(tmplep);
      if( ctrl.debug ) cout << endl;
    }
    
    if( ctrl.debug ) { cout << "\tnElectron: " << electronArr->GetEntries() << endl; }
    for(int i=0; i<electronArr->GetEntries(); i++) 
      {
        const Electron *const_ele = (Electron*)((*electronArr)[i]);
        Electron *ele = const_cast<Electron*>(const_ele);

        // evaluate MVA *before* regression correction (but cut on mvaVal *after*)
        double mvaVal = getElectronIDMVAval(ctrl, ele, vtx);

        electron_momentum_correction.correct_electron_momentum(ctrl, ele,puEnergyDensity->At(0)->RhoKt6PFJets(),vtxArr->GetEntries());

        SelectionStatus elesel;
        if( ctrl.debug ) cout << "--> status before anything: " << hex << elesel.getStatus() << dec << endl;

        elesel |= (*ElectronPreSelector)(ctrl,ele,vtx);
        if( ctrl.debug ) cout << "--> status after presel: " << hex << elesel.getStatus() << dec << endl;
        if( !(elesel.getStatus() & SelectionStatus::PRESELECTION) ) continue;

        // elesel |= (*ElectronIDSelector)(ctrl,ele,vtx);
        elesel |= passElectronIDMVA(ctrl, mvaVal, ele);
        if( ctrl.debug ) cout << "--> status after ID: " << hex << elesel.getStatus() << dec << endl;

        if( !(ctrl.doFSR) ) { 
          elesel |= (*ElectronIsoSelector)(ctrl,ele,vtx,pfCandidates,puEnergyDensity,eraEle,photonsToVeto);
          if( ctrl.debug ) cout << "--> status after iso: " << hex << elesel.getStatus() << dec << endl;
        }
        
        if( ctrl.debug )
          cout << "\tscEt: " << ele->SCluster()->Et()
               << "\tscEta: " << ele->SCluster()->Eta()
               << "\tstatus: " << hex << elesel.getStatus() << dec
               << endl;

        SimpleLepton tmplep;
        tmplep.vec.SetPtEtaPhiM( ele->Pt(), ele->Eta(), ele->Phi(), ELECTRON_MASS );
        tmplep.type    = 11;
        tmplep.index   = i;
        tmplep.charge  = ele->Charge();
        tmplep.ip3dSig = ele->Ip3dPVSignificance();
        tmplep.is4l    = false;
        tmplep.isEB    = ele->IsEB();
        tmplep.scID    = ele->SCluster()->GetUniqueID();
        tmplep.isTight = elesel.tight();
        tmplep.isLoose = elesel.loose();
        tmplep.status  = elesel;
        tmplep.fsrRecoveryAttempted = false; // NOTE: this is *used* inside FSR.cc
        lepvec.push_back(tmplep);
        if( ctrl.debug ) cout << endl;
      }

    //********************************************************
    // Dump Stuff
    //********************************************************
    sort( lepvec.begin(), lepvec.end(), SimpleLepton::lep_pt_sort );
    int nmu=0, nele=0;
    for( int i=0; i<lepvec.size(); i++ ) {
      if(ctrl.debug) { 
        bitset<16> tmpbits(lepvec[i].status.getStatus());
        cout << "lepvec :: evt: " << setw(8) << info->EvtNum() << " index: " << setw(3) << i << " type: " << setw(4) << lepvec[i].type 
             << " pt: " << setw(11) << lepvec[i].vec.Pt() << " eta: " << setw(11) << lepvec[i].vec.Eta() << " status: " << tmpbits;
        if(!ctrl.doFSR)
          cout << "\tpf: " << lepvec[i].status.isoPF04 << "\tch: " << lepvec[i].status.chisoPF04 << "\tga: " << lepvec[i].status.gaisoPF04 
               << "\tne: " << lepvec[i].status.neisoPF04 << endl;
        if( abs(lepvec[i].type) == 11 ) {
          const Electron *tmpele = (Electron*)((*electronArr)[lepvec[i].index]);
          cout << "\tSCeta: " << tmpele->SCluster()->Eta()
               << "\tidMVA: " << lepvec[i].status.idMVA;
        }
        cout << endl;
      }
      if( abs(lepvec[i].type) == 11 ) nele++;
      else nmu++;
    }
    if( ctrl.debug ) { 
      cout << "postsel nlep: " << lepvec.size()
           << "\tnmuon: " << nmu
           << "\tnelectron: " << nele
           << endl;
    }
    
    //********************************************************
    // Step 2: Lepton Cross-Cleaning: remove electrons that are next to 'good'ish muons
    //********************************************************
    vector<vector<SimpleLepton>::iterator> electrons_to_erase; 
    for (vector<SimpleLepton>::iterator it1=lepvec.begin(); 
         it1 != lepvec.end(); it1++ ) {
      if ( abs(it1->type) != 11 ) continue;
      TVector3 evec = it1->vec.Vect();
      
      bool erase_this_electron=false;
      for (vector<SimpleLepton>::iterator it2=lepvec.begin(); 
           it2 != lepvec.end(); it2++ ) {
        if ( it2 == it1 ) continue;
        if ( abs(it2->type) != 13 ) continue;
        // from ANL:
        // ??-  if( !it2->status.passPre() || (!it2->isPFMuon && !it2->isGlobalMuon)  )       continue;
        if( !it2->status.passPre() ) continue;
        if( !it2->isTight )       continue; // NOTE: isTight is set to it2->isPFMuon || it2->isGlobalMuon

        TVector3 mvec = it2->vec.Vect();
        if ( evec.DrEtaPhi(mvec) < 0.05 ) {
          erase_this_electron=true;
          break;
        }
      }
      if( erase_this_electron ) {
        if( ctrl.debug ) cout << "erasing electron with pt  " <<  it1->vec.Pt() << endl;
        electrons_to_erase.push_back(it1);
      }
    }
    for( int i=0; i<electrons_to_erase.size(); i++ ) {
      lepvec.erase(electrons_to_erase[i]);
    }

    //********************************************************
    // Step 3: Good Leptons (remove non-id-and-pre'd leptons)
    //********************************************************
    vector<double> pt_of_leptons_to_erase; 
    for (int i=0; i<lepvec.size(); i++ ) {
      bool already_pushed=false;
      if( !(lepvec[i].status.looseIDAndPre()) ) { 
        pt_of_leptons_to_erase.push_back(lepvec[i].vec.Pt());
        already_pushed = true;
        if(ctrl.debug)
          cout << "pushing failed lepton type: " << lepvec[i].type
               << "\tpt: " << lepvec[i].vec.Pt()
               << "\teta: " << lepvec[i].vec.Eta()
               << endl;
        failingLeptons.push_back(lepvec[i]); // these should pass preselection
      } else {
        passingLeptons.push_back(lepvec[i]);
      }
      if( !already_pushed && fabs(lepvec[i].ip3dSig)>4 )  
        pt_of_leptons_to_erase.push_back(lepvec[i].vec.Pt()); 
    }
   for( int i=0; i<pt_of_leptons_to_erase.size(); i++ ) {
      for( vector<SimpleLepton>::iterator it=lepvec.begin();
           it != lepvec.end(); it++ ) {
        SimpleLepton flep = *it;
        if( flep.vec.Pt() != pt_of_leptons_to_erase[i] ) continue;
        if(ctrl.debug) cout << "erasing lepton : "
                            << flep.vec.Pt() << "\t"
                            << flep.type << "\t"
                            << endl;
        lepvec.erase(it);
        break;
      }
    }
    if( ctrl.debug ) cout << "good leptons : " << lepvec.size() << endl;

    //********************************************************
    //  Ghost Removal
    //********************************************************
    vector<vector<SimpleLepton>::iterator> leptons_to_erase; 
    for (vector<SimpleLepton>::iterator it1=lepvec.begin(); it1 != lepvec.end(); it1++ ) {
      TVector3 lep1 = it1->vec.Vect();
      
      bool erase_this_lepton=false;
      for (vector<SimpleLepton>::iterator it2=it1+1; it2 != lepvec.end(); it2++ ) {
        TVector3 lep2 = it2->vec.Vect();
        
        if ( lep1.DrEtaPhi(lep2) <= 0.02 ) {
          erase_this_lepton=true;
          break;
        }
      }
      if( erase_this_lepton ) {
        if( ctrl.debug ) cout << "erasing ghost with pt  " <<  it1->vec.Pt() << endl;
        leptons_to_erase.push_back(it1);
      }
    }
    for( int i=0; i<leptons_to_erase.size(); i++ ) {
      lepvec.erase(leptons_to_erase[i]);
    }

    //********************************************************
    // Step 4: Z candidate preselection
    //********************************************************
    vector<pair<int,int> > ZCandidates;
    vector<pair<SimpleLepton,SimpleLepton> > ZCandidatesLeptons;
    vector<pair<SelectionStatus,SelectionStatus> > ZCandidatesSelectionStatus;
    for(int i = 0; i < lepvec.size(); ++i) {
      for(int j = i+1; j < lepvec.size(); ++j) {
        if( abs(lepvec[i].type) != abs(lepvec[j].type) )      continue;
        if( lepvec[i].charge == lepvec[j].charge )            continue;

        TLorentzVector Zvec = (lepvec[i].vec+lepvec[j].vec);

        vector<SimpleLepton> lepvec_i = lepvec;
        vector<SimpleLepton> lepvec_j = lepvec;

        if( ctrl.doFSR ) {
          if(ctrl.debug) cout << endl << "----------------> FSR ("<<i<<","<<j<<") <----------------------" << endl;
          photonsToVeto.clear();
          float old_pt_i = lepvec[i].vec.Pt();
          float old_pt_j = lepvec[j].vec.Pt();
          float old_M  =  Zvec.M(); 
          
          if(ctrl.debug) cout << "i: " << i << endl;
          if( abs(lepvec[i].type) == 13 ) {
            const Muon *mu = (Muon*)((*muonArr)[lepvec[i].index]);      
            Muon * newmu = const_cast<Muon *>(mu);
            if( recover_typeI_Photon( ctrl, ret, newmu, i, lepvec_i, pfCandidates, electronArr, &Zvec, photonsToVeto ) )
              if(ctrl.debug) cout << "  FSR TYPEI :: oldpt: " << old_pt_i << "\tnewpt: " << lepvec_i[i].vec.Pt() << "\tindex: " << i << endl;
            if( recover_typeII_Photon( ctrl, ret, newmu, i, lepvec_i, pfCandidates ) )
              if(ctrl.debug) cout << "  FSR TYPEII :: oldpt: " << old_pt_i << "\tnewpt: " << lepvec_i[i].vec.Pt() << "\tindex: " << i << endl;
          } else { 
            const Electron *el = (Electron*)((*electronArr)[lepvec[i].index]);      
            Electron* newel = const_cast<Electron*>(el);
            if( recover_typeI_Photon( ctrl, ret, newel, i, lepvec_i, pfCandidates, electronArr, &Zvec, photonsToVeto ) )
              if(ctrl.debug) cout << "  FSR TYPEI :: oldpt: " << old_pt_i << "\tnewpt: " << lepvec_i[i].vec.Pt() << "\tindex: " << i << endl;
          }

          if(ctrl.debug) cout << "j: " << j << endl;
          ChargedParticle *cp;
          if( abs(lepvec[j].type) == 13 ) {
            const Muon *mu = (Muon*)((*muonArr)[lepvec[j].index]);      
            Muon * newmu = const_cast<Muon *>(mu);
            if( recover_typeI_Photon( ctrl, ret, newmu, j, lepvec_j, pfCandidates, electronArr, &Zvec, photonsToVeto ) )
              if(ctrl.debug) cout << "  FSR TYPEI :: oldpt: " << old_pt_j << "\tnewpt: " << lepvec_j[j].vec.Pt() << "\tindex: " << j << endl;
            if( recover_typeII_Photon( ctrl, ret, newmu, j, lepvec_j, pfCandidates ) )
              if(ctrl.debug) cout << "  FSR TYPEII :: oldpt: " << old_pt_j << "\tnewpt: " << lepvec_j[j].vec.Pt() << "\tindex: " << j << endl;
          } else { 
            const Electron *el = (Electron*)((*electronArr)[lepvec[j].index]);      
            Electron* newel = const_cast<Electron*>(el);
            if( recover_typeI_Photon( ctrl, ret, newel, j, lepvec_j, pfCandidates, electronArr, &Zvec, photonsToVeto ) )
              if(ctrl.debug) cout << "  FSR TYPEI :: oldpt: " << old_pt_j << "\tnewpt: " << lepvec_j[j].vec.Pt() << "\tindex: " << j << endl;
          }

          //
          // now fix isolation
          //
          if( abs(lepvec[i].type) == 11 ) {
            const Electron *el = (Electron*)((*electronArr)[lepvec_i[i].index]);      
            lepvec_i[i].status |= 
              (*ElectronIsoSelector)(ctrl,el,vtx,pfCandidates,puEnergyDensity,eraEle,photonsToVeto);      

          } else { 
            const Muon *mu = (Muon*)((*muonArr)[lepvec_i[i].index]);      
            lepvec_i[i].status |= 
              (*MuonIsoSelector)(ctrl,mu,vtx,pfCandidates,puEnergyDensity,eraMu,photonsToVeto);   
          }
          updateSimpleLepton(lepvec_i[i]); 

          if( abs(lepvec[j].type) == 11 ) {
            const Electron *el = (Electron*)((*electronArr)[lepvec_j[j].index]);      
            lepvec_j[j].status |= 
              (*ElectronIsoSelector)(ctrl,el,vtx,pfCandidates,puEnergyDensity,eraEle,photonsToVeto);      
          } else { 
            const Muon *mu = (Muon*)((*muonArr)[lepvec_j[j].index]);      
            lepvec_j[j].status |= 
              (*MuonIsoSelector)(ctrl,mu,vtx,pfCandidates,puEnergyDensity,eraMu,photonsToVeto);   
          }
          updateSimpleLepton(lepvec_j[j]); 

          float new_M    =   (lepvec_i[i].vec+lepvec_j[j].vec).M();
          float new_pt_i = lepvec_i[i].vec.Pt(); 
          float new_pt_j = lepvec_j[j].vec.Pt(); 
          if( ctrl.debug ) {
            cout << "  ----> post FSR Z:";
            cout << "      oldM: " << old_M << "\tnewM:" << new_M << endl;
            cout << "      old_pt_i: " << old_pt_i << "\tnew_pt_i:" << new_pt_i << endl;
            cout << "      old_pt_j: " << old_pt_j << "\tnew_pt_j:" << new_pt_j << endl;
          }
          
        } // doFSR
          
        if( !(lepvec_i[i].status.loose()) || !(lepvec_j[j].status.loose()) ) {
          if(ctrl.debug) {
            bitset<16> tmpbits_i(lepvec_i[i].status.getStatus()),tmpbits_j(lepvec_j[j].status.getStatus());
            cout << "  leptons fail selection (i: " << tmpbits_i << ", j: " << tmpbits_j << ")" << endl;
          }
          continue;
        }
        ZCandidates.push_back(pair<int,int> (i,j) );
        ZCandidatesLeptons.push_back(pair<SimpleLepton,SimpleLepton> (lepvec_i[i],lepvec_j[j]) );
        if( ctrl.debug ) cout << "Z candidate ("<< i << "," << j << ")"
                              << "\tmass: " << (lepvec_i[i].vec+lepvec_j[j].vec).M() << endl;
      }
    }
    if( ZCandidates.size() > 0 ) { 
      ret.status.selectionBits.flip(PASS_ZCANDIDATE);
      if( ctrl.debug ) cout << "event has >0  Z candidates (" << ZCandidates.size() << " of em)" << endl;
    } else {
      ret.status.setStatus(SelectionStatus::FAIL);
      return ret;
    }

    //
    // !!!!!!!!!!!!!! Z1 SELECTED HERE
    //
    int best_Z1_index;
    float best_Z1_mass = 9999.;
    TLorentzVector Z1vec;
    for( int i=0; i<ZCandidates.size(); i++ ) {
      //       TLorentzVector tmpZ1vec = (lepvec[ZCandidates[i].first].vec) + 
      //        (lepvec[ZCandidates[i].second].vec); 
      TLorentzVector tmpZ1vec = (ZCandidatesLeptons[i].first.vec) + 
        (ZCandidatesLeptons[i].second.vec); 
      if( fabs(tmpZ1vec.M()-Z_MASS) < fabs(best_Z1_mass-Z_MASS) ) {
        best_Z1_index=i;
        best_Z1_mass=tmpZ1vec.M();
        Z1vec = tmpZ1vec;
      }
    }
    // update the leptons, damn FSR ...
    //     ret.Z1leptons.push_back(lepvec[ZCandidates[best_Z1_index].first]);
    //     ret.Z1leptons.push_back(lepvec[ZCandidates[best_Z1_index].second]);
    ret.Z1leptons.push_back(ZCandidatesLeptons[best_Z1_index].first);
    ret.Z1leptons.push_back(ZCandidatesLeptons[best_Z1_index].second);


    cout << "best mZ1:  " << best_Z1_mass << " from ("  << 
      ZCandidates[best_Z1_index].first << "," << ZCandidates[best_Z1_index].second << ")" << endl;

    //******************************************************************************
    // Step 6.3 : require Z1 with  40<m<120
    //******************************************************************************
    if( Z1vec.M() > 40. && Z1vec.M() < 120. ) { 
      ret.status.selectionBits.flip(PASS_GOODZ1);
    } else {
      ret.status.setStatus(SelectionStatus::FAIL);
      return ret;
    }

    //******************************************************************************
    // Step 6.3 : 4 good leptons 
    //******************************************************************************
    if( lepvec.size() >= 4 ) { 
      if( ctrl.debug) cout << "pass4L: "  << info->EvtNum() << endl;
      ret.status.selectionBits.flip(PASS_4L);
    } else {
      ret.status.setStatus(SelectionStatus::FAIL);
      return ret;
    }
    int nEl=0, nMu=0;
    for( int i=0; i<4; i++ ) {
      if(abs(lepvec[i].type) == 11 ) nEl++;
      if(abs(lepvec[i].type) == 13 ) nMu++;
    }

    //********************************************************
    // Step 5: ZZ candidates
    //********************************************************
    int nZZCandidates=0;
    vector<pair<int,int> > ZZCandidates;
    int Z2type;
    if(ctrl.debug) cout << "looping over " << ZCandidates.size() << " Z candidates: " << endl;
    for(int z2index=0; z2index<ZCandidates.size(); ++z2index) {
      int z1index = best_Z1_index;
      if ( z2index == z1index )                                         { if(ctrl.debug) cout << "  Z2 candidate failing z2i = z1i" << endl; continue; }
      if( ZCandidates[z1index].first  == ZCandidates[z2index].first )   { if(ctrl.debug) cout << "  Z2 candidate failing icheck1" << endl; continue; }
      if( ZCandidates[z1index].first  == ZCandidates[z2index].second )  { if(ctrl.debug) cout << "  Z2 candidate failing icheck2" << endl; continue; }
      if( ZCandidates[z1index].second == ZCandidates[z2index].first )   { if(ctrl.debug) cout << "  Z2 candidate failing icheck3" << endl; continue; }
      if( ZCandidates[z1index].second == ZCandidates[z2index].second )  { if(ctrl.debug) cout << "  Z2 candidate failing icheck4" << endl; continue; }

      ZZCandidates.push_back(pair<int,int> (z1index,z2index));
      Z2type = abs(ZCandidatesLeptons[z2index].first.type);
    }
    if( ZZCandidates.size() > 0 ) { 
      if( ctrl.debug) cout << "passZZ: "  << info->EvtNum() << endl;
      ret.status.selectionBits.flip(PASS_ZZCANDIDATE);
      if( ZZCandidates.size() > 1 ) Z2type=999;
      if( ctrl.debug ) cout << "nZZcandidates: " << ZZCandidates.size() << endl;
      if( ctrl.debug ) { 
        cout << "evt:  " << info->EvtNum() << "\tnZZcandidates: " << ZZCandidates.size() 
             << "\tZ1f: " << abs(ret.Z1leptons[0].type) << "\tZ2f: " << Z2type << endl;
        cout << "------------------------------------------------------- ( un-fsrd )" << endl;
        for( int l=0; l<lepvec.size(); l++ ) lepvec[l].print();
        cout << "-------------------------------------------------------" << endl;
      }
    } else {
      ret.status.setStatus(SelectionStatus::FAIL);
      return ret;
    }

    //
    // !!!!!!!!!!!!!! Z2 SELECTED HERE
    //
    int best_Z2_index;
    float best_Z2_pt = -1.;
    for( int i=0; i<ZZCandidates.size(); i++ ) {
      int z2index = ZZCandidates[i].second;
      TLorentzVector Z2 = (ZCandidatesLeptons[z2index].first.vec) + 
        (ZCandidatesLeptons[z2index].second.vec); 
      if( Z2.Pt() > best_Z2_pt ) {
        best_Z2_index=z2index;
        best_Z2_pt=Z2.Pt();
      }
    }
    // update the leptons, damn FSR ...
    //    ret.Z2leptons.push_back(lepvec[ZCandidates[best_Z2_index].first]);
    //    ret.Z2leptons.push_back(lepvec[ZCandidates[best_Z2_index].second]);
    ret.Z2leptons.push_back(ZCandidatesLeptons[best_Z2_index].first);
    ret.Z2leptons.push_back(ZCandidatesLeptons[best_Z2_index].second);

    int theZ1type = abs(ZCandidatesLeptons[best_Z1_index].first.type);
    int theZ2type = abs(ZCandidatesLeptons[best_Z2_index].first.type);

    if(ctrl.debug) cout << "best mZ2:  " << (ret.Z2leptons[0].vec+ret.Z2leptons[1].vec).M() << " from ("
                        << ZCandidates[best_Z2_index].first << "," << ZCandidates[best_Z2_index].second << ")" << endl;

    increment("sfOsHiPt",theZ1type,theZ2type);

    if(ctrl.debug) cout << "ZZ :: evt: " << info->EvtNum()
                        << "\tmZ1: " << (ret.Z1leptons[0].vec+ret.Z1leptons[1].vec).M()  
                        << "\tmZ2: " << (ret.Z2leptons[0].vec+ret.Z2leptons[1].vec).M()  
                        << endl;


    //******************************************************************************
    // Step 6.4 : require Z2 with  4<m<120
    //******************************************************************************
    TLorentzVector Z2vec = (ZCandidatesLeptons[best_Z2_index].first.vec) + 
      (ZCandidatesLeptons[best_Z2_index].second.vec); 
    if( Z2vec.M() > 4 && Z2vec.M() < 120 ) { 
      ret.status.selectionBits.flip(PASS_GOODZ2);
      increment("4<mZ2<120",theZ1type,theZ2type);
    } else {
      ret.status.setStatus(SelectionStatus::FAIL);
      return ret;
    }


    //******************************************************************************
    // Step 6.1 : any two leptons 20/10
    //******************************************************************************
    vector<SimpleLepton> zzleptons;
    zzleptons.push_back( ZCandidatesLeptons[best_Z1_index].first );
    zzleptons.push_back( ZCandidatesLeptons[best_Z1_index].second );
    zzleptons.push_back( ZCandidatesLeptons[best_Z2_index].first );
    zzleptons.push_back( ZCandidatesLeptons[best_Z2_index].second ); 
    int nlep_above_10=0,nlep_above_20=0;
    for( int i=0; i<zzleptons.size(); i++ ) {
      if( zzleptons[i].vec.Pt() > 10 ) nlep_above_10++;
      if( zzleptons[i].vec.Pt() > 20 ) nlep_above_20++;
    }
    if( nlep_above_10 > 1 && nlep_above_20 > 0 ) { 
      ret.status.selectionBits.flip(PASS_ZZ_20_10);
      increment("pt>20,10",theZ1type,theZ2type);
      if( ctrl.debug ) cout << "passess 20/10 ..." << endl;
    } else {
      ret.status.setStatus(SelectionStatus::FAIL);
      return ret;
    }
    
    
    //******************************************************************************
    // Step 6.5 : resonance killing (4/4)
    //******************************************************************************
    bool resonance = false;
    for( int i=0; i<zzleptons.size(); i++ ) {
      for( int j=i+1; j<zzleptons.size(); j++ ) {
        if( zzleptons[i].charge == zzleptons[j].charge ) continue; // 4/4
        if( (zzleptons[i].vec+zzleptons[j].vec).M() <= 4. ) {
          resonance = true;
          break;
        }
      }
    }
    if( !resonance ) { 
      ret.status.selectionBits.flip(PASS_RESONANCE);
      increment("mll>4",theZ1type,theZ2type);
      if( ctrl.debug ) cout << "\tpasses resonance killing ... " << endl;
    } else {
      ret.status.setStatus(SelectionStatus::FAIL);
      return ret;
    }
    
    
    //******************************************************************************
    // Step 6.6 : m(4l) > 70 , m(4l) > 100
    //******************************************************************************
    TLorentzVector zzvec = (ZCandidatesLeptons[best_Z1_index].first.vec) +
      (ZCandidatesLeptons[best_Z1_index].second.vec) +
      (ZCandidatesLeptons[best_Z2_index].first.vec) +
      (ZCandidatesLeptons[best_Z2_index].second.vec);
    
    if(ctrl.debug)
      cout << "forming zz from: "
           << setw(9) << (ZCandidatesLeptons[best_Z1_index].first.vec).Pt()
           << setw(9) << (ZCandidatesLeptons[best_Z1_index].second.vec).Pt()
           << setw(9) << (ZCandidatesLeptons[best_Z2_index].first.vec).Pt()
           << setw(9) << (ZCandidatesLeptons[best_Z2_index].second.vec).Pt() << endl;

    if( zzvec.M() > 70. ) {
      if(ctrl.debug) cout << "passes mzz > 70, mzz: " << zzvec.M() << endl;
      ret.status.selectionBits.flip(PASS_m4l_GT_70);
      increment("m4l>70",theZ1type,theZ2type);
    } else {
      ret.status.setStatus(SelectionStatus::FAIL);
      return ret;
    }

    if( (ret.Z2leptons[0].vec+ret.Z2leptons[1].vec).M() > 12 ){
      if(ctrl.debug) cout << "passes mZ2 > 12" << endl;
      increment("mZ2>12",theZ1type,theZ2type);
    } else {
      ret.status.setStatus(SelectionStatus::FAIL);
      return ret;
    }

    if( zzvec.M() > 100 ) { 
      if(ctrl.debug) cout << "passes mzz > 100, mzz: " << zzvec.M() << endl;
      increment("m4l>100",theZ1type,theZ2type);
    } else {
      cout << "NOTE: failed m4l>100 (" << zzvec.M() << "), but saving to ntuple anyway" << endl;
    }

    //***************************************************************
    // finish
    //***************************************************************
    
    TLorentzVector theZ1 = (ZCandidatesLeptons[best_Z1_index].first.vec) + 
        (ZCandidatesLeptons[best_Z1_index].second.vec); 
    TLorentzVector theZ2 = (ZCandidatesLeptons[best_Z2_index].first.vec) + 
        (ZCandidatesLeptons[best_Z2_index].second.vec); 
    TLorentzVector theZZ = theZ1 + theZ2;
    
    if(ret.fsrPhotons.size() > 0) {
      if(ctrl.debug)
        cout << "  fsr photon! pt: "
             << setw(12) << ret.fsrPhotons[0].vec.Pt()
             << " eta: " << setw(12) << ret.fsrPhotons[0].vec.Eta()
             << " phi: " << setw(12) << ret.fsrPhotons[0].vec.Phi() << endl;
    }
    if(ret.fsrPhotons.size() > 1) {
      if(ctrl.debug)
        cout << "  fsr photon! pt: "
             << setw(12) << ret.fsrPhotons[1].vec.Pt()
             << " eta: " << setw(12) << ret.fsrPhotons[1].vec.Eta()
             << " phi: " << setw(12) << ret.fsrPhotons[1].vec.Phi() << endl;
    }
    if( ctrl.debug ) cout  << "run: " << info->RunNum()   
                           << "\tevt: " << info->EvtNum()
                           << "\tZ1channel: " << theZ1type
                           << "\tZ2channel: " << theZ2type
                           << "\tmZ1: " << theZ1.M()
                           << "\tmZ2: " << theZ2.M()
                           << "\tm4l: " << theZZ.M()
                           << endl;
    
    ret.status.setStatus(SelectionStatus::EVTPASS);
    
    return ret;
}
//----------------------------------------------------------------------------
void updateSimpleLepton(SimpleLepton &tmplep) 
//----------------------------------------------------------------------------
{ 
  tmplep.isTight = tmplep.status.tight();
  tmplep.isLoose = tmplep.status.loose();
}
Revision:	1.21
Committed:	Tue Oct 23 11:39:21 2012 UTC (12 years, 6 months ago) by dkralph
Content type:	text/plain
Branch:	MAIN
Changes since 1.20:	+242 -507 lines
Log Message:	* empty log message *