Selection/src/Selection.cc

#include "Selection.h"
#include "PassHLT.h"

#include "SiMVAElectronSelection.h"

#include "HZZCiCElectronSelection.h"
#include "HZZLikelihoodElectronSelection.h"
#include "HZZBDTElectronSelection.h"
#include "RunLumiRangeMap.h"
#include "EfficiencyWeightsInterface.h"

RunLumiRangeMap rlrm;
TH1D * hpu;


// #include "TH2D.h"
// extern TH2D * h_mu_eff_idiso_2011a;
// extern TH2D * h_mu_eff_idiso_s11;
// extern TH2D * h_mu_effdata_trigger_doublemu7_2011a;
// extern TH2D * h_mu_effdata_trigger_doublemu13_8_leading_2011a;
// extern TH2D * h_mu_effdata_trigger_doublemu13_8_trailing_2011a;

// unsigned getGenChannel(mithep::TGenInfo * ginfo) { 
//   int gchannel=-1;
//   if( abs(ginfo->id_1_a) == EGenType::kElectron && abs(ginfo->id_1_b) == EGenType::kElectron )  gchannel=0;
//   else if( abs(ginfo->id_1_a) == EGenType::kMuon && abs(ginfo->id_1_b) == EGenType::kMuon ) gchannel=1;
//   else if( (abs(ginfo->id_1_a) == EGenType::kElectron && abs(ginfo->id_1_b) == EGenType::kMuon) || 
//         (abs(ginfo->id_1_a) == EGenType::kMuon && abs(ginfo->id_1_b) == EGenType::kElectron) ) gchannel=2; 
  
//   return gchannel;
// };
 
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-178078_7TeV_PromptReco_Collisions11_JSON.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"));  
  TFile * puf = new TFile("data/PileupReweighting.Summer11DYmm_To_Run2011A.root");
  hpu = (TH1D*)(puf->Get("puWeights"));
};

void initRunLumiRangeMap(RunLumiRangeMap &rlrm) {
  cout << "adding JSONS ... "  << endl;
  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-178078_7TeV_PromptReco_Collisions11_JSON.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
//                             TTree * 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::TGenInfo *ginfo ,     // input gen info
                               mithep::TEventInfo *info,     // input event info
                               TClonesArray *electronArr,    // input electrons
                               TClonesArray *muonArr,        // input muons
                               double eventweight,           // weight
                               TTree * passtuple ) { 

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

};

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


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

  unsigned evtfail = 0x0;
  unsigned gchannel=0xdeaddead;

//   if( ctrl.mc && ginfo != NULL ) { 
//     gchannel = getGenChannel(ginfo);
//   }


  if( ctrl.debug ) {
    cout << "Run: " << info->runNum 
         << "\tEvt: " << info->evtNum
         << "\tLumi: " << info->lumiSec
         << endl;
  }

  unsigned npu; double npuw;
  if( !ctrl.mc ) {
    // not accounting for overlap atm
    RunLumiRangeMap::RunLumiPairType rl(info->runNum, info->lumiSec);      
    if( !(rlrm.HasRunLumi(rl)) )  {
      if( ctrl.debug ) cout << "\tfails JSON" << endl;
      evtfail |= (1<<EVTFAIL_JSON);
      return evtfail;
    }
  } else {
    npu = info->nPU;
    npuw = hpu->GetBinContent(hpu->FindBin(npu)); 
  }
  
  
  //********************************************************
  // 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;
    if( ctrl.muSele == "ksWW" ) 
      muonfail = passKSMuonSelection(mu);
    else 
      muonfail = passMuonSelectionZZ(mu); 
    if( ctrl.debug ) { 
      cout << "muon:: pt: " << mu->pt 
           << "\teta: " << mu->eta 
           << "\tmask: 0x" << hex << muonfail << dec 
           << endl;
    }

