Selection/src/Selection.cc

#include "Selection.h"
#include "PassHLT.h"
#include "HZZCiCElectronSelection.h"
#include "HZZLikelihoodElectronSelection.h"
#include "RunLumiRangeMap.h"

RunLumiRangeMap rlrm;

void initRunLumiRangeMap() {
  rlrm.AddJSONFile(std::string("./data/Cert_136033-149442_7TeV_Apr21ReReco_Collisions10_JSON.txt"));
  rlrm.AddJSONFile(std::string("./data/Cert_160404-173244_7TeV_PromptReco_Collisions11_JSON_v2.txt"));
  rlrm.AddJSONFile(std::string("./data/Cert_160404-163869_7TeV_May10ReReco_Collisions11_JSON_v3.txt"));  
  rlrm.AddJSONFile(std::string("./data/Cert_170249-172619_7TeV_ReReco5Aug_Collisions11_JSON.txt"));  
};

unsigned fails_HZZ4L_selection(ControlFlags &ctrl,           // input control 
                               mithep::TEventInfo *info,     // input event inof
                               TClonesArray *electronArr,    // input electrons
                               TClonesArray *muonArr,        // input muons
                               double eventweight,           // weight
                               TNtuple * passtuple ) { 

  fails_HZZ4L_selection( ctrl, info, electronArr, muonArr, eventweight, passtuple, NULL ); 

};

unsigned fails_HZZ4L_selection(ControlFlags &ctrl,           // input control 
                               mithep::TEventInfo *info,     // input event inof
                               TClonesArray *electronArr,    // input electrons
                               TClonesArray *muonArr,        // input muons
                               double eventweight,           // weight
                               LabVectors  *l ) { 

  fails_HZZ4L_selection( ctrl, info, electronArr, muonArr, eventweight, NULL, l ); 

};


