Selection/src/Selection.cc

#include "Selection.h"
#include "PassHLT.h"
#include "HZZCiCElectronSelection.h"
#include "HZZLikelihoodElectronSelection.h"
#include "HZZBDTElectronSelection.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=0;
    CICStruct ciccuts = getCiCCuts(ctrl.eleSeleScheme); 
    unsigned  failsCIC=0;
    if(ctrl.eleSele=="cic") {
      failsCIC = failsCicSelection(ctrl, ele, ciccuts, ctrl.kinematics);
      FAIL = failsCIC;
    }
    LikStruct likcuts;
    unsigned failsLike=0;
    if(ctrl.eleSele=="lik") {
      likcuts = getLikCuts(ctrl.eleSeleScheme);
      failsLike = failsLikelihoodSelection(ele, likcuts, ctrl.kinematics);
      FAIL = failsLike;
    }
    unsigned  failsBDT=0;
    if(ctrl.eleSele=="bdt") {
      failsBDT = failsBDTSelection(ele,ctrl.kinematics,ctrl.eleSeleScheme);
      FAIL = failsBDT;
    }

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