Selection/src/ReferenceSelection.cc

#include "ReferenceSelection.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;
extern TrigInfo ti;


//we need to store the fsr photons because they are used to veto jets in the vbf selection
//the leptons in the lepton pair include the fsr photons
struct ZCandidate
{
  pair<SimpleLepton,SimpleLepton> lepton_pair;
  unsigned int which_lepton_has_fsr; //0 means neither lepton has an fsr photon, 1 means the first lepton in the pair has the fsr photon, 2 means second lepton in the pair has the fsr photon
  PFCandidate fsr_photon;
};

//----------------------------------------------------------------------------------------
void updateSimpleLepton(SimpleLepton &tmplep); 
//--------------------------------------------------------------------------------------------------
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*>)
                                          )
//--------------------------------------------------------------------------------------------------
{
  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);
    if(ctrl.debug) cout << "mu corr: " << endl;
    for(int i=0; i<muonArr->GetEntries(); i++) {
      const Muon *const_mu = (Muon*)((*muonArr)[i]);
      Muon *mu = const_cast<Muon*>(const_mu);
      double ptBefore = mu->Pt();
      correct_muon_momentum(ctrl,muCorr,mu,info->RunNum());
      if(ctrl.debug) cout << "    " << setw(12) << ptBefore << " --> " << setw(12) << mu->Pt() << " (" << fabs(ptBefore-mu->Pt())/ptBefore << ")" << endl;
    }
  }

  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(counts,"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;
      musel |= (*MuonPreSelector)(ctrl,mu,vtx,pfCandidates);

      if(ctrl.fakeScheme == "none") {  // regular selection
        if( !musel.passPre() ) continue;
      } else { // fake denominator selection
        SelectionStatus denomSel;
        denomSel |= muonReferencePreSelection(ctrl,mu,vtx,pfCandidates);
        if(ctrl.muSele != "none") { // don't apply *any* preselection, to get a super hi stat fake sample for BDT training
          if( !denomSel.passPre() ) continue;
        }
      }

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

      musel |= (*MuonIDSelector)(ctrl,mu,vtx,pfCandidates );

      // NOTE: if we do FSR this is *changed* later on
      musel |=  (*MuonIsoSelector)(ctrl,mu,vtx,pfCandidates,puEnergyDensity,eraMu,photonsToVeto);
      if(ctrl.debug) cout << "mu status after iso (before fsr) " << musel.getStatus() << 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);

        if(ctrl.debug) cout << setprecision(8) << "el:: scEt " << setw(12) << ele->SCluster()->Et() << " P: "
                            << setw(12) << ele->P() << " pT: " << setw(12) << ele->Pt() << " scEta " << setw(12) << ele->SCluster()->Eta() << setprecision(5) << endl;

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

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

        if(ctrl.debug) cout << setprecision(8) << "  corr el: scEt " << setw(12) << ele->SCluster()->Et() << " P: "
                            << setw(12) << ele->P() << " pT: " << setw(12) << ele->Pt() << " scEta " << setw(12) << ele->SCluster()->Eta() << setprecision(5) << endl;

        SelectionStatus elesel;
        elesel |= (*ElectronPreSelector)(ctrl,ele,vtx);

        if(ctrl.fakeScheme == "none") { // regular selection
          if( !elesel.passPre() ) continue;
        } else { // fake denominator selection
          SelectionStatus denomSel;
          denomSel |= electronReferencePreSelection(ctrl,ele,vtx);
          if(ctrl.eleSele != "none") { // don't apply *any* preselection, to get a super hi stat fake sample for BDT training
            if( !denomSel.passPre() ) continue;
          }
        }


        elesel |= passElectronIDMVA(ctrl, mvaVal, ele);

        // NOTE: if we do FSR this is *changed* later on
        elesel |= (*ElectronIsoSelector)(ctrl,ele,vtx,pfCandidates,puEnergyDensity,eraEle,photonsToVeto);
        if( ctrl.debug ) cout << "el status after iso (no fsr) " << 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.perr = combination_perr;
        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;
        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]);
    }

    // fill passing and failing leptons for fakes
    // (NOTE: isolation hasn't been corrected post-fsr, so this isn't strictly correct)
    // also note: sip cut not applied
    for (int i=0; i<lepvec.size(); i++ ) {
      if( lepvec[i].isLoose )
        passingLeptons.push_back(lepvec[i]);
      else
        failingLeptons.push_back(lepvec[i]);
    }

    //********************************************************
    // Step 3: Good Leptons (remove non-id-and-pre'd leptons, and apply |ip3ds| < 4)
    //********************************************************
    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;                          // indices in lepvec of the two leptons
    vector<ZCandidate > ZCandidatesLeptons; // fsr-corrected leptons that form the Z
    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);

          // copies of lepvec that have 1) the photons added to the momenta and 2) the fsrRecoveryAttempted flag set
        vector<SimpleLepton> lepvec_i = lepvec;
        vector<SimpleLepton> lepvec_j = lepvec;

        ZCandidate cand;
        cand.which_lepton_has_fsr = 0; //default value for when neither photon has an fsr photon

        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(); 

          const ChargedParticle *lepton(NULL);
          if(ctrl.debug) cout << "i: " << i << endl;
          if(abs(lepvec[i].type) == 13) lepton = dynamic_cast<const ChargedParticle *>((*muonArr)[lepvec[i].index]);
          else                          lepton = dynamic_cast<const ChargedParticle *>((*electronArr)[lepvec[i].index]);
          pair<TLorentzVector,int> photon_i = findFsrPhoton(ctrl, ret, lepton, i, lepvec_i, pfCandidates, electronArr, &Zvec);
          if(photon_i.second>=0 && ctrl.debug) cout << "  FSR :: oldpt: " << old_pt_i << "\tnewpt: " << lepvec_i[i].vec.Pt() << "\tindex: " << i << endl;

          lepton = NULL;
          if(ctrl.debug) cout << "j: " << j << endl;
          if( abs(lepvec[j].type) == 13 ) lepton = dynamic_cast<const ChargedParticle *>((*muonArr)[lepvec[j].index]);
          else                            lepton = dynamic_cast<const ChargedParticle *>((*electronArr)[lepvec[j].index]);
          pair<TLorentzVector,int> photon_j = findFsrPhoton(ctrl, ret, lepton, j, lepvec_j, pfCandidates, electronArr, &Zvec);
          if(photon_j.second>=0 && ctrl.debug) cout << "  FSR :: oldpt: " << old_pt_j << "\tnewpt: " << lepvec_j[j].vec.Pt() << "\tindex: " << j << endl;

          bool useI(false),useJ(false);
          if(photon_i.second >= 0) { // photon for i
            if(photon_j.second < 0) {
              useI = true;
            } else { // both have photons
              TLorentzVector lep_i(lepvec_i[i].vec),lep_j(lepvec_j[j].vec);
              float oldM   = (lep_i + lep_j).M();
              float newM_i = (lep_i + photon_i.first + lep_j).M();
              float newM_j = (lep_i + lep_j + photon_j.first).M();
              if(ctrl.debug) cout << "  two FSRs "  << setw(12) << oldM << setw(12) << newM_i << setw(12) << newM_j << endl;
              if(fabs(newM_i - Z_MASS) < fabs(newM_j - Z_MASS)) { // if mass with photon i is closer to z pole than with photon j
                useI = true;
              } else {
                useJ = true;
              }
            }
          } else if(photon_j.second >= 0) { // photon for j
            useJ = true;
          }

          assert(!(useI && useJ));
          if(useI) {
            cand.fsr_photon =  *((const PFCandidate*) pfCandidates->At(photon_i.second));
            cand.which_lepton_has_fsr = 1;
            lepvec_i[i].vec += photon_i.first;
            lepvec_i[i].fsrRecoveryAttempted = true;
            photonsToVeto.push_back((const PFCandidate*)(pfCandidates->At(photon_i.second)));
          } else if(useJ) {
            cand.fsr_photon = *((const PFCandidate*) pfCandidates->At(photon_j.second));
            cand.which_lepton_has_fsr = 2;
            lepvec_j[j].vec += photon_j.first;
            lepvec_j[j].fsrRecoveryAttempted = true;
            photonsToVeto.push_back((const PFCandidate*)(pfCandidates->At(photon_j.second)));
          }       
          
          // recompute isolation excluding the fsr photons
          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 && ret.fsrPhotons.size() > 0) {
            cout << "  ----> post FSR Z:";
            cout << "      oldM: "     << setprecision(8) << old_M << "\tnewM:"        << setprecision(8) << new_M << endl;
            cout << "      old_pt_i: " << setprecision(8) << old_pt_i << "\tnew_pt_i:" << setprecision(8) << new_pt_i << endl;
            cout << "      old_pt_j: " << setprecision(8) << old_pt_j << "\tnew_pt_j:" << setprecision(8) << new_pt_j << endl;
          }
          
        } // doFSR

        // apply isolation criteria with recomputed isolation
        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) );
        cand.lepton_pair = pair<SimpleLepton,SimpleLepton> (lepvec_i[i],lepvec_j[j]);
        
        ZCandidatesLeptons.push_back(cand );
        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;
    }

    //
    // Select Z1
    //
    int best_Z1_index;
    float best_Z1_mass = 9999.;
    TLorentzVector Z1vec;
    for( int i=0; i<ZCandidates.size(); i++ ) {
      TLorentzVector tmpZ1vec = (ZCandidatesLeptons[i].lepton_pair.first.vec) + (ZCandidatesLeptons[i].lepton_pair.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;
      }
    }
    ret.Z1leptons.push_back(ZCandidatesLeptons[best_Z1_index].lepton_pair.first);
    ret.Z1leptons.push_back(ZCandidatesLeptons[best_Z1_index].lepton_pair.second);
    if(ZCandidatesLeptons[best_Z1_index].which_lepton_has_fsr != 0){
      SimpleLepton photon;
      photon.vec.SetPtEtaPhiM( ZCandidatesLeptons[best_Z1_index].fsr_photon.Pt(), ZCandidatesLeptons[best_Z1_index].fsr_photon.Eta(), ZCandidatesLeptons[best_Z1_index].fsr_photon.Phi(), 0);
      ret.fsrPhotons.push_back(photon);
    }
    
    if(ctrl.debug)
      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 << "four leps" << endl;
      ret.status.selectionBits.flip(PASS_4L);
    } else {
      ret.status.setStatus(SelectionStatus::FAIL);
      return ret;
    }

    //********************************************************
    // Step 5: ZZ candidates
    //********************************************************
    int nZZCandidates=0;
    vector<pair<int,int> > ZZCandidates;
    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 fails z2i = z1i" << endl; continue; }
      if( ZCandidates[z1index].first  == ZCandidates[z2index].first )   { if(ctrl.debug) cout << "  Z2 fails icheck1" << endl; continue; }
      if( ZCandidates[z1index].first  == ZCandidates[z2index].second )  { if(ctrl.debug) cout << "  Z2 fails icheck2" << endl; continue; }
      if( ZCandidates[z1index].second == ZCandidates[z2index].first )   { if(ctrl.debug) cout << "  Z2 fails icheck3" << endl; continue; }
      if( ZCandidates[z1index].second == ZCandidates[z2index].second )  { if(ctrl.debug) cout << "  Z2 fails icheck4" << endl; continue; }

      ZZCandidates.push_back(pair<int,int> (z1index,z2index));
    }
    if(ZZCandidates.size() > 0) { 
      ret.status.selectionBits.flip(PASS_ZZCANDIDATE);
      if(ctrl.debug) { 
        cout << ZZCandidates.size() << " zz cands" << endl;
        cout << "-------------------------------------------------------" << endl;
        // for( int l=0; l<lepvec.size(); l++ ) lepvec[l].print();
        for(unsigned ican=0; ican<ZCandidatesLeptons.size(); ican++) {
          ZCandidatesLeptons[ican].lepton_pair.first.print();
          ZCandidatesLeptons[ican].lepton_pair.second.print();
        }
        cout << "-------------------------------------------------------" << endl;
        if(ctrl.debug)
          for(unsigned iph=0; iph<ret.fsrPhotons.size(); iph++)
            cout << "  fsr photon "
                 << setw(12) << ret.fsrPhotons[iph].vec.Pt()
                 << " eta: " << setw(12) << ret.fsrPhotons[iph].vec.Eta()
                 << " phi: " << setw(12) << ret.fsrPhotons[iph].vec.Phi() << endl;
      }
    } else {
      if(ctrl.debug) cout << "no zz cands" << endl;
      ret.status.setStatus(SelectionStatus::FAIL);
      return ret;
    }

    //
    // Select z2
    //
    int best_Z2_index;
    float best_Z2_pt_sum = -1.;
    for( int i=0; i<ZZCandidates.size(); i++ ) {
      int z2index = ZZCandidates[i].second;
      double pt_sum = ZCandidatesLeptons[z2index].lepton_pair.first.vec.Pt() + 
        ZCandidatesLeptons[z2index].lepton_pair.second.vec.Pt(); 
      if( pt_sum > best_Z2_pt_sum ) {
        best_Z2_index=z2index;
        best_Z2_pt_sum=pt_sum;
      }
    }
    ret.Z2leptons.push_back(ZCandidatesLeptons[best_Z2_index].lepton_pair.first); // push back the *fsr-corrected* leptons
    ret.Z2leptons.push_back(ZCandidatesLeptons[best_Z2_index].lepton_pair.second);
    if(ZCandidatesLeptons[best_Z2_index].which_lepton_has_fsr != 0){
      SimpleLepton photon;
      photon.vec.SetPtEtaPhiM( ZCandidatesLeptons[best_Z2_index].fsr_photon.Pt(), ZCandidatesLeptons[best_Z2_index].fsr_photon.Eta(), ZCandidatesLeptons[best_Z2_index].fsr_photon.Phi(), 0);
      ret.fsrPhotons.push_back(photon);
    }


    int theZ1type = abs(ZCandidatesLeptons[best_Z1_index].lepton_pair.first.type);
    int theZ2type = abs(ZCandidatesLeptons[best_Z2_index].lepton_pair.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(counts,"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].lepton_pair.first.vec) + 
      (ZCandidatesLeptons[best_Z2_index].lepton_pair.second.vec); 
    if( Z2vec.M() > 4 && Z2vec.M() < 120 ) { 
      ret.status.selectionBits.flip(PASS_GOODZ2);
      increment(counts,"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].lepton_pair.first );
    zzleptons.push_back( ZCandidatesLeptons[best_Z1_index].lepton_pair.second );
    zzleptons.push_back( ZCandidatesLeptons[best_Z2_index].lepton_pair.first );
    zzleptons.push_back( ZCandidatesLeptons[best_Z2_index].lepton_pair.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(counts,"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)
    //******************************************************************************

    //the resonance killing is done based on the leptons without fsr recovered to them
    vector<SimpleLepton > zzleptons_nofsr;

    //add all of the leptons to a vector and remove the fsr photon if there is one
    if(ZCandidatesLeptons[best_Z1_index].which_lepton_has_fsr == 1){
      SimpleLepton no_fsr_lepton = ZCandidatesLeptons[best_Z1_index].lepton_pair.first;
      TLorentzVector photon_vec;
      photon_vec.SetPtEtaPhiM(ZCandidatesLeptons[best_Z1_index].fsr_photon.Pt(),ZCandidatesLeptons[best_Z1_index].fsr_photon.Eta(),ZCandidatesLeptons[best_Z1_index].fsr_photon.Phi(),0);
      no_fsr_lepton.vec = no_fsr_lepton.vec - photon_vec;
      zzleptons_nofsr.push_back( no_fsr_lepton );
      ret.Z1leptons_without_fsr.push_back(no_fsr_lepton);
    }
    else {
      zzleptons_nofsr.push_back( ZCandidatesLeptons[best_Z1_index].lepton_pair.first);
      ret.Z1leptons_without_fsr.push_back(ZCandidatesLeptons[best_Z1_index].lepton_pair.first);
    }
    
    if(ZCandidatesLeptons[best_Z1_index].which_lepton_has_fsr == 2){
      SimpleLepton no_fsr_lepton = ZCandidatesLeptons[best_Z1_index].lepton_pair.second;
      TLorentzVector photon_vec;
      photon_vec.SetPtEtaPhiM(ZCandidatesLeptons[best_Z1_index].fsr_photon.Pt(),ZCandidatesLeptons[best_Z1_index].fsr_photon.Eta(),ZCandidatesLeptons[best_Z1_index].fsr_photon.Phi(),0);
      no_fsr_lepton.vec = no_fsr_lepton.vec - photon_vec;
      zzleptons_nofsr.push_back( no_fsr_lepton );
      ret.Z1leptons_without_fsr.push_back(no_fsr_lepton);
    }
    else{
       zzleptons_nofsr.push_back( ZCandidatesLeptons[best_Z1_index].lepton_pair.second);
       ret.Z1leptons_without_fsr.push_back(ZCandidatesLeptons[best_Z1_index].lepton_pair.second);
    }

    if(ZCandidatesLeptons[best_Z2_index].which_lepton_has_fsr == 1){
      SimpleLepton no_fsr_lepton = ZCandidatesLeptons[best_Z2_index].lepton_pair.first;
      TLorentzVector photon_vec;
      photon_vec.SetPtEtaPhiM(ZCandidatesLeptons[best_Z2_index].fsr_photon.Pt(),ZCandidatesLeptons[best_Z2_index].fsr_photon.Eta(),ZCandidatesLeptons[best_Z2_index].fsr_photon.Phi(),0);
      no_fsr_lepton.vec = no_fsr_lepton.vec - photon_vec;
      zzleptons_nofsr.push_back( no_fsr_lepton );
      ret.Z2leptons_without_fsr.push_back(no_fsr_lepton);
    }
    else {
      zzleptons_nofsr.push_back( ZCandidatesLeptons[best_Z2_index].lepton_pair.first);
      ret.Z2leptons_without_fsr.push_back(ZCandidatesLeptons[best_Z2_index].lepton_pair.first);
    }

    if(ZCandidatesLeptons[best_Z2_index].which_lepton_has_fsr == 2){
      SimpleLepton no_fsr_lepton = ZCandidatesLeptons[best_Z2_index].lepton_pair.second;
      TLorentzVector photon_vec;
      photon_vec.SetPtEtaPhiM(ZCandidatesLeptons[best_Z2_index].fsr_photon.Pt(),ZCandidatesLeptons[best_Z2_index].fsr_photon.Eta(),ZCandidatesLeptons[best_Z2_index].fsr_photon.Phi(),0);
      no_fsr_lepton.vec = no_fsr_lepton.vec - photon_vec;
      zzleptons_nofsr.push_back( no_fsr_lepton );
      ret.Z2leptons_without_fsr.push_back(no_fsr_lepton);
    } 
    else{
      zzleptons_nofsr.push_back( ZCandidatesLeptons[best_Z2_index].lepton_pair.second); 
      ret.Z2leptons_without_fsr.push_back( ZCandidatesLeptons[best_Z2_index].lepton_pair.second);
    }

    bool resonance = false;
    for( int i=0; i<zzleptons_nofsr.size(); i++ ) {
      for( int j=i+1; j<zzleptons_nofsr.size(); j++ ) {
        if( zzleptons_nofsr[i].charge == zzleptons_nofsr[j].charge ) continue; // 4/4
        
        if( (zzleptons_nofsr[i].vec+zzleptons_nofsr[j].vec).M() <= 4. ) {
          resonance = true;
          break;
        }
      }
    }
    if( !resonance ) { 
      ret.status.selectionBits.flip(PASS_RESONANCE);
      increment(counts,"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].lepton_pair.first.vec) +
      (ZCandidatesLeptons[best_Z1_index].lepton_pair.second.vec) +
      (ZCandidatesLeptons[best_Z2_index].lepton_pair.first.vec) +
      (ZCandidatesLeptons[best_Z2_index].lepton_pair.second.vec);
    
    if(ctrl.debug)
      cout << "forming zz from: "
           << setw(9) << (ZCandidatesLeptons[best_Z1_index].lepton_pair.first.vec).Pt()
           << setw(9) << (ZCandidatesLeptons[best_Z1_index].lepton_pair.second.vec).Pt()
           << setw(9) << (ZCandidatesLeptons[best_Z2_index].lepton_pair.first.vec).Pt()
           << setw(9) << (ZCandidatesLeptons[best_Z2_index].lepton_pair.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(counts,"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(counts,"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(counts,"m4l>100",theZ1type,theZ2type);
    } else {
      ret.status.setStatus(SelectionStatus::FAIL);
      //return ret;
      // cout << "NOTE: failed m4l>100 (" << zzvec.M() << "), but saving to ntuple anyway" << endl;
    }

    //***************************************************************
    // finish
    //***************************************************************
    
    TLorentzVector theZ1 = (ZCandidatesLeptons[best_Z1_index].lepton_pair.first.vec) + 
        (ZCandidatesLeptons[best_Z1_index].lepton_pair.second.vec); 
    TLorentzVector theZ2 = (ZCandidatesLeptons[best_Z2_index].lepton_pair.first.vec) + 
        (ZCandidatesLeptons[best_Z2_index].lepton_pair.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.25
Committed:	Sun Jan 20 06:55:17 2013 UTC (12 years, 3 months ago) by anlevin
Content type:	text/plain
Branch:	MAIN
Changes since 1.24:	+2 -1 lines
Log Message:	added mass errors calculation