Selection/src/Selection.cc

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

mithep::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
                               TNtuple * passtuple ) { 

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

};

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

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

};


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

  unsigned evtfail = 0x0;
  double eventweight = info->eventweight;

  if( !ctrl.mc ) {
    // not accounting for overlap atm
    mithep::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;
  
  //    
  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 ) {
        cout << "-----> isMuonOverlap! " << endl;
        isMuonOverlap = kTRUE; 
        break;
      }
    }

    CICStruct tightcuts = getTightCuts(); 
    unsigned failsCIC = failsCicSelection(ele, tightcuts);
    if( ctrl.debug ){ 
      cout << "CIC category: " << cicCategory(ele) 
           << "\tfailmask : 0x" << hex << failsCIC  << dec
           << "\tscEt: " << ele->scEt 
           << "\tscEta: " << ele->scEta 
           << endl;
    }
    if ( !failsCIC && !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();
      cout << "Z1 selection, tmpZ1Mass: " << tmpZ1Mass << endl;
      if( tmpZ1Mass > 60 ) {
        if (fabs(tmpZ1Mass - 91.1876) < fabs(BestZ1Mass - 91.1876)) {
          BestZ1Mass = tmpZ1Mass;
          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; 
  }
  cout << "\tgot a Z1 ... run: " << info->runNum << "\tevt: " << info->evtNum  << endl;
  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;
  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;

            cout << "dilepton.M() : " << dilepton.M() << endl;
            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;
              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;
        }
        cout << "\tgot a Z2 ..." << endl;
        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->PileupEnergyDensity;
        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 ) { 
                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 
        //***************************************************************
        if ( Z1Candidate.M() > 120 || 
             Z2Candidate.M() < 20  ||
             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);
        }

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