unsigned fails_HZZ4L_selection(ControlFlags &ctrl,           // input control 
                               mithep::TEventInfo *info,     // input event inof
                               TClonesArray *electronArr,    // input electrons
                               TClonesArray *muonArr,        // input muons
                               double eventweight,           // weight
                               TNtuple * passtuple,
                               LabVectors * l) {       // output ntuple 

  unsigned evtfail = 0x0;

  if( !ctrl.mc ) {
    // not accounting for overlap atm
    RunLumiRangeMap::RunLumiPairType rl(info->runNum, info->lumiSec);      
    if( !(rlrm.HasRunLumi(rl)) )  {
      evtfail |= (1<<EVTFAIL_JSON);
      return evtfail;
    }
  }
  
  
  if( ctrl.debug ) {
    cout << "Run: " << info->runNum 
         << "\tEvt: " << info->evtNum
         << "\tLumi: " << info->lumiSec
         << endl;
  }
  
  
  //********************************************************
  // Trigger
  //********************************************************
  if( !ctrl.mc ) { 
    //  if( !(passHLT(info->triggerBits, info->runNum, channel) )  ) { 
    if( !(passHLT(info->triggerBits, info->runNum, 999) )  ) { 
      evtfail |= (1<<EVTFAIL_TRIGGER);
      return evtfail;
    }   
  } else { 
    if( !(passHLTMC(info->triggerBits)) ) { 
      evtfail |= (1<<EVTFAIL_TRIGGER);
      return evtfail;
    }       
    //    cout << "MC trigger bits: " << hex << info->triggerBits << dec << endl;
  }

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

  //********************************************************
  // Lepton Selection
  //********************************************************
  vector<SimpleLepton> lepvec;
  
  //    
  if( ctrl.debug ) cout << "\tnMuons: " << muonArr->GetEntries() << endl;
  //----------------------------------------------------
  for(Int_t i=0; i<muonArr->GetEntries(); i++) {
    const mithep::TMuon *mu = (mithep::TMuon*)((*muonArr)[i]);      
    unsigned muonfail = passMuonSelectionZZ(mu);
    if( ctrl.debug ) { 
      cout << "muon:: pt: " << mu->pt 
           << "\teta: " << mu->eta 
           << "\tmask: 0x" << hex << muonfail << dec 
           << endl;
    }
    if ( !muonfail ) {
      SimpleLepton tmplep;
      tmplep.vec.SetPtEtaPhiM(mu->pt, 
                              mu->eta, 
                              mu->phi,
                              105.658369e-3);
      tmplep.type    = 13;
      tmplep.index   = i;
      tmplep.charge  = mu->q;
      tmplep.isoTrk  = mu->trkIso03;
      tmplep.isoEcal = mu->emIso03;
      tmplep.isoHcal = mu->hadIso03;
      tmplep.ip3dSig = mu->ip3dSig;
      tmplep.is4l    = false;
      tmplep.isEB    = (fabs(mu->eta) < 1.479 ? 1 : 0 ); 
      lepvec.push_back(tmplep);
      if( ctrl.debug ) { cout << "muon passes ... " << endl;}
    }
  }
  
  if( ctrl.debug ) { cout << "\tnElectron: " << electronArr->GetEntries() << endl; }

  //----------------------------------------------------
  for(Int_t i=0; i<electronArr->GetEntries(); i++) {
    const mithep::TElectron *ele = (mithep::TElectron*)((*electronArr)[i]);

    Bool_t isMuonOverlap = kFALSE;
    for (int k=0; k<lepvec.size(); ++k) {
      TVector3 tmplep;
      tmplep.SetPtEtaPhi(ele->pt, ele->eta, ele->phi);
      if ( lepvec[k].type == 13 && lepvec[k].vec.Vect().DrEtaPhi(tmplep) < 0.1 ) {
        if( ctrl.debug ) cout << "-----> isMuonOverlap! " << endl;
        isMuonOverlap = kTRUE; 
        break;
      }
    }

    unsigned FAIL;
    CICStruct tightcuts = getTightCuts(); 
    unsigned  failsCIC  = failsCicSelection(ctrl, ele, tightcuts, ctrl.kinematics);
    LikStruct likcuts = getLikCuts(ctrl.lik_cut);
    unsigned  failsLike = failsLikelihoodSelection(ele, likcuts, ctrl.kinematics);

    if ( ctrl.cic ) { FAIL = failsCIC;}
    else            { FAIL = failsLike; }

    if( ctrl.debug ){ 
      cout << "CIC category: " << cicCategory(ele) 
           << "\tlikelihood: " << ele->likelihood
           << "\tFAIL:  0x"     << hex << FAIL      << dec
           << "\tfailsCIC:  0x" << hex << failsCIC  << dec
           << "\tfailsLike: 0x" << hex << failsLike << dec
           << "\tscEt: " << ele->scEt 
           << "\tscEta: " << ele->scEta
           << "\tncluster: " << ele->ncluster
           << endl;
    }
    if ( !FAIL && !isMuonOverlap ) {
      SimpleLepton tmplep;
      tmplep.vec.SetPtEtaPhiM( ele->pt,
                               ele->eta,
                               ele->phi,
                               0.51099892e-3 );
      tmplep.type    = 11;
      tmplep.index   = i;
      tmplep.charge  = ele->q;
      tmplep.isoTrk  = ele->trkIso03;
      tmplep.isoEcal = ele->emIso03;
      tmplep.isoHcal = ele->hadIso03;
      tmplep.ip3dSig = ele->ip3dSig;
      tmplep.is4l    = false;
      tmplep.isEB    = ele->isEB;
      lepvec.push_back(tmplep);
      if( ctrl.debug ) { cout << "\telectron passes ... " << endl; }
    }
  }
  
  sort( lepvec.begin(), lepvec.end(), SimpleLepton::lep_pt_sort );

  int nmu=0, nele=0;
  for( int i=0; i<lepvec.size(); i++ ) {
    if( abs(lepvec[i].type) == 11 ) nele++;
    else nmu++;
  }
  if( ctrl.debug ) { 
    cout << "postsel nlep: " << lepvec.size()
         << "\tnmuon: " << nmu
         << "\tnelectron: " << nele
         << endl;
  }

  //******************************************************************************
  //Z1 Selection
  //******************************************************************************
  int Z1LeptonPlusIndex = -1;
  int Z1LeptonMinusIndex = -1;
  double BestZ1Mass = -999;
  if( ctrl.debug ) { cout << "looking for a Z1 ..." << endl; }
  for(int i = 0; i < lepvec.size(); ++i) {
    for(int j = i+1; j < lepvec.size(); ++j) {
      if( ctrl.debug ) { cout << "\tconsidering leptons " << i << " & " << j << endl; }
      if (!(lepvec[i].vec.Pt() > 20.0 || lepvec[j].vec.Pt() > 20.0)) continue; 
      if( ctrl.debug ) { cout << "\tat least one is > 20 GeV" << endl; }
      if (!(lepvec[i].vec.Pt() > 10.0 && lepvec[j].vec.Pt() > 10.0)) continue; 
      if( ctrl.debug ) { cout << "\tthe other  is > 10 GeV" << endl; }
      if (lepvec[i].charge == lepvec[j].charge) continue;          
      if( ctrl.debug ) { cout << "\tthey're opposite charge" << endl; }
      if (fabs(lepvec[i].type) != fabs(lepvec[j].type)) continue;  
      if( ctrl.debug ) { cout << "\tthey're same flavor" << endl; }
      
      //Make Z1 hypothesis
      TLorentzVector leptonPlus, leptonMinus;
      if ( lepvec[i].charge > 0 ) {
        leptonPlus  = lepvec[i].vec;
        leptonMinus = lepvec[j].vec;
      } else {
        leptonPlus  = lepvec[j].vec;
        leptonMinus = lepvec[i].vec;
      }
      
      float tmpZ1Mass = (leptonPlus+leptonMinus).M();
      if( ctrl.debug ) cout << "Z1 selection, tmpZ1Mass: " << tmpZ1Mass << endl;
      if( tmpZ1Mass > 60 ) {
        if (fabs(tmpZ1Mass - 91.1876) < fabs(BestZ1Mass - 91.1876)) {
          BestZ1Mass = tmpZ1Mass;
          if( ctrl.debug ) cout << "Z1 selection, new BestZ1Mass: " << BestZ1Mass 
               << "\tdM: " << fabs(BestZ1Mass - 91.1876)
               << endl;
          if (lepvec[i].charge > 0) {
            Z1LeptonPlusIndex = i;
            Z1LeptonMinusIndex = j;
          } else {
            Z1LeptonPlusIndex = j;
            Z1LeptonMinusIndex = i;
          }
        }
      }
    }
  }
  // stop if no Z1 candidate is found
  if( BestZ1Mass < 0 ) { 
    evtfail |= (1<<EVTFAIL_Z1);
    return evtfail; 
  }
  if( ctrl.debug ) cout << "\tgot a Z1 ... run: " << info->runNum << "\tevt: " << info->evtNum  << endl;
  if( ctrl.debug ) cout << "\tZ1 plusindex: " << Z1LeptonPlusIndex << "\tminusindex: " << Z1LeptonMinusIndex << endl;
  TLorentzVector Z1LeptonPlus  = lepvec[Z1LeptonPlusIndex].vec;
  TLorentzVector Z1LeptonMinus = lepvec[Z1LeptonMinusIndex].vec;
  TLorentzVector Z1Candidate   =  Z1LeptonPlus + Z1LeptonMinus;
  if( l != NULL ) { 
    l->vecz1 = Z1Candidate; 
    l->vecl1p = Z1LeptonPlus; 
    l->vecl1m = Z1LeptonMinus; 
  }
  
  //******************************************************************************
  // Z1 + l
  //******************************************************************************
  if( lepvec.size() < 3 ) { 
    evtfail |= (1<<EVTFAIL_Z1_PLUSL);
    return evtfail; 
  }
        
  //******************************************************************************
  // 4l/Z2 Selection
  //******************************************************************************
  Int_t Z2LeptonPlusIndex = -1;
  Int_t Z2LeptonMinusIndex = -1;
  Double_t BestZ2Mass = -1;
  if( ctrl.debug ) cout << "looking for a Z2 ... out of " << lepvec.size() << " leptons" <<endl;
  for(int i = 0; i < lepvec.size(); ++i) {
    for(int j = i+1; j < lepvec.size(); ++j) {
      //            cout << "i: " << i << "\tj: " << j << endl;
      if (i == Z1LeptonPlusIndex || i == Z1LeptonMinusIndex) { 
        //            cout << "\ti matches a Z1 index, skipping ..." << endl;
        continue; //skip Z1 leptons
      }
      if (j == Z1LeptonPlusIndex || j == Z1LeptonMinusIndex) { 
        //            cout << "\tj matches a Z1 index, skipping ..." << endl;
              continue; //skip Z1 leptons
            }
            if (lepvec[i].charge == lepvec[j].charge) { 
              //              cout << "\ti and j are same sign, skipping ..." << endl;
              continue;         //require opp sign
            }
            if (fabs(lepvec[i].type) != fabs(lepvec[j].type)) { 
              //              cout << "\ti and j are not same flavor, skipping ..." << endl;
              continue; //require same flavor
            }
            
            
            //Make Z2 hypothesis
            TLorentzVector leptonPlus, leptonMinus;
            
            if (lepvec[i].charge > 0 ) {
              leptonPlus  = lepvec[i].vec;
              leptonMinus = lepvec[j].vec;
            } else {
              leptonPlus  = lepvec[j].vec;
              leptonMinus = lepvec[i].vec;
            }

            TLorentzVector dilepton = leptonPlus+leptonMinus;
            TLorentzVector fourLepton = Z1Candidate + dilepton;

            if( ctrl.debug ) cout << "dilepton.M() : " << dilepton.M() << endl;
            if( ctrl.debug ) cout << "fourLepton.M() : " << fourLepton.M() << endl;

            if (!(dilepton.M() > 12.0)) continue;
            if (!(fourLepton.M() > 100.0)) continue;

            
            //for 4e and 4mu, require at least 1 of the other opp sign lepton pairs have mass > 12
            if (fabs(lepvec[i].type) == fabs(lepvec[Z1LeptonPlusIndex].type)) {
              TLorentzVector pair1 = Z1LeptonPlus+leptonMinus;
              TLorentzVector pair2 = Z1LeptonMinus+leptonPlus;
              if( ctrl.debug ) cout << "pair1: " << pair1.M() << "\tpair2: "<< pair2.M() << endl;
              if (!(pair1.M() > 12 || pair2.M() > 12)) continue;
            }
            
            //Disambiguiation is done by choosing the pair with the largest ptMax and largest ptMin
            if (Z2LeptonPlusIndex < 0) {
              if (lepvec[i].charge > 0) {
                Z2LeptonPlusIndex = i;
                Z2LeptonMinusIndex = j;
              } else {
                Z2LeptonPlusIndex = j;
                Z2LeptonMinusIndex = i;
              }
            } else {
              Double_t BestPairPtMax = lepvec[Z2LeptonPlusIndex].vec.Pt();               
              Double_t BestPairPtMin = lepvec[Z2LeptonMinusIndex].vec.Pt(); 
              if (lepvec[Z2LeptonMinusIndex].vec.Pt() > BestPairPtMax) {
                BestPairPtMax = lepvec[Z2LeptonMinusIndex].vec.Pt();
                BestPairPtMin = lepvec[Z2LeptonPlusIndex].vec.Pt();
              }

              Double_t CurrentPairPtMax = lepvec[i].vec.Pt();               
              Double_t CurrentPairPtMin = lepvec[j].vec.Pt(); 
              if (lepvec[j].vec.Pt() > CurrentPairPtMax) {
                CurrentPairPtMax = lepvec[j].vec.Pt();
                CurrentPairPtMin = lepvec[i].vec.Pt();
              }

              if (CurrentPairPtMax > BestPairPtMax) {
                if (lepvec[i].charge > 0) {
                  Z2LeptonPlusIndex = i;
                  Z2LeptonMinusIndex = j;
                } else {
                  Z2LeptonPlusIndex = j;
                  Z2LeptonMinusIndex = i;
                }
              } else if (CurrentPairPtMax  == BestPairPtMax) {
                if (CurrentPairPtMin > BestPairPtMin) {
                  if (lepvec[i].charge > 0) {
                    Z2LeptonPlusIndex = i;
                    Z2LeptonMinusIndex = j;
                  } else {
                    Z2LeptonPlusIndex = j;
                    Z2LeptonMinusIndex = i;
                  }                  
                }
              }
            }            
          }
        }
 
        // stop if no Z2 candidate is found
        if (Z2LeptonPlusIndex == -1) { 
          evtfail |= ( 1<<EVTFAIL_4L );
          return evtfail; 
          //      h_evtfail->Fill( evtfail );
          //      cout << "evtfail: " << hex << evtfail << dec << endl;
          //      continue;
        }
        if( ctrl.debug ) cout << "\tgot a Z2 ..." << endl;
        if( ctrl.debug ) cout << "\tZ2 plusindex: " << Z2LeptonPlusIndex
                              << "\tminusindex: " << Z2LeptonMinusIndex << endl;
        TLorentzVector  Z2LeptonPlus  = lepvec[Z2LeptonPlusIndex].vec;
        TLorentzVector  Z2LeptonMinus = lepvec[Z2LeptonMinusIndex].vec;
        TLorentzVector  Z2Candidate   = Z2LeptonPlus+Z2LeptonMinus;
        TLorentzVector  ZZSystem      = Z1Candidate + Z2Candidate;
        if( l != NULL ) { 
          l->vecz2 = Z2Candidate; 
          l->vecl2p = Z2LeptonPlus; 
          l->vecl2m = Z2LeptonMinus; 
          l->vec4l = ZZSystem; 
        }
        lepvec[Z1LeptonPlusIndex].is4l = true;
        lepvec[Z1LeptonMinusIndex].is4l = true;
        lepvec[Z2LeptonPlusIndex].is4l = true;
        lepvec[Z2LeptonMinusIndex].is4l = true;

        //***************************************************************
        // Isolation
        //***************************************************************
        bool failiso=false;

        /*
        int i,j;
        i=Z1LeptonPlusIndex;
        j=Z1LeptonMinusIndex;
        float RIso1 = (lepvec[i].isoTrk+lepvec[i].isoEcal+lepvec[i].isoHcal)/lepvec[i].vec.Pt();
        float RIso2 = (lepvec[j].isoTrk+lepvec[j].isoEcal+lepvec[j].isoHcal)/lepvec[j].vec.Pt();
        float comboIso12 = RIso1 + RIso2;
        i=Z2LeptonPlusIndex;
        j=Z2LeptonMinusIndex;
        float RIso3 = (lepvec[i].isoTrk+lepvec[i].isoEcal+lepvec[i].isoHcal)/lepvec[i].vec.Pt();
        float RIso4 = (lepvec[j].isoTrk+lepvec[j].isoEcal+lepvec[j].isoHcal)/lepvec[j].vec.Pt();
        float comboIso34 = RIso3 + RIso4;
        if( comboIso12 > 0.35 || comboIso34 > 0.35 ) { 
              failiso = true;
        }
        */

        float rho = info->rho;
        for( int i=0; i<lepvec.size(); i++ ) { 
          if( !(lepvec[i].is4l) ) continue;
            float effArea_ecal_i, effArea_hcal_i;
            if( lepvec[i].isEB ) { 
              if( lepvec[i].type == 11 ) {
                effArea_ecal_i = 0.101;
                effArea_hcal_i = 0.021;
              } else { 
                effArea_ecal_i = 0.074;
                effArea_hcal_i = 0.022;
              }
            } else {
              if( lepvec[i].type == 11 ) {
                effArea_ecal_i = 0.046;
                effArea_hcal_i = 0.040;
              } else {
                effArea_ecal_i = 0.045;
                effArea_hcal_i = 0.030;
              }
            }
            float isoEcal_corr_i = lepvec[i].isoEcal - (effArea_ecal_i*rho);
            float isoHcal_corr_i = lepvec[i].isoHcal - (effArea_hcal_i*rho);
            for( int j=i+1; j<lepvec.size(); j++ ) { 
              if( !(lepvec[j].is4l) ) continue;
              float effArea_ecal_j, effArea_hcal_j;
              if( lepvec[j].isEB ) { 
                if( lepvec[j].type == 11 ) {
                  effArea_ecal_j = 0.101;
                  effArea_hcal_j = 0.021;
                } else { 
                  effArea_ecal_j = 0.074;
                  effArea_hcal_j = 0.022;
                }
              } else {
                if( lepvec[j].type == 11 ) {
                  effArea_ecal_j = 0.046;
                  effArea_hcal_j = 0.040;
                } else {
                  effArea_ecal_j = 0.045;
                  effArea_hcal_j = 0.030;
                }
              }
              float isoEcal_corr_j = lepvec[j].isoEcal - (effArea_ecal_j*rho);
              float isoHcal_corr_j = lepvec[j].isoHcal - (effArea_hcal_j*rho);
              float RIso_i = (lepvec[i].isoTrk+isoEcal_corr_i+isoHcal_corr_i)/lepvec[i].vec.Pt();
              float RIso_j = (lepvec[j].isoTrk+isoEcal_corr_j+isoHcal_corr_j)/lepvec[j].vec.Pt();           
              float comboIso = RIso_i + RIso_j; 
              if( info->evtNum == 1038911933 ) {
                float tmpdR = lepvec[i].vec.DrEtaPhi(lepvec[j].vec);
                cout << "i: " << i 
                     << "\tdR: " << tmpdR
                     << "\trho: " << rho
                     << "\tRIso_i: " << RIso_i  
                     << "\ttkrel: " << lepvec[i].isoTrk/lepvec[i].vec.Pt() 
                     << "\tecalrel: " << lepvec[i].isoEcal/lepvec[i].vec.Pt() 
                     << "\tecalrelcor: " << isoEcal_corr_i/lepvec[i].vec.Pt() 
                     << "\thcalrel: " << lepvec[i].isoHcal/lepvec[i].vec.Pt() 
                     << "\thcalrelcor: " << isoHcal_corr_i/lepvec[i].vec.Pt() 
                     << "\tpt_i: " << lepvec[i].vec.Pt()
                     << "\tj: " << j 
                     << "\tRIso_j: " << RIso_j  
                     << "\ttkrel: "  << lepvec[j].isoTrk/lepvec[j].vec.Pt() 
                     << "\tecalrel: " << lepvec[j].isoEcal/lepvec[j].vec.Pt() 
                     << "\tecalrelcor: " << isoEcal_corr_j/lepvec[j].vec.Pt() 
                     << "\thcalrel: " << lepvec[j].isoHcal/lepvec[j].vec.Pt() 
                     << "\thcalrelcor: " << isoHcal_corr_j/lepvec[j].vec.Pt() 
                     << "\tpt_j: " << lepvec[j].vec.Pt()
                     << "\tcombo: " << comboIso 
                     << endl;
                cout.flush();
              }
              if( comboIso > 0.35 ) { 
                if( ctrl.debug ) cout << "combo failing for indices: " << i << "," << j << endl;
                failiso = true;
                //            break;
              }
            }
        }
        if( failiso ) { 
          evtfail |= ( 1<<EVTFAIL_ISOLATION );
          return evtfail; 
          //h_evtfail->Fill( evtfail, eventweight );
          //      h_evtfail->Fill( evtfail );
          //      cout << "evtfail: " << hex << evtfail << dec << endl;
          //      continue;
        }

        //***************************************************************
        // IP significance
        //***************************************************************
        bool failip = false;
        for( int i=0; i<lepvec.size(); i++ ) { 
          if( !(lepvec[i].is4l) ) continue;
          if( lepvec[i].ip3dSig > 4 ) { 
            failip=true;
            break;
          }
        }
        if( failip ) {
          evtfail |= (1<<EVTFAIL_IP );
          return evtfail; 
          //h_evtfail->Fill( evtfail, eventweight );
          //      h_evtfail->Fill( evtfail );
          //      cout << "evtfail: " << hex << evtfail << dec << endl;
          //      continue;
        }

        //***************************************************************
        // remaining kinematic cuts 
        //***************************************************************
        double Z2massCut=0;
        if      ( ctrl.kinematics == "loose" ) Z2massCut = 12;
        else if ( ctrl.kinematics == "tight" ) Z2massCut = 20;
        else { cout <<  "error! kinematic tightness not defined!" << endl; assert(0); }

        if ( Z1Candidate.M() > 120        || 
             Z2Candidate.M() < Z2massCut  ||
             Z2Candidate.M() > 120        || 
             !(lepvec[Z1LeptonPlusIndex].vec.Pt() > 20.0 || lepvec[Z1LeptonMinusIndex].vec.Pt() > 20.0) ||
             !(lepvec[Z1LeptonPlusIndex].vec.Pt() > 10.0 && lepvec[Z1LeptonMinusIndex].vec.Pt() > 10.0)
             ) {
          evtfail |= (1<<EVTFAIL_KINEMATICS );
          return evtfail; 
          //h_evtfail->Fill( evtfail, eventweight );
          //      h_evtfail->Fill( evtfail );
          //      cout << "evtfail: " << hex << evtfail << dec << endl;
          //      continue;
        }
        
        int channel;
        if( lepvec[Z1LeptonMinusIndex].type == 11 && lepvec[Z2LeptonMinusIndex].type == 11 ) channel=0; 
        if( lepvec[Z1LeptonMinusIndex].type == 13 && lepvec[Z2LeptonMinusIndex].type == 13 ) channel=1; 
        if( (lepvec[Z1LeptonMinusIndex].type == 11 && lepvec[Z2LeptonMinusIndex].type == 13) || 
            (lepvec[Z1LeptonMinusIndex].type == 13 && lepvec[Z2LeptonMinusIndex].type == 11)) channel=2; 
        

        if( passtuple != NULL ) { 
          passtuple->Fill( info->runNum, 
                           info->evtNum,
                           info->lumiSec, 
                           channel, 
                           Z1Candidate.M(),
                           Z2Candidate.M(),
                           ZZSystem.M(),
                           ZZSystem.Pt(),
                           eventweight);
        }

        if( ctrl.debug ) cout  << "run: " << info->runNum   
              << "\tevt: " << info->evtNum
              << "\tZ1channel: " << lepvec[Z1LeptonMinusIndex].type
              << "\tZ2channel: " << lepvec[Z2LeptonMinusIndex].type
              << "\tmZ1: " << Z1Candidate.M()
              << "\tmZ2: " << Z2Candidate.M()
              << "\tm4l: " << ZZSystem.M()
              << "\tevtfail: " << hex << evtfail << dec 
              << "\ttrigbits: " << hex << info->triggerBits << dec
          //          << "\ttree: " << inputFiles[q][f] 
              << endl;

        return evtfail;
          
}