#ifdef Z2FO
    if ( isMuFO(mu) ) {
#else
    if ( !muonfail ) {
#endif
      SimpleLepton tmplep;

      float pt = mu->pt;
      if( ctrl.do_escale_up   ) pt*=(1.01); 
      if( ctrl.do_escale_down ) pt*=(0.99);

      tmplep.vec.SetPtEtaPhiM(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.isoPF03 = mu->pfIso03;
      tmplep.isoPF04 = mu->pfIso04;
      tmplep.ip3dSig = mu->ip3dSig;
      tmplep.is4l    = false;
      tmplep.isEB    = (fabs(mu->eta) < 1.479 ? 1 : 0 ); 
      tmplep.isTight = (muonfail > 0 ? 0 : 1 );
      tmplep.isLoose = (muonfail > 0 ? 0 : 1 );
      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]);

    if( !(isEleFO(ctrl,ele) ) ) continue;
    if( ele->pt < 7 ) continue; //move this to ID

    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].isLoose && lepvec[k].type == 13 && lepvec[k].vec.Vect().DrEtaPhi(tmplep) < 0.1 ) {
        if( ctrl.debug ) cout << "-----> isMuonOverlap! " << endl;
        isMuonOverlap = kTRUE; 
        break;
      }
    }

    unsigned FAIL=0, isEleTight=0, isEleLoose=0;

    unsigned  failsCIC=0;
    CICStruct ciccuts_tight, ciccuts_medium, ciccuts_loose;
    if(ctrl.eleSele=="cic") {
      if( ctrl.eleSeleScheme == "mediumloose" ) {
        ciccuts_medium = getCiCCuts("medium");
        unsigned failsCICMedium = failsCicSelection(ctrl, ele, ciccuts_medium, ctrl.kinematics);
        ciccuts_loose = getCiCCuts("loose"); 
        unsigned failsCICLoose = failsCicSelection(ctrl, ele, ciccuts_loose, ctrl.kinematics);
        failsCIC = ( failsCICLoose > 0 && failsCICMedium > 0 );
        if( !failsCICMedium ) isEleTight=1;
      }
      else {
        ciccuts_tight = getCiCCuts(ctrl.eleSeleScheme); 
        failsCIC = failsCicSelection(ctrl, ele, ciccuts_tight, ctrl.kinematics);
        if( !failsCIC ) isEleTight=1;
      }
      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") {
      if( ctrl.eleSeleScheme == "mediumloose" ) {
        unsigned failsBDTMedium = failsBDTSelection(ctrl,"medium",ele);
        unsigned failsBDTLoose  = failsBDTSelection(ctrl,"loose",ele);
        failsBDT = ( failsBDTLoose > 0 && failsBDTMedium > 0 );
        if( !failsBDTMedium ) isEleTight=1;
        if( !failsBDTLoose ) isEleLoose=1;
      } else {
        failsBDT = failsBDTSelection(ctrl,"tight",ele);
        if( !failsBDT ) isEleTight=1;
      }
      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;

      float pt = ele->pt;
      if( ctrl.do_escale_up   ) pt*=(1.02);
      if( ctrl.do_escale_down ) pt*=(0.98);

      tmplep.vec.SetPtEtaPhiM( 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.isoPF03 = ele->pfIso03;
      tmplep.isoPF04 = ele->pfIso04;
      tmplep.ip3dSig = ele->ip3dSig;
      tmplep.is4l    = false;
      tmplep.isTight = isEleTight;
      tmplep.isEB    = ele->isEB;
      lepvec.push_back(tmplep);
      if( ctrl.debug ) { cout << "\telectron passes ... " << endl; }
    }
  }
  
  sort( lepvec.begin(), lepvec.end(), SimpleLepton::lep_pt_sort );

  for( int i=0; i<lepvec.size(); i++ ) {
    //
    //  check for a matched btagged jet
    //
    lepvec[i].tche = -1;
    lepvec[i].tchp = -1;
    lepvec[i].csv = -1;
    lepvec[i].csvMva = -1;

    if( jetArr != NULL ) {
      for(int k=0; k<jetArr->GetEntries(); k++) {
        const mithep::TJet *jet = (mithep::TJet*)((*jetArr)[k]);      
        TVector3 jvec;
        jvec.SetPtEtaPhi(jet->pt, jet->eta, jet->phi);
        if( jvec.DrEtaPhi( lepvec[i].vec.Vect()) < 0.2 ) { 
          lepvec[i].tche = jet->tche;
          lepvec[i].tchp = jet->tchp;
          lepvec[i].csv = jet->csv;
          lepvec[i].csvMva = jet->csvMva;
        }
      }
    }
  }

  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) {
    if( !(lepvec[i].isLoose ) ) continue; 
    for(int j = i+1; j < lepvec.size(); ++j) {
      if( !(lepvec[j].isLoose ) ) continue; 
      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) {

#ifndef Z2FO
    if( ctrl.eleSeleScheme == "mediumloose" && 
        !(lepvec[i].isTight) ) continue;
#endif

    for(int j = i+1; j < lepvec.size(); ++j) {

#ifndef Z2FO
      if( ctrl.eleSeleScheme == "mediumloose" && 
          !(lepvec[j].isTight) ) continue;
#endif

      //            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;

        if( ctrl.isoScheme == "pf" ) { 

          for( int i=0; i<lepvec.size(); i++ ) { 

            if( !(lepvec[i].is4l) ) continue;

            if( abs(lepvec[i].type) == 11 ) { 
              float reliso = lepvec[i].isoPF04/lepvec[i].vec.Pt();
              if( lepvec[i].isEB && lepvec[i].vec.Pt() > 20 && reliso > PFISO_ELE_LOOSE_EB_HIGHPT ) {
                failiso = true;
                break;
              }
              if( lepvec[i].isEB && lepvec[i].vec.Pt() < 20 && reliso > PFISO_ELE_LOOSE_EB_LOWPT ) {
                failiso = true;
                break;
              }
              if( !(lepvec[i].isEB) && lepvec[i].vec.Pt() > 20 && reliso > PFISO_ELE_LOOSE_EE_HIGHPT ) {
                failiso = true;
                break;
              }
              if( !(lepvec[i].isEB) && lepvec[i].vec.Pt() < 20 && reliso > PFISO_ELE_LOOSE_EE_LOWPT ) {
                failiso = true;
                break;
              }
            }

            if( abs(lepvec[i].type) == 13 ) { 

              /*
              float reliso = lepvec[i].isoPF03/lepvec[i].vec.Pt();
              if( lepvec[i].isEB && lepvec[i].vec.Pt() > 20 && reliso > PFISO_MU_LOOSE_EB_HIGHPT ) {  //0.13
                failiso = true;
                break;
              }
              if( lepvec[i].isEB && lepvec[i].vec.Pt() < 20 && reliso > PFISO_MU_LOOSE_EB_LOWPT ) { //0.06
                failiso = true;
                break;
              }
              if( !(lepvec[i].isEB) && lepvec[i].vec.Pt() > 20 && reliso > PFISO_MU_LOOSE_EE_HIGHPT ) { //0.09
                failiso = true;
                break;
              }
              if( !(lepvec[i].isEB) && lepvec[i].vec.Pt() < 20 && reliso > PFISO_MU_LOOSE_EE_LOWPT ) { //0.05
                failiso = true;
                break;
              }
              */
            }
          }  
        } else if ( ctrl.isoScheme == "pairwise" ) {
          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;
        }
        
        unsigned 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; 
        

        double w_offline=-1, werr_offline=0;
        double w_online=-1, werr_online=0;

        if( ctrl.mc ) { 

          vector< pair <double,double> > wlegs; // pair here is eff & err
          vector< pair <float,float> > mvec;  // pair here is eta & pt
          // now deal with medium vs loose
          //      vector< pair <float,float> > evec;  // pair here is eta & pt
          vector< pair< bool, pair <float,float> > > evec;  // pair here is eta & pt

          for( int k=0; k<lepvec.size(); k++ ) {
            if( !(lepvec[k].is4l) ) continue;
            if( abs(lepvec[k].type) == 13 ) {
              mvec.push_back( std::pair<float,float> (fabs(lepvec[k].vec.Eta()), lepvec[k].vec.Pt()) );
              wlegs.push_back( muonPerLegOfflineEfficiencyWeight( fabs(lepvec[k].vec.Eta()), 
                                                                  lepvec[k].vec.Pt() ) );
            } else { 

              // now deal with medium vs loose
              //              evec.push_back( std::pair<float,float> (fabs(lepvec[k].vec.Eta()), lepvec[k].vec.Pt()) );

              std::pair<float,float> tmppair(fabs(lepvec[k].vec.Eta()), lepvec[k].vec.Pt()); 
              evec.push_back( std::pair<bool, std::pair<float,float> > (lepvec[k].isTight, tmppair) );

              wlegs.push_back( elePerLegOfflineEfficiencyWeight(  fabs(lepvec[k].vec.Eta()), 
                                                                 lepvec[k].vec.Pt() ) );
            }
          }

          pair<double,double> offpair = getOfflineEfficiencyWeight( wlegs ); 
          w_offline    = offpair.first;
          werr_offline = offpair.second;

          pair<double,double> onpair  = getOnlineEfficiencyWeight( mvec, evec );
          w_online    = onpair.first;
          werr_online = onpair.second;

        } // if mc


        if( l != NULL ) {
          l->vec4l  = ZZSystem;
          l->vecz1  = Z1Candidate;
          l->vecz2  = Z2Candidate;
          l->vecl1p = lepvec[Z1LeptonPlusIndex].vec;
          l->vecl1m = lepvec[Z1LeptonMinusIndex].vec;
          l->vecl2p = lepvec[Z2LeptonPlusIndex].vec;
          l->vecl2m = lepvec[Z2LeptonMinusIndex].vec;
        }
        if( passtuple != NULL ) { 
          unsigned run   = info->runNum;
          unsigned evt   = info->evtNum;
          unsigned lumi  = info->lumiSec;
          unsigned chan  = channel;
          double   w     = eventweight;
          float mZ1      = Z1Candidate.M() ;
          float mZ2      = Z2Candidate.M() ;
          float m4l      = ZZSystem.M() ;
          float pt4l     = ZZSystem.Pt() ;
          unsigned tZ1   = abs(lepvec[Z1LeptonMinusIndex].type);
          unsigned tZ2   = abs(lepvec[Z2LeptonMinusIndex].type);
          float l3tche   = lepvec[Z2LeptonMinusIndex].tche;
          float l4tche   = lepvec[Z2LeptonPlusIndex].tche;
          float l3tchp   = lepvec[Z2LeptonMinusIndex].tchp;
          float l4tchp   = lepvec[Z2LeptonPlusIndex].tchp;
          float l3csv    = lepvec[Z2LeptonMinusIndex].csv;
          float l4csv    = lepvec[Z2LeptonPlusIndex].csv;
          float l3csvMva = lepvec[Z2LeptonMinusIndex].csvMva;
          float l4csvMva = lepvec[Z2LeptonPlusIndex].csvMva;


          passtuple->SetBranchAddress("channel",  &channel);
          passtuple->SetBranchAddress("run",  &run);
          passtuple->SetBranchAddress("evt",  &evt);
          passtuple->SetBranchAddress("lumi", &lumi);
          passtuple->SetBranchAddress("mZ1",  &mZ1);
          passtuple->SetBranchAddress("mZ2",  &mZ2);
          passtuple->SetBranchAddress("tZ1",  &tZ1);
          passtuple->SetBranchAddress("tZ2",  &tZ2);
          passtuple->SetBranchAddress("m4l",  &m4l);
          passtuple->SetBranchAddress("pt4l", &pt4l);
          passtuple->SetBranchAddress("w",    &w);
          if( ctrl.mc ) {
            passtuple->SetBranchAddress("woff", &w_offline );
            passtuple->SetBranchAddress("werroff", &werr_offline );
            passtuple->SetBranchAddress("won", &w_online );
            passtuple->SetBranchAddress("werron", &werr_online );
            passtuple->SetBranchAddress("npu", &npu );
            passtuple->SetBranchAddress("npuw", &npuw );
            passtuple->SetBranchAddress("gchannel",  &gchannel);
          }
          if( ctrl.btag ) {
          passtuple->SetBranchAddress("l3tche",    &l3tche);
          passtuple->SetBranchAddress("l4tche",    &l4tche);
          passtuple->SetBranchAddress("l3tchp",    &l3tchp);
          passtuple->SetBranchAddress("l4tchp",    &l4tchp);
          passtuple->SetBranchAddress("l3csv",     &l3csv);
          passtuple->SetBranchAddress("l4csv",     &l4csv);
          passtuple->SetBranchAddress("l3csvMva",  &l3csvMva);
          passtuple->SetBranchAddress("l4csvMva",  &l4csvMva);
          }

          passtuple->Fill( ); 
        }

        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.17
Committed:	Mon Nov 7 06:53:01 2011 UTC (13 years, 6 months ago) by khahn
Content type:	text/plain
Branch:	MAIN
Changes since 1.16:	+20 -6 lines
Log Message:	ifdef for Z2FO