Revision:	1.3
Committed:	Thu Sep 22 18:01:19 2011 UTC (13 years, 7 months ago) by dkralph
Content type:	text/plain
Branch:	MAIN
Changes since 1.2:	+41 -23 lines
Log Message:	Added likelihood electron ID.
#	User	Rev	Content
1	khahn	1.1	#include "Selection.h"
2			#include "PassHLT.h"
3			#include "HZZCiCElectronSelection.h"
4	dkralph	1.3	#include "HZZLikelihoodElectronSelection.h"
5	khahn	1.1	#include "RunLumiRangeMap.h"
6
7	dkralph	1.2	RunLumiRangeMap rlrm;
8	khahn	1.1
9			void initRunLumiRangeMap() {
10			rlrm.AddJSONFile(std::string("./data/Cert_136033-149442_7TeV_Apr21ReReco_Collisions10_JSON.txt"));
11			rlrm.AddJSONFile(std::string("./data/Cert_160404-173244_7TeV_PromptReco_Collisions11_JSON_v2.txt"));
12			rlrm.AddJSONFile(std::string("./data/Cert_160404-163869_7TeV_May10ReReco_Collisions11_JSON_v3.txt"));
13			rlrm.AddJSONFile(std::string("./data/Cert_170249-172619_7TeV_ReReco5Aug_Collisions11_JSON.txt"));
14			};
15
16			unsigned fails_HZZ4L_selection(ControlFlags &ctrl, // input control
17			mithep::TEventInfo *info, // input event inof
18			TClonesArray *electronArr, // input electrons
19			TClonesArray *muonArr, // input muons
20	dkralph	1.2	double eventweight, // weight
21	khahn	1.1	TNtuple * passtuple ) {
22
23	dkralph	1.2	fails_HZZ4L_selection( ctrl, info, electronArr, muonArr, eventweight, passtuple, NULL );
24	khahn	1.1
25			};
26
27			unsigned fails_HZZ4L_selection(ControlFlags &ctrl, // input control
28			mithep::TEventInfo *info, // input event inof
29			TClonesArray *electronArr, // input electrons
30			TClonesArray *muonArr, // input muons
31	dkralph	1.2	double eventweight, // weight
32	khahn	1.1	LabVectors *l ) {
33
34	dkralph	1.2	fails_HZZ4L_selection( ctrl, info, electronArr, muonArr, eventweight, NULL, l );
35	khahn	1.1
36			};
37
38
39			unsigned fails_HZZ4L_selection(ControlFlags &ctrl, // input control
40			mithep::TEventInfo *info, // input event inof
41			TClonesArray *electronArr, // input electrons
42			TClonesArray *muonArr, // input muons
43	dkralph	1.2	double eventweight, // weight
44	khahn	1.1	TNtuple * passtuple,
45	dkralph	1.3	LabVectors * l) { // output ntuple
46	khahn	1.1
47			unsigned evtfail = 0x0;
48
49			if( !ctrl.mc ) {
50			// not accounting for overlap atm
51	dkralph	1.2	RunLumiRangeMap::RunLumiPairType rl(info->runNum, info->lumiSec);
52	khahn	1.1	if( !(rlrm.HasRunLumi(rl)) ) {
53			evtfail \|= (1<<EVTFAIL_JSON);
54			return evtfail;
55			}
56			}
57
58
59			if( ctrl.debug ) {
60			cout << "Run: " << info->runNum
61			<< "\tEvt: " << info->evtNum
62			<< "\tLumi: " << info->lumiSec
63			<< endl;
64			}
65
66
67			//********************************************************
68			// Trigger
69			//********************************************************
70			if( !ctrl.mc ) {
71			// if( !(passHLT(info->triggerBits, info->runNum, channel) ) ) {
72			if( !(passHLT(info->triggerBits, info->runNum, 999) ) ) {
73			evtfail \|= (1<<EVTFAIL_TRIGGER);
74			return evtfail;
75			}
76			} else {
77			if( !(passHLTMC(info->triggerBits)) ) {
78			evtfail \|= (1<<EVTFAIL_TRIGGER);
79			return evtfail;
80			}
81			// cout << "MC trigger bits: " << hex << info->triggerBits << dec << endl;
82			}
83
84			if( ctrl.debug ) {
85			cout << "presel nlep: " << muonArr->GetEntries() + electronArr->GetEntries()
86			<< "\tnmuon: " << muonArr->GetEntries()
87			<< "\tnelectron: " << electronArr->GetEntries()
88			<< endl;
89			}
90
91			//********************************************************
92			// Lepton Selection
93			//********************************************************
94			vector<SimpleLepton> lepvec;
95
96			//
97	dkralph	1.3	if( ctrl.debug ) cout << "\tnMuons: " << muonArr->GetEntries() << endl;
98	khahn	1.1	//----------------------------------------------------
99			for(Int_t i=0; i<muonArr->GetEntries(); i++) {
100			const mithep::TMuon mu = (mithep::TMuon)((*muonArr)[i]);
101			unsigned muonfail = passMuonSelectionZZ(mu);
102			if( ctrl.debug ) {
103			cout << "muon:: pt: " << mu->pt
104			<< "\teta: " << mu->eta
105			<< "\tmask: 0x" << hex << muonfail << dec
106			<< endl;
107			}
108			if ( !muonfail ) {
109			SimpleLepton tmplep;
110			tmplep.vec.SetPtEtaPhiM(mu->pt,
111			mu->eta,
112			mu->phi,
113			105.658369e-3);
114			tmplep.type = 13;
115			tmplep.index = i;
116			tmplep.charge = mu->q;
117			tmplep.isoTrk = mu->trkIso03;
118			tmplep.isoEcal = mu->emIso03;
119			tmplep.isoHcal = mu->hadIso03;
120			tmplep.ip3dSig = mu->ip3dSig;
121			tmplep.is4l = false;
122			tmplep.isEB = (fabs(mu->eta) < 1.479 ? 1 : 0 );
123			lepvec.push_back(tmplep);
124			if( ctrl.debug ) { cout << "muon passes ... " << endl;}
125			}
126			}
127
128			if( ctrl.debug ) { cout << "\tnElectron: " << electronArr->GetEntries() << endl; }
129	dkralph	1.3
130	khahn	1.1	//----------------------------------------------------
131			for(Int_t i=0; i<electronArr->GetEntries(); i++) {
132			const mithep::TElectron ele = (mithep::TElectron)((*electronArr)[i]);
133
134			Bool_t isMuonOverlap = kFALSE;
135			for (int k=0; k<lepvec.size(); ++k) {
136			TVector3 tmplep;
137			tmplep.SetPtEtaPhi(ele->pt, ele->eta, ele->phi);
138			if ( lepvec[k].type == 13 && lepvec[k].vec.Vect().DrEtaPhi(tmplep) < 0.1 ) {
139	dkralph	1.3	if( ctrl.debug ) cout << "-----> isMuonOverlap! " << endl;
140	khahn	1.1	isMuonOverlap = kTRUE;
141			break;
142			}
143			}
144
145	dkralph	1.3	unsigned FAIL;
146	khahn	1.1	CICStruct tightcuts = getTightCuts();
147	dkralph	1.3	unsigned failsCIC = failsCicSelection(ctrl, ele, tightcuts, ctrl.kinematics);
148			LikStruct likcuts = getLikCuts(ctrl.lik_cut);
149			unsigned failsLike = failsLikelihoodSelection(ele, likcuts, ctrl.kinematics);
150
151			if ( ctrl.cic ) { FAIL = failsCIC;}
152			else { FAIL = failsLike; }
153
154	khahn	1.1	if( ctrl.debug ){
155			cout << "CIC category: " << cicCategory(ele)
156	dkralph	1.3	<< "\tlikelihood: " << ele->likelihood
157			<< "\tFAIL: 0x" << hex << FAIL << dec
158			<< "\tfailsCIC: 0x" << hex << failsCIC << dec
159			<< "\tfailsLike: 0x" << hex << failsLike << dec
160	khahn	1.1	<< "\tscEt: " << ele->scEt
161	dkralph	1.3	<< "\tscEta: " << ele->scEta
162			<< "\tncluster: " << ele->ncluster
163	khahn	1.1	<< endl;
164			}
165	dkralph	1.3	if ( !FAIL && !isMuonOverlap ) {
166	khahn	1.1	SimpleLepton tmplep;
167			tmplep.vec.SetPtEtaPhiM( ele->pt,
168			ele->eta,
169			ele->phi,
170			0.51099892e-3 );
171			tmplep.type = 11;
172			tmplep.index = i;
173			tmplep.charge = ele->q;
174			tmplep.isoTrk = ele->trkIso03;
175			tmplep.isoEcal = ele->emIso03;
176			tmplep.isoHcal = ele->hadIso03;
177			tmplep.ip3dSig = ele->ip3dSig;
178			tmplep.is4l = false;
179			tmplep.isEB = ele->isEB;
180			lepvec.push_back(tmplep);
181			if( ctrl.debug ) { cout << "\telectron passes ... " << endl; }
182			}
183			}
184
185			sort( lepvec.begin(), lepvec.end(), SimpleLepton::lep_pt_sort );
186
187			int nmu=0, nele=0;
188			for( int i=0; i<lepvec.size(); i++ ) {
189			if( abs(lepvec[i].type) == 11 ) nele++;
190			else nmu++;
191			}
192			if( ctrl.debug ) {
193			cout << "postsel nlep: " << lepvec.size()
194			<< "\tnmuon: " << nmu
195			<< "\tnelectron: " << nele
196			<< endl;
197			}
198
199			//******************************************************************************
200			//Z1 Selection
201			//******************************************************************************
202			int Z1LeptonPlusIndex = -1;
203			int Z1LeptonMinusIndex = -1;
204			double BestZ1Mass = -999;
205			if( ctrl.debug ) { cout << "looking for a Z1 ..." << endl; }
206			for(int i = 0; i < lepvec.size(); ++i) {
207			for(int j = i+1; j < lepvec.size(); ++j) {
208			if( ctrl.debug ) { cout << "\tconsidering leptons " << i << " & " << j << endl; }
209			if (!(lepvec[i].vec.Pt() > 20.0 \|\| lepvec[j].vec.Pt() > 20.0)) continue;
210			if( ctrl.debug ) { cout << "\tat least one is > 20 GeV" << endl; }
211			if (!(lepvec[i].vec.Pt() > 10.0 && lepvec[j].vec.Pt() > 10.0)) continue;
212			if( ctrl.debug ) { cout << "\tthe other is > 10 GeV" << endl; }
213			if (lepvec[i].charge == lepvec[j].charge) continue;
214			if( ctrl.debug ) { cout << "\tthey're opposite charge" << endl; }
215			if (fabs(lepvec[i].type) != fabs(lepvec[j].type)) continue;
216			if( ctrl.debug ) { cout << "\tthey're same flavor" << endl; }
217
218			//Make Z1 hypothesis
219			TLorentzVector leptonPlus, leptonMinus;
220			if ( lepvec[i].charge > 0 ) {
221			leptonPlus = lepvec[i].vec;
222			leptonMinus = lepvec[j].vec;
223			} else {
224			leptonPlus = lepvec[j].vec;
225			leptonMinus = lepvec[i].vec;
226			}
227
228			float tmpZ1Mass = (leptonPlus+leptonMinus).M();
229	dkralph	1.3	if( ctrl.debug ) cout << "Z1 selection, tmpZ1Mass: " << tmpZ1Mass << endl;
230	khahn	1.1	if( tmpZ1Mass > 60 ) {
231			if (fabs(tmpZ1Mass - 91.1876) < fabs(BestZ1Mass - 91.1876)) {
232			BestZ1Mass = tmpZ1Mass;
233	dkralph	1.3	if( ctrl.debug ) cout << "Z1 selection, new BestZ1Mass: " << BestZ1Mass
234	khahn	1.1	<< "\tdM: " << fabs(BestZ1Mass - 91.1876)
235			<< endl;
236			if (lepvec[i].charge > 0) {
237			Z1LeptonPlusIndex = i;
238			Z1LeptonMinusIndex = j;
239			} else {
240			Z1LeptonPlusIndex = j;
241			Z1LeptonMinusIndex = i;
242			}
243			}
244			}
245			}
246			}
247			// stop if no Z1 candidate is found
248			if( BestZ1Mass < 0 ) {
249			evtfail \|= (1<<EVTFAIL_Z1);
250			return evtfail;
251			}
252	dkralph	1.3	if( ctrl.debug ) cout << "\tgot a Z1 ... run: " << info->runNum << "\tevt: " << info->evtNum << endl;
253			if( ctrl.debug ) cout << "\tZ1 plusindex: " << Z1LeptonPlusIndex << "\tminusindex: " << Z1LeptonMinusIndex << endl;
254	khahn	1.1	TLorentzVector Z1LeptonPlus = lepvec[Z1LeptonPlusIndex].vec;
255			TLorentzVector Z1LeptonMinus = lepvec[Z1LeptonMinusIndex].vec;
256			TLorentzVector Z1Candidate = Z1LeptonPlus + Z1LeptonMinus;
257			if( l != NULL ) {
258			l->vecz1 = Z1Candidate;
259			l->vecl1p = Z1LeptonPlus;
260			l->vecl1m = Z1LeptonMinus;
261			}
262
263			//******************************************************************************
264			// Z1 + l
265			//******************************************************************************
266			if( lepvec.size() < 3 ) {
267			evtfail \|= (1<<EVTFAIL_Z1_PLUSL);
268			return evtfail;
269			}
270
271			//******************************************************************************
272			// 4l/Z2 Selection
273			//******************************************************************************
274			Int_t Z2LeptonPlusIndex = -1;
275			Int_t Z2LeptonMinusIndex = -1;
276			Double_t BestZ2Mass = -1;
277	dkralph	1.3	if( ctrl.debug ) cout << "looking for a Z2 ... out of " << lepvec.size() << " leptons" <<endl;
278	khahn	1.1	for(int i = 0; i < lepvec.size(); ++i) {
279			for(int j = i+1; j < lepvec.size(); ++j) {
280			// cout << "i: " << i << "\tj: " << j << endl;
281			if (i == Z1LeptonPlusIndex \|\| i == Z1LeptonMinusIndex) {
282			// cout << "\ti matches a Z1 index, skipping ..." << endl;
283			continue; //skip Z1 leptons
284			}
285			if (j == Z1LeptonPlusIndex \|\| j == Z1LeptonMinusIndex) {
286			// cout << "\tj matches a Z1 index, skipping ..." << endl;
287			continue; //skip Z1 leptons
288			}
289			if (lepvec[i].charge == lepvec[j].charge) {
290			// cout << "\ti and j are same sign, skipping ..." << endl;
291			continue; //require opp sign
292			}
293			if (fabs(lepvec[i].type) != fabs(lepvec[j].type)) {
294			// cout << "\ti and j are not same flavor, skipping ..." << endl;
295			continue; //require same flavor
296			}
297
298
299			//Make Z2 hypothesis
300			TLorentzVector leptonPlus, leptonMinus;
301
302			if (lepvec[i].charge > 0 ) {
303			leptonPlus = lepvec[i].vec;
304			leptonMinus = lepvec[j].vec;
305			} else {
306			leptonPlus = lepvec[j].vec;
307			leptonMinus = lepvec[i].vec;
308			}
309
310			TLorentzVector dilepton = leptonPlus+leptonMinus;
311			TLorentzVector fourLepton = Z1Candidate + dilepton;
312
313	dkralph	1.3	if( ctrl.debug ) cout << "dilepton.M() : " << dilepton.M() << endl;
314			if( ctrl.debug ) cout << "fourLepton.M() : " << fourLepton.M() << endl;
315	khahn	1.1
316			if (!(dilepton.M() > 12.0)) continue;
317			if (!(fourLepton.M() > 100.0)) continue;
318
319
320			//for 4e and 4mu, require at least 1 of the other opp sign lepton pairs have mass > 12
321			if (fabs(lepvec[i].type) == fabs(lepvec[Z1LeptonPlusIndex].type)) {
322			TLorentzVector pair1 = Z1LeptonPlus+leptonMinus;
323			TLorentzVector pair2 = Z1LeptonMinus+leptonPlus;
324	dkralph	1.3	if( ctrl.debug ) cout << "pair1: " << pair1.M() << "\tpair2: "<< pair2.M() << endl;
325	khahn	1.1	if (!(pair1.M() > 12 \|\| pair2.M() > 12)) continue;
326			}
327
328			//Disambiguiation is done by choosing the pair with the largest ptMax and largest ptMin
329			if (Z2LeptonPlusIndex < 0) {
330			if (lepvec[i].charge > 0) {
331			Z2LeptonPlusIndex = i;
332			Z2LeptonMinusIndex = j;
333			} else {
334			Z2LeptonPlusIndex = j;
335			Z2LeptonMinusIndex = i;
336			}
337			} else {
338			Double_t BestPairPtMax = lepvec[Z2LeptonPlusIndex].vec.Pt();
339			Double_t BestPairPtMin = lepvec[Z2LeptonMinusIndex].vec.Pt();
340			if (lepvec[Z2LeptonMinusIndex].vec.Pt() > BestPairPtMax) {
341			BestPairPtMax = lepvec[Z2LeptonMinusIndex].vec.Pt();
342			BestPairPtMin = lepvec[Z2LeptonPlusIndex].vec.Pt();
343			}
344
345			Double_t CurrentPairPtMax = lepvec[i].vec.Pt();
346			Double_t CurrentPairPtMin = lepvec[j].vec.Pt();
347			if (lepvec[j].vec.Pt() > CurrentPairPtMax) {
348			CurrentPairPtMax = lepvec[j].vec.Pt();
349			CurrentPairPtMin = lepvec[i].vec.Pt();
350			}
351
352			if (CurrentPairPtMax > BestPairPtMax) {
353			if (lepvec[i].charge > 0) {
354			Z2LeptonPlusIndex = i;
355			Z2LeptonMinusIndex = j;
356			} else {
357			Z2LeptonPlusIndex = j;
358			Z2LeptonMinusIndex = i;
359			}
360			} else if (CurrentPairPtMax == BestPairPtMax) {
361			if (CurrentPairPtMin > BestPairPtMin) {
362			if (lepvec[i].charge > 0) {
363			Z2LeptonPlusIndex = i;
364			Z2LeptonMinusIndex = j;
365			} else {
366			Z2LeptonPlusIndex = j;
367			Z2LeptonMinusIndex = i;
368			}
369			}
370			}
371			}
372			}
373			}
374
375			// stop if no Z2 candidate is found
376			if (Z2LeptonPlusIndex == -1) {
377			evtfail \|= ( 1<<EVTFAIL_4L );
378			return evtfail;
379			// h_evtfail->Fill( evtfail );
380			// cout << "evtfail: " << hex << evtfail << dec << endl;
381			// continue;
382			}
383	dkralph	1.3	if( ctrl.debug ) cout << "\tgot a Z2 ..." << endl;
384			if( ctrl.debug ) cout << "\tZ2 plusindex: " << Z2LeptonPlusIndex
385			<< "\tminusindex: " << Z2LeptonMinusIndex << endl;
386	khahn	1.1	TLorentzVector Z2LeptonPlus = lepvec[Z2LeptonPlusIndex].vec;
387			TLorentzVector Z2LeptonMinus = lepvec[Z2LeptonMinusIndex].vec;
388			TLorentzVector Z2Candidate = Z2LeptonPlus+Z2LeptonMinus;
389			TLorentzVector ZZSystem = Z1Candidate + Z2Candidate;
390			if( l != NULL ) {
391			l->vecz2 = Z2Candidate;
392			l->vecl2p = Z2LeptonPlus;
393			l->vecl2m = Z2LeptonMinus;
394			l->vec4l = ZZSystem;
395			}
396			lepvec[Z1LeptonPlusIndex].is4l = true;
397			lepvec[Z1LeptonMinusIndex].is4l = true;
398			lepvec[Z2LeptonPlusIndex].is4l = true;
399			lepvec[Z2LeptonMinusIndex].is4l = true;
400
401			//***************************************************************
402			// Isolation
403			//***************************************************************
404			bool failiso=false;
405
406			/*
407			int i,j;
408			i=Z1LeptonPlusIndex;
409			j=Z1LeptonMinusIndex;
410			float RIso1 = (lepvec[i].isoTrk+lepvec[i].isoEcal+lepvec[i].isoHcal)/lepvec[i].vec.Pt();
411			float RIso2 = (lepvec[j].isoTrk+lepvec[j].isoEcal+lepvec[j].isoHcal)/lepvec[j].vec.Pt();
412			float comboIso12 = RIso1 + RIso2;
413			i=Z2LeptonPlusIndex;
414			j=Z2LeptonMinusIndex;
415			float RIso3 = (lepvec[i].isoTrk+lepvec[i].isoEcal+lepvec[i].isoHcal)/lepvec[i].vec.Pt();
416			float RIso4 = (lepvec[j].isoTrk+lepvec[j].isoEcal+lepvec[j].isoHcal)/lepvec[j].vec.Pt();
417			float comboIso34 = RIso3 + RIso4;
418			if( comboIso12 > 0.35 \|\| comboIso34 > 0.35 ) {
419			failiso = true;
420			}
421			*/
422
423	dkralph	1.2	float rho = info->rho;
424	khahn	1.1	for( int i=0; i<lepvec.size(); i++ ) {
425			if( !(lepvec[i].is4l) ) continue;
426			float effArea_ecal_i, effArea_hcal_i;
427			if( lepvec[i].isEB ) {
428			if( lepvec[i].type == 11 ) {
429			effArea_ecal_i = 0.101;
430			effArea_hcal_i = 0.021;
431			} else {
432			effArea_ecal_i = 0.074;
433			effArea_hcal_i = 0.022;
434			}
435			} else {
436			if( lepvec[i].type == 11 ) {
437			effArea_ecal_i = 0.046;
438			effArea_hcal_i = 0.040;
439			} else {
440			effArea_ecal_i = 0.045;
441			effArea_hcal_i = 0.030;
442			}
443			}
444			float isoEcal_corr_i = lepvec[i].isoEcal - (effArea_ecal_i*rho);
445			float isoHcal_corr_i = lepvec[i].isoHcal - (effArea_hcal_i*rho);
446			for( int j=i+1; j<lepvec.size(); j++ ) {
447			if( !(lepvec[j].is4l) ) continue;
448			float effArea_ecal_j, effArea_hcal_j;
449			if( lepvec[j].isEB ) {
450			if( lepvec[j].type == 11 ) {
451			effArea_ecal_j = 0.101;
452			effArea_hcal_j = 0.021;
453			} else {
454			effArea_ecal_j = 0.074;
455			effArea_hcal_j = 0.022;
456			}
457			} else {
458			if( lepvec[j].type == 11 ) {
459			effArea_ecal_j = 0.046;
460			effArea_hcal_j = 0.040;
461			} else {
462			effArea_ecal_j = 0.045;
463			effArea_hcal_j = 0.030;
464			}
465			}
466			float isoEcal_corr_j = lepvec[j].isoEcal - (effArea_ecal_j*rho);
467			float isoHcal_corr_j = lepvec[j].isoHcal - (effArea_hcal_j*rho);
468			float RIso_i = (lepvec[i].isoTrk+isoEcal_corr_i+isoHcal_corr_i)/lepvec[i].vec.Pt();
469			float RIso_j = (lepvec[j].isoTrk+isoEcal_corr_j+isoHcal_corr_j)/lepvec[j].vec.Pt();
470			float comboIso = RIso_i + RIso_j;
471			if( info->evtNum == 1038911933 ) {
472			float tmpdR = lepvec[i].vec.DrEtaPhi(lepvec[j].vec);
473			cout << "i: " << i
474			<< "\tdR: " << tmpdR
475			<< "\trho: " << rho
476			<< "\tRIso_i: " << RIso_i
477			<< "\ttkrel: " << lepvec[i].isoTrk/lepvec[i].vec.Pt()
478			<< "\tecalrel: " << lepvec[i].isoEcal/lepvec[i].vec.Pt()
479			<< "\tecalrelcor: " << isoEcal_corr_i/lepvec[i].vec.Pt()
480			<< "\thcalrel: " << lepvec[i].isoHcal/lepvec[i].vec.Pt()
481			<< "\thcalrelcor: " << isoHcal_corr_i/lepvec[i].vec.Pt()
482			<< "\tpt_i: " << lepvec[i].vec.Pt()
483			<< "\tj: " << j
484			<< "\tRIso_j: " << RIso_j
485			<< "\ttkrel: " << lepvec[j].isoTrk/lepvec[j].vec.Pt()
486			<< "\tecalrel: " << lepvec[j].isoEcal/lepvec[j].vec.Pt()
487			<< "\tecalrelcor: " << isoEcal_corr_j/lepvec[j].vec.Pt()
488			<< "\thcalrel: " << lepvec[j].isoHcal/lepvec[j].vec.Pt()
489			<< "\thcalrelcor: " << isoHcal_corr_j/lepvec[j].vec.Pt()
490			<< "\tpt_j: " << lepvec[j].vec.Pt()
491			<< "\tcombo: " << comboIso
492			<< endl;
493			cout.flush();
494			}
495			if( comboIso > 0.35 ) {
496	dkralph	1.3	if( ctrl.debug ) cout << "combo failing for indices: " << i << "," << j << endl;
497	khahn	1.1	failiso = true;
498			// break;
499			}
500			}
501			}
502			if( failiso ) {
503			evtfail \|= ( 1<<EVTFAIL_ISOLATION );
504			return evtfail;
505			//h_evtfail->Fill( evtfail, eventweight );
506			// h_evtfail->Fill( evtfail );
507			// cout << "evtfail: " << hex << evtfail << dec << endl;
508			// continue;
509			}
510
511			//***************************************************************
512			// IP significance
513			//***************************************************************
514			bool failip = false;
515			for( int i=0; i<lepvec.size(); i++ ) {
516			if( !(lepvec[i].is4l) ) continue;
517			if( lepvec[i].ip3dSig > 4 ) {
518			failip=true;
519			break;
520			}
521			}
522			if( failip ) {
523			evtfail \|= (1<<EVTFAIL_IP );
524			return evtfail;
525			//h_evtfail->Fill( evtfail, eventweight );
526			// h_evtfail->Fill( evtfail );
527			// cout << "evtfail: " << hex << evtfail << dec << endl;
528			// continue;
529			}
530
531			//***************************************************************
532			// remaining kinematic cuts
533			//***************************************************************
534	dkralph	1.3	double Z2massCut=0;
535			if ( ctrl.kinematics == "loose" ) Z2massCut = 12;
536			else if ( ctrl.kinematics == "tight" ) Z2massCut = 20;
537			else { cout << "error! kinematic tightness not defined!" << endl; assert(0); }
538
539			if ( Z1Candidate.M() > 120 \|\|
540			Z2Candidate.M() < Z2massCut \|\|
541			Z2Candidate.M() > 120 \|\|
542	khahn	1.1	!(lepvec[Z1LeptonPlusIndex].vec.Pt() > 20.0 \|\| lepvec[Z1LeptonMinusIndex].vec.Pt() > 20.0) \|\|
543			!(lepvec[Z1LeptonPlusIndex].vec.Pt() > 10.0 && lepvec[Z1LeptonMinusIndex].vec.Pt() > 10.0)
544			) {
545			evtfail \|= (1<<EVTFAIL_KINEMATICS );
546			return evtfail;
547			//h_evtfail->Fill( evtfail, eventweight );
548			// h_evtfail->Fill( evtfail );
549			// cout << "evtfail: " << hex << evtfail << dec << endl;
550			// continue;
551			}
552
553			int channel;
554			if( lepvec[Z1LeptonMinusIndex].type == 11 && lepvec[Z2LeptonMinusIndex].type == 11 ) channel=0;
555			if( lepvec[Z1LeptonMinusIndex].type == 13 && lepvec[Z2LeptonMinusIndex].type == 13 ) channel=1;
556			if( (lepvec[Z1LeptonMinusIndex].type == 11 && lepvec[Z2LeptonMinusIndex].type == 13) \|\|
557			(lepvec[Z1LeptonMinusIndex].type == 13 && lepvec[Z2LeptonMinusIndex].type == 11)) channel=2;
558
559
560
561			if( passtuple != NULL ) {
562			passtuple->Fill( info->runNum,
563			info->evtNum,
564			info->lumiSec,
565			channel,
566			Z1Candidate.M(),
567			Z2Candidate.M(),
568			ZZSystem.M(),
569			ZZSystem.Pt(),
570			eventweight);
571			}
572
573	dkralph	1.3	if( ctrl.debug ) cout << "run: " << info->runNum
574	khahn	1.1	<< "\tevt: " << info->evtNum
575			<< "\tZ1channel: " << lepvec[Z1LeptonMinusIndex].type
576			<< "\tZ2channel: " << lepvec[Z2LeptonMinusIndex].type
577			<< "\tmZ1: " << Z1Candidate.M()
578			<< "\tmZ2: " << Z2Candidate.M()
579			<< "\tm4l: " << ZZSystem.M()
580			<< "\tevtfail: " << hex << evtfail << dec
581			<< "\ttrigbits: " << hex << info->triggerBits << dec
582			// << "\ttree: " << inputFiles[q][f]
583			<< endl;
584
585			return evtfail;
586
587			}
588
589
590
591
592