Selection/src/SelectionZ4L.cc

#include "SelectionStatus.h" 
#include "EventData.h" 
#include "SimpleLepton.h" 
#include "EfficiencyWeightsInterface.h"

#include "HZZCiCElectronSelection.h"
//#include "IsolationSelection.h"
#include "PassHLT.h"
#include "SelectionZ4L.h"

#include "ExternData.h"
#include "SelectionDefs.h"


// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
EventData apply_Z4L_selection(ControlFlags &ctrl,           // input control 
                                mithep::TEventInfo *info,     // input event info
                                TClonesArray *electronArr,    // input electrons
                                CICStruct &cicCuts,
                                SelectionStatus (*ElectronPreSelector)( ControlFlags &, 
                                                                        const mithep::TElectron*),
                                SelectionStatus (*ElectronIDSelector)( ControlFlags &, 
                                                                       const mithep::TElectron*, 
                                                                       CICStruct &cicCuts),
                                TClonesArray *muonArr,         // input muons
                                SelectionStatus (*MuonPreSelector)( ControlFlags &, const mithep::TMuon*),
                                SelectionStatus (*MuonIDSelector)( ControlFlags &, const mithep::TMuon*),
                                // additional post-ID selector 
                                bool (*PairwiseIsoSelector)( ControlFlags &, 
                                                             vector<SimpleLepton> &,
                                                             float ))
  // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
{       

  EventData ret;
  unsigned evtfail = 0x0;
  TRandom3 r;
        
  if( ctrl.debug ) {
    cout << "Run: " << info->runNum 
         << "\tEvt: " << info->evtNum
         << "\tLumi: " << info->lumiSec
         << endl;
  }

  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;
      ret.status.setStatus(0);
      return ret;
    }
  }
  
  
  //********************************************************
  // Trigger
  //********************************************************
  if( !ctrl.mc ) { 
    //if( !(passHLT(info->triggerBits, info->runNum, channel) )  ) { 
    if( !(passHLT(info->triggerBits, info->runNum, 999) )  ) { 
      if( ctrl.debug ) cout << "\tfails trigger" << endl;
      evtfail |= (1<<EVTFAIL_TRIGGER);
      ret.status.setStatus(0);
      return ret;
    }   
  } 
  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]);      
      
      SelectionStatus musel;
      if(ctrl.debug) cout << "musel.status  before anything: " << musel.getStatus() << endl;
      musel |= (*MuonPreSelector)(ctrl,mu);
      if(ctrl.debug) cout << "musel.status  after presel: " << musel.getStatus() << endl;
      musel |= (*MuonIDSelector)(ctrl,mu);
      if(ctrl.debug) cout << "musel.status  after ID: " << musel.getStatus() << endl;


      if( ctrl.debug ) { 
        cout << "muon:: pt: " << mu->pt 
             << "\teta: " << mu->eta 
             << "\tisorel: " <<  mu->pfIso03/mu->pt
             << "\tstatus: " << hex << musel.getStatus() << dec
             << endl;
      }
      
      if ( musel.passID() ) {
        SimpleLepton tmplep;
        float pt = mu->pt;
        tmplep.vecorig->SetPtEtaPhiM(pt, 
                                    mu->eta, 
                                    mu->phi,
                                    MUON_MASS);
        
        if( ctrl.do_escale_up ) { 
          pt=scale_smear_muon_Up(pt, 1,  r);
        }
        if( ctrl.do_escale_down ) { 
          pt=scale_smear_muon_Down(pt, 1,  r);
        }
        
        tmplep.vec->SetPtEtaPhiM(pt, 
                                mu->eta, 
                                mu->phi,
                                MUON_MASS);
        
        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 = musel.tightID();
        tmplep.isLoose = musel.looseID();
        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].isLoose && lepvec[k].type == 13 && lepvec[k].vec->Vect().DrEtaPhi(tmplep) < 0.05 ) {
            if( ctrl.debug ) cout << "-----> isMuonOverlap! " << endl;
            isMuonOverlap = kTRUE; 
            break;
          }
        }
        
        SelectionStatus elesel;
        if( ctrl.debug ) cout << "--> status before anything: " << hex << elesel.getStatus() << dec << endl;
        elesel |= (*ElectronPreSelector)(ctrl,ele);
        if( ctrl.debug ) cout << "--> status after presel: " << hex << elesel.getStatus() << dec << endl;
        elesel |= (*ElectronIDSelector)(ctrl,ele,cicCuts);
        if( ctrl.debug ) cout << "--> status after ID: " << hex << elesel.getStatus() << dec << endl;

        if( ctrl.debug ){ 
          cout << "\tscEt: " << ele->scEt 
               << "\tscEta: " << ele->scEta
               << "\tncluster: " << ele->ncluster
               << "\tisorel: " <<  ele->pfIso04/ele->pt
               << "\tstatus: " << hex << elesel.getStatus() << dec
               << endl;
        }

        if ( elesel.passIDAndPre() && !isMuonOverlap ) 
          {
            SimpleLepton tmplep;
           
                        
            float pt = ele->pt;
                        
            tmplep.vecorig->SetPtEtaPhiM( pt,
                                          ele->eta,
                                          ele->phi,
                                          ELECTRON_MASS );
            
            if( ctrl.do_escale ) { 
              pt=scale_smear_electron(pt, ele->isEB, r);
            }
            if( ctrl.do_escale_up ) { 
              pt=scale_smear_electron_Up(pt, ele->isEB,  r);
            }
            if( ctrl.do_escale_down ) { 
              pt=scale_smear_electron_Down(pt, ele->isEB,  r);
            }

            
            tmplep.vec->SetPtEtaPhiM( pt,
                                      ele->eta,
                                      ele->phi,
                                      ELECTRON_MASS );
            
            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.isEB    = ele->isEB;
            tmplep.scID    = ele->scID;
            tmplep.isTight = elesel.tightID();
            tmplep.isLoose = elesel.looseID();
            lepvec.push_back(tmplep);
            if( ctrl.debug ) { cout << "\telectron passes ... " << endl; }
          }
      }
    
    //********************************************************
    // Dump Stuff
    //********************************************************
    sort( lepvec.begin(), lepvec.end(), SimpleLepton::lep_pt_sort );
    int nmu=0, nele=0;
    for( int i=0; i<lepvec.size(); i++ ) {
      if(ctrl.debug) cout << "lepvec :: index: " << i 
                          << "\tpt: " << lepvec[i].vec->Pt()
                          << "\ttype: " << lepvec[i].type 
                          << endl;
      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].isTight) ) continue; 
      for(int j = i+1; j < lepvec.size(); ++j) {
        if( !(lepvec[j].isTight) ) 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 - Z_MASS) < fabs(BestZ1Mass - Z_MASS)) {
            BestZ1Mass = tmpZ1Mass;
            if( ctrl.debug ) cout << "Z1 selection, new BestZ1Mass: " << BestZ1Mass 
                                  << "\tdM: " << fabs(BestZ1Mass - Z_MASS)
                                  << 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);
      //ret.status = evtfail;
      ret.status.setStatus(0);
      return ret;
    }
    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;
    
    //******************************************************************************
    // Z1 + l
    //******************************************************************************
    if( lepvec.size() < 3 ) { 
      evtfail |= (1<<EVTFAIL_Z1_PLUSL);
      //ret.status = evtfail;
      ret.status.setStatus(0);
      return ret;
    }
    
    //******************************************************************************
    // 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) {
      
      if( ctrl.debug)  cout << "i: " << i 
                            << "\tpt: " << lepvec[i].vec->Pt() 
                            << "\ttype: " << lepvec[i].type 
                            << endl;
      
      if( !(lepvec[i].isTight) ) { 
        if( ctrl.debug)  cout << "it's not tight, skipping ... " << endl;
        continue;
      }
      
      for(int j = i+1; j < lepvec.size(); ++j) {
        if( ctrl.debug)  cout << "\t\tj: " << j 
                              << "\tpt: " << lepvec[j].vec->Pt() 
                              << "\ttype: " << lepvec[j].type 
                              << endl;
        
        if( !(lepvec[j].isTight) ) { 
          if( ctrl.debug)  cout << "it's not tight, skipping ... " << endl;
          continue;
        }
        
        
        if (i == Z1LeptonPlusIndex || i == Z1LeptonMinusIndex) { 
          if( ctrl.debug)  cout << "\ti matches a Z1 index, skipping ..." << endl;
          continue; //skip Z1 leptons
        }
        if (j == Z1LeptonPlusIndex || j == Z1LeptonMinusIndex) { 
          if( ctrl.debug)  cout << "\tj matches a Z1 index, skipping ..." << endl;
          continue; //skip Z1 leptons
        }
        if (lepvec[i].charge == lepvec[j].charge) { 
          if( ctrl.debug)  cout << "\ti and j are same sign, skipping ..." << endl;
          continue;         //require opp sign
        }
        if (fabs(lepvec[i].type) != fabs(lepvec[j].type)) { 
          if( ctrl.debug) 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() > 2.0)) continue;
        if (!(fourLepton.M() > 60.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() > 2 && pair2.M() > 2)) 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 );
      //      ret.status = evtfail;
      ret.status.setStatus(0);
      return ret;
    }
    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;
    lepvec[Z1LeptonPlusIndex].is4l = true;
    lepvec[Z1LeptonMinusIndex].is4l = true;
    lepvec[Z2LeptonPlusIndex].is4l = true;
    lepvec[Z2LeptonMinusIndex].is4l = true;

                
    //***************************************************************
    // remaining cuts :  kinematic & iso 
    //***************************************************************

    if( !(PairwiseIsoSelector(ctrl, lepvec, info->rho))  ) { 
      ret.status.setStatus(0);
      return ret;
    }

    double Z2massCut= 2;
  
   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 );
      //      ret.status = evtfail;
      ret.status.setStatus(0);
      return ret;
    }
  
    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; 
  
  
    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;
  

    //***************************************************************
    // finish
    //***************************************************************

    if( !evtfail ) { 
      ret.status.setStatus(SelectionStatus::EVTPASS);
      ret.Z1leptons.push_back(lepvec[Z1LeptonMinusIndex]);
      ret.Z1leptons.push_back(lepvec[Z1LeptonPlusIndex]);
      ret.Z2leptons.push_back(lepvec[Z2LeptonMinusIndex]);
      ret.Z2leptons.push_back(lepvec[Z2LeptonPlusIndex]);
    }
                
                
    return ret;
}


Revision:	1.4
Committed:	Wed Feb 29 09:43:13 2012 UTC (13 years, 2 months ago) by anlevin
Content type:	text/plain
Branch:	MAIN
CVS Tags:	synced_FSR_2, synced_FSR, synched2, synched
Changes since 1.3:	+21 -23 lines
Log Message:	changed Z4l selection code
#	User	Rev	Content
1	khahn	1.1	#include "SelectionStatus.h"
2			#include "EventData.h"
3			#include "SimpleLepton.h"
4			#include "EfficiencyWeightsInterface.h"
5
6			#include "HZZCiCElectronSelection.h"
7	anlevin	1.3	//#include "IsolationSelection.h"
8	khahn	1.1	#include "PassHLT.h"
9			#include "SelectionZ4L.h"
10
11			#include "ExternData.h"
12			#include "SelectionDefs.h"
13
14
15			// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
16			// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
17			// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
18	anlevin	1.2	EventData apply_Z4L_selection(ControlFlags &ctrl, // input control
19	khahn	1.1	mithep::TEventInfo *info, // input event info
20			TClonesArray *electronArr, // input electrons
21			CICStruct &cicCuts,
22			SelectionStatus (*ElectronPreSelector)( ControlFlags &,
23			const mithep::TElectron*),
24			SelectionStatus (*ElectronIDSelector)( ControlFlags &,
25			const mithep::TElectron*,
26			CICStruct &cicCuts),
27			TClonesArray *muonArr, // input muons
28			SelectionStatus (MuonPreSelector)( ControlFlags &, const mithep::TMuon),
29			SelectionStatus (MuonIDSelector)( ControlFlags &, const mithep::TMuon),
30			// additional post-ID selector
31			bool (*PairwiseIsoSelector)( ControlFlags &,
32			vector<SimpleLepton> &,
33			float ))
34			// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
35			{
36
37			EventData ret;
38			unsigned evtfail = 0x0;
39			TRandom3 r;
40	anlevin	1.4
41	khahn	1.1	if( ctrl.debug ) {
42			cout << "Run: " << info->runNum
43			<< "\tEvt: " << info->evtNum
44			<< "\tLumi: " << info->lumiSec
45			<< endl;
46			}
47
48			if( !ctrl.mc ) {
49			// not accounting for overlap atm
50			RunLumiRangeMap::RunLumiPairType rl(info->runNum, info->lumiSec);
51			if( !(rlrm.HasRunLumi(rl)) ) {
52			if( ctrl.debug ) cout << "\tfails JSON" << endl;
53			ret.status.setStatus(0);
54			return ret;
55			}
56			}
57
58
59			//********************************************************
60			// Trigger
61			//********************************************************
62			if( !ctrl.mc ) {
63			//if( !(passHLT(info->triggerBits, info->runNum, channel) ) ) {
64			if( !(passHLT(info->triggerBits, info->runNum, 999) ) ) {
65			if( ctrl.debug ) cout << "\tfails trigger" << endl;
66			evtfail \|= (1<<EVTFAIL_TRIGGER);
67			ret.status.setStatus(0);
68			return ret;
69			}
70			}
71			if( ctrl.debug ) {
72			cout << "presel nlep: " << muonArr->GetEntries() + electronArr->GetEntries()
73			<< "\tnmuon: " << muonArr->GetEntries()
74			<< "\tnelectron: " << electronArr->GetEntries()
75			<< endl;
76			}
77
78			//********************************************************
79			// Lepton Selection
80			//********************************************************
81			vector<SimpleLepton> lepvec;
82	anlevin	1.4
83	khahn	1.1	//
84			if( ctrl.debug ) cout << "\tnMuons: " << muonArr->GetEntries() << endl;
85			//----------------------------------------------------
86			for(Int_t i=0; i<muonArr->GetEntries(); i++)
87			{
88			const mithep::TMuon mu = (mithep::TMuon)((*muonArr)[i]);
89
90			SelectionStatus musel;
91			if(ctrl.debug) cout << "musel.status before anything: " << musel.getStatus() << endl;
92			musel \|= (*MuonPreSelector)(ctrl,mu);
93			if(ctrl.debug) cout << "musel.status after presel: " << musel.getStatus() << endl;
94			musel \|= (*MuonIDSelector)(ctrl,mu);
95			if(ctrl.debug) cout << "musel.status after ID: " << musel.getStatus() << endl;
96
97
98			if( ctrl.debug ) {
99			cout << "muon:: pt: " << mu->pt
100			<< "\teta: " << mu->eta
101			<< "\tisorel: " << mu->pfIso03/mu->pt
102			<< "\tstatus: " << hex << musel.getStatus() << dec
103			<< endl;
104			}
105
106			if ( musel.passID() ) {
107			SimpleLepton tmplep;
108			float pt = mu->pt;
109			tmplep.vecorig->SetPtEtaPhiM(pt,
110			mu->eta,
111			mu->phi,
112			MUON_MASS);
113
114			if( ctrl.do_escale_up ) {
115			pt=scale_smear_muon_Up(pt, 1, r);
116			}
117			if( ctrl.do_escale_down ) {
118			pt=scale_smear_muon_Down(pt, 1, r);
119			}
120
121			tmplep.vec->SetPtEtaPhiM(pt,
122			mu->eta,
123			mu->phi,
124			MUON_MASS);
125
126			tmplep.type = 13;
127			tmplep.index = i;
128			tmplep.charge = mu->q;
129			tmplep.isoTrk = mu->trkIso03;
130			tmplep.isoEcal = mu->emIso03;
131			tmplep.isoHcal = mu->hadIso03;
132			tmplep.isoPF03 = mu->pfIso03;
133			tmplep.isoPF04 = mu->pfIso04;
134			tmplep.ip3dSig = mu->ip3dSig;
135			tmplep.is4l = false;
136			tmplep.isEB = (fabs(mu->eta) < 1.479 ? 1 : 0 );
137	anlevin	1.4	tmplep.isTight = musel.tightID();
138			tmplep.isLoose = musel.looseID();
139	khahn	1.1	lepvec.push_back(tmplep);
140			if( ctrl.debug ) { cout << "muon passes ... " << endl;}
141			}
142			// }
143			}
144
145
146
147			//
148			if( ctrl.debug ) { cout << "\tnElectron: " << electronArr->GetEntries() << endl; }
149			// --------------------------------------------------------------------------------
150			for(Int_t i=0; i<electronArr->GetEntries(); i++)
151			{
152			const mithep::TElectron ele = (mithep::TElectron)((*electronArr)[i]);
153
154			Bool_t isMuonOverlap = kFALSE;
155			for (int k=0; k<lepvec.size(); ++k) {
156			TVector3 tmplep;
157			tmplep.SetPtEtaPhi(ele->pt, ele->eta, ele->phi);
158	anlevin	1.4	if ( lepvec[k].isLoose && lepvec[k].type == 13 && lepvec[k].vec->Vect().DrEtaPhi(tmplep) < 0.05 ) {
159	khahn	1.1	if( ctrl.debug ) cout << "-----> isMuonOverlap! " << endl;
160			isMuonOverlap = kTRUE;
161			break;
162			}
163			}
164
165			SelectionStatus elesel;
166			if( ctrl.debug ) cout << "--> status before anything: " << hex << elesel.getStatus() << dec << endl;
167			elesel \|= (*ElectronPreSelector)(ctrl,ele);
168			if( ctrl.debug ) cout << "--> status after presel: " << hex << elesel.getStatus() << dec << endl;
169			elesel \|= (*ElectronIDSelector)(ctrl,ele,cicCuts);
170			if( ctrl.debug ) cout << "--> status after ID: " << hex << elesel.getStatus() << dec << endl;
171
172			if( ctrl.debug ){
173			cout << "\tscEt: " << ele->scEt
174			<< "\tscEta: " << ele->scEta
175			<< "\tncluster: " << ele->ncluster
176			<< "\tisorel: " << ele->pfIso04/ele->pt
177			<< "\tstatus: " << hex << elesel.getStatus() << dec
178			<< endl;
179			}
180
181	anlevin	1.4	if ( elesel.passIDAndPre() && !isMuonOverlap )
182	khahn	1.1	{
183			SimpleLepton tmplep;
184	anlevin	1.4
185
186	khahn	1.1	float pt = ele->pt;
187	anlevin	1.4
188	khahn	1.1	tmplep.vecorig->SetPtEtaPhiM( pt,
189			ele->eta,
190			ele->phi,
191			ELECTRON_MASS );
192
193			if( ctrl.do_escale ) {
194			pt=scale_smear_electron(pt, ele->isEB, r);
195			}
196			if( ctrl.do_escale_up ) {
197			pt=scale_smear_electron_Up(pt, ele->isEB, r);
198			}
199			if( ctrl.do_escale_down ) {
200			pt=scale_smear_electron_Down(pt, ele->isEB, r);
201			}
202
203
204			tmplep.vec->SetPtEtaPhiM( pt,
205			ele->eta,
206			ele->phi,
207			ELECTRON_MASS );
208
209			tmplep.type = 11;
210			tmplep.index = i;
211			tmplep.charge = ele->q;
212			tmplep.isoTrk = ele->trkIso03;
213			tmplep.isoEcal = ele->emIso03;
214			tmplep.isoHcal = ele->hadIso03;
215			tmplep.isoPF03 = ele->pfIso03;
216			tmplep.isoPF04 = ele->pfIso04;
217			tmplep.ip3dSig = ele->ip3dSig;
218			tmplep.is4l = false;
219			tmplep.isEB = ele->isEB;
220			tmplep.scID = ele->scID;
221			tmplep.isTight = elesel.tightID();
222			tmplep.isLoose = elesel.looseID();
223			lepvec.push_back(tmplep);
224			if( ctrl.debug ) { cout << "\telectron passes ... " << endl; }
225			}
226			}
227
228			//********************************************************
229			// Dump Stuff
230			//********************************************************
231			sort( lepvec.begin(), lepvec.end(), SimpleLepton::lep_pt_sort );
232			int nmu=0, nele=0;
233			for( int i=0; i<lepvec.size(); i++ ) {
234			if(ctrl.debug) cout << "lepvec :: index: " << i
235			<< "\tpt: " << lepvec[i].vec->Pt()
236			<< "\ttype: " << lepvec[i].type
237			<< endl;
238			if( abs(lepvec[i].type) == 11 ) nele++;
239			else nmu++;
240			}
241			if( ctrl.debug ) {
242			cout << "postsel nlep: " << lepvec.size()
243			<< "\tnmuon: " << nmu
244			<< "\tnelectron: " << nele
245			<< endl;
246	anlevin	1.4	}
247
248	khahn	1.1	//******************************************************************************
249			// Z1 Selection
250			//******************************************************************************
251			int Z1LeptonPlusIndex = -1;
252			int Z1LeptonMinusIndex = -1;
253			double BestZ1Mass = -999;
254			if( ctrl.debug ) { cout << "looking for a Z1 ..." << endl; }
255			for(int i = 0; i < lepvec.size(); ++i) {
256			if( !(lepvec[i].isTight) ) continue;
257			for(int j = i+1; j < lepvec.size(); ++j) {
258			if( !(lepvec[j].isTight) ) continue;
259			if( ctrl.debug ) { cout << "\tconsidering leptons " << i << " & " << j << endl; }
260			if (!(lepvec[i].vec->Pt() > 20.0 \|\| lepvec[j].vec->Pt() > 20.0)) continue;
261			if( ctrl.debug ) { cout << "\tat least one is > 20 GeV" << endl; }
262			if (!(lepvec[i].vec->Pt() > 10.0 && lepvec[j].vec->Pt() > 10.0)) continue;
263			if( ctrl.debug ) { cout << "\tthe other is > 10 GeV" << endl; }
264			if (lepvec[i].charge == lepvec[j].charge) continue;
265			if( ctrl.debug ) { cout << "\tthey're opposite charge" << endl; }
266			if (fabs(lepvec[i].type) != fabs(lepvec[j].type)) continue;
267			if( ctrl.debug ) { cout << "\tthey're same flavor" << endl; }
268
269			//Make Z1 hypothesis
270			TLorentzVector leptonPlus, leptonMinus;
271			if ( lepvec[i].charge > 0 ) {
272			leptonPlus = lepvec[i].vec;
273			leptonMinus = lepvec[j].vec;
274			} else {
275			leptonPlus = lepvec[j].vec;
276			leptonMinus = lepvec[i].vec;
277			}
278
279			float tmpZ1Mass = (leptonPlus + leptonMinus).M();
280			if( ctrl.debug ) cout << "Z1 selection, tmpZ1Mass: " << tmpZ1Mass << endl;
281			if( tmpZ1Mass > 60 ) {
282			if (fabs(tmpZ1Mass - Z_MASS) < fabs(BestZ1Mass - Z_MASS)) {
283			BestZ1Mass = tmpZ1Mass;
284			if( ctrl.debug ) cout << "Z1 selection, new BestZ1Mass: " << BestZ1Mass
285			<< "\tdM: " << fabs(BestZ1Mass - Z_MASS)
286			<< endl;
287			if (lepvec[i].charge > 0) {
288			Z1LeptonPlusIndex = i;
289			Z1LeptonMinusIndex = j;
290			} else {
291			Z1LeptonPlusIndex = j;
292			Z1LeptonMinusIndex = i;
293			}
294			}
295			}
296			}
297			}
298			// stop if no Z1 candidate is found
299			if( BestZ1Mass < 0 ) {
300			evtfail \|= (1<<EVTFAIL_Z1);
301			//ret.status = evtfail;
302			ret.status.setStatus(0);
303			return ret;
304			}
305			if( ctrl.debug ) cout << "\tgot a Z1 ... run: " << info->runNum << "\tevt: " << info->evtNum << endl;
306			if( ctrl.debug ) cout << "\tZ1 plusindex: " << Z1LeptonPlusIndex << "\tminusindex: " << Z1LeptonMinusIndex << endl;
307			TLorentzVector Z1LeptonPlus = *(lepvec[Z1LeptonPlusIndex].vec);
308			TLorentzVector Z1LeptonMinus = *(lepvec[Z1LeptonMinusIndex].vec);
309			TLorentzVector Z1Candidate = Z1LeptonPlus + Z1LeptonMinus;
310
311			//******************************************************************************
312			// Z1 + l
313			//******************************************************************************
314			if( lepvec.size() < 3 ) {
315			evtfail \|= (1<<EVTFAIL_Z1_PLUSL);
316			//ret.status = evtfail;
317			ret.status.setStatus(0);
318			return ret;
319			}
320
321			//******************************************************************************
322			// 4l/Z2 Selection
323			//******************************************************************************
324			Int_t Z2LeptonPlusIndex = -1;
325			Int_t Z2LeptonMinusIndex = -1;
326			Double_t BestZ2Mass = -1;
327			if( ctrl.debug ) cout << "looking for a Z2 ... out of " << lepvec.size() << " leptons" <<endl;
328			for(int i = 0; i < lepvec.size(); ++i) {
329
330			if( ctrl.debug) cout << "i: " << i
331			<< "\tpt: " << lepvec[i].vec->Pt()
332			<< "\ttype: " << lepvec[i].type
333			<< endl;
334
335	anlevin	1.3	if( !(lepvec[i].isTight) ) {
336	khahn	1.1	if( ctrl.debug) cout << "it's not tight, skipping ... " << endl;
337			continue;
338			}
339
340			for(int j = i+1; j < lepvec.size(); ++j) {
341			if( ctrl.debug) cout << "\t\tj: " << j
342			<< "\tpt: " << lepvec[j].vec->Pt()
343			<< "\ttype: " << lepvec[j].type
344			<< endl;
345
346	anlevin	1.3	if( !(lepvec[j].isTight) ) {
347	khahn	1.1	if( ctrl.debug) cout << "it's not tight, skipping ... " << endl;
348			continue;
349			}
350
351
352			if (i == Z1LeptonPlusIndex \|\| i == Z1LeptonMinusIndex) {
353			if( ctrl.debug) cout << "\ti matches a Z1 index, skipping ..." << endl;
354			continue; //skip Z1 leptons
355			}
356			if (j == Z1LeptonPlusIndex \|\| j == Z1LeptonMinusIndex) {
357			if( ctrl.debug) cout << "\tj matches a Z1 index, skipping ..." << endl;
358			continue; //skip Z1 leptons
359			}
360			if (lepvec[i].charge == lepvec[j].charge) {
361			if( ctrl.debug) cout << "\ti and j are same sign, skipping ..." << endl;
362			continue; //require opp sign
363			}
364			if (fabs(lepvec[i].type) != fabs(lepvec[j].type)) {
365			if( ctrl.debug) cout << "\ti and j are not same flavor, skipping ..." << endl;
366			continue; //require same flavor
367			}
368
369			//Make Z2 hypothesis
370			TLorentzVector leptonPlus, leptonMinus;
371
372			if (lepvec[i].charge > 0 ) {
373			leptonPlus = lepvec[i].vec;
374			leptonMinus = lepvec[j].vec;
375			} else {
376			leptonPlus = lepvec[j].vec;
377			leptonMinus = lepvec[i].vec;
378			}
379
380			TLorentzVector dilepton = leptonPlus + leptonMinus;
381			TLorentzVector fourLepton = Z1Candidate + dilepton;
382
383			if( ctrl.debug ) cout << "dilepton.M() : " << dilepton.M() << endl;
384			if( ctrl.debug ) cout << "fourLepton.M() : " << fourLepton.M() << endl;
385
386	anlevin	1.4	if (!(dilepton.M() > 2.0)) continue;
387			if (!(fourLepton.M() > 60.0)) continue;
388	khahn	1.1
389			//for 4e and 4mu, require at least 1 of the other opp sign lepton pairs have mass > 12
390			if (fabs(lepvec[i].type) == fabs(lepvec[Z1LeptonPlusIndex].type)) {
391			TLorentzVector pair1 = Z1LeptonPlus + *leptonMinus;
392			TLorentzVector pair2 = Z1LeptonMinus + *leptonPlus;
393			if( ctrl.debug ) cout << "pair1: " << pair1.M() << "\tpair2: "<< pair2.M() << endl;
394	anlevin	1.4	if (!(pair1.M() > 2 && pair2.M() > 2)) continue;
395	khahn	1.1	}
396	anlevin	1.4
397
398	khahn	1.1
399			//Disambiguiation is done by choosing the pair with the largest ptMax and largest ptMin
400			if (Z2LeptonPlusIndex < 0) {
401			if (lepvec[i].charge > 0) {
402			Z2LeptonPlusIndex = i;
403			Z2LeptonMinusIndex = j;
404			} else {
405			Z2LeptonPlusIndex = j;
406			Z2LeptonMinusIndex = i;
407			}
408			} else {
409			Double_t BestPairPtMax = lepvec[Z2LeptonPlusIndex].vec->Pt();
410			Double_t BestPairPtMin = lepvec[Z2LeptonMinusIndex].vec->Pt();
411			if (lepvec[Z2LeptonMinusIndex].vec->Pt() > BestPairPtMax) {
412			BestPairPtMax = lepvec[Z2LeptonMinusIndex].vec->Pt();
413			BestPairPtMin = lepvec[Z2LeptonPlusIndex].vec->Pt();
414			}
415
416			Double_t CurrentPairPtMax = lepvec[i].vec->Pt();
417			Double_t CurrentPairPtMin = lepvec[j].vec->Pt();
418			if (lepvec[j].vec->Pt() > CurrentPairPtMax) {
419			CurrentPairPtMax = lepvec[j].vec->Pt();
420			CurrentPairPtMin = lepvec[i].vec->Pt();
421			}
422
423			if (CurrentPairPtMax > BestPairPtMax) {
424			if (lepvec[i].charge > 0) {
425			Z2LeptonPlusIndex = i;
426			Z2LeptonMinusIndex = j;
427			} else {
428			Z2LeptonPlusIndex = j;
429			Z2LeptonMinusIndex = i;
430			}
431			} else if (CurrentPairPtMax == BestPairPtMax) {
432			if (CurrentPairPtMin > BestPairPtMin) {
433			if (lepvec[i].charge > 0) {
434			Z2LeptonPlusIndex = i;
435			Z2LeptonMinusIndex = j;
436			} else {
437			Z2LeptonPlusIndex = j;
438			Z2LeptonMinusIndex = i;
439			}
440			}
441			}
442			}
443			}
444			}
445
446			// stop if no Z2 candidate is found
447			if (Z2LeptonPlusIndex == -1) {
448			evtfail \|= ( 1<<EVTFAIL_4L );
449			// ret.status = evtfail;
450			ret.status.setStatus(0);
451			return ret;
452			}
453			if( ctrl.debug ) cout << "\tgot a Z2 ..." << endl;
454			if( ctrl.debug ) cout << "\tZ2 plusindex: " << Z2LeptonPlusIndex
455			<< "\tminusindex: " << Z2LeptonMinusIndex << endl;
456			TLorentzVector Z2LeptonPlus = *(lepvec[Z2LeptonPlusIndex].vec);
457			TLorentzVector Z2LeptonMinus = *(lepvec[Z2LeptonMinusIndex].vec);
458			TLorentzVector Z2Candidate = Z2LeptonPlus+Z2LeptonMinus;
459			TLorentzVector ZZSystem = Z1Candidate + Z2Candidate;
460			lepvec[Z1LeptonPlusIndex].is4l = true;
461			lepvec[Z1LeptonMinusIndex].is4l = true;
462			lepvec[Z2LeptonPlusIndex].is4l = true;
463			lepvec[Z2LeptonMinusIndex].is4l = true;
464
465	anlevin	1.4
466	khahn	1.1
467			//***************************************************************
468			// remaining cuts : kinematic & iso
469			//***************************************************************
470
471			if( !(PairwiseIsoSelector(ctrl, lepvec, info->rho)) ) {
472			ret.status.setStatus(0);
473			return ret;
474			}
475
476	anlevin	1.4	double Z2massCut= 2;
477	khahn	1.1
478	anlevin	1.4	if ( Z1Candidate.M() > 120 \|\|
479	khahn	1.1	Z2Candidate.M() < Z2massCut \|\|
480			Z2Candidate.M() > 120 \|\|
481			!(lepvec[Z1LeptonPlusIndex].vec->Pt() > 20.0 \|\| lepvec[Z1LeptonMinusIndex].vec->Pt() > 20.0) \|\|
482			!(lepvec[Z1LeptonPlusIndex].vec->Pt() > 10.0 && lepvec[Z1LeptonMinusIndex].vec->Pt() > 10.0)
483			) {
484			evtfail \|= (1<<EVTFAIL_KINEMATICS );
485			// ret.status = evtfail;
486			ret.status.setStatus(0);
487			return ret;
488			}
489
490			unsigned channel;
491			if( lepvec[Z1LeptonMinusIndex].type == 11 && lepvec[Z2LeptonMinusIndex].type == 11 ) channel=0;
492			if( lepvec[Z1LeptonMinusIndex].type == 13 && lepvec[Z2LeptonMinusIndex].type == 13 ) channel=1;
493			if( (lepvec[Z1LeptonMinusIndex].type == 11 && lepvec[Z2LeptonMinusIndex].type == 13) \|\|
494			(lepvec[Z1LeptonMinusIndex].type == 13 && lepvec[Z2LeptonMinusIndex].type == 11)) channel=2;
495
496
497
498
499			if( ctrl.debug ) cout << "run: " << info->runNum
500			<< "\tevt: " << info->evtNum
501			<< "\tZ1channel: " << lepvec[Z1LeptonMinusIndex].type
502			<< "\tZ2channel: " << lepvec[Z2LeptonMinusIndex].type
503			<< "\tmZ1: " << Z1Candidate.M()
504			<< "\tmZ2: " << Z2Candidate.M()
505			<< "\tm4l: " << ZZSystem.M()
506			<< "\tevtfail: " << hex << evtfail << dec
507			<< "\ttrigbits: " << hex << info->triggerBits << dec
508			// << "\ttree: " << inputFiles[q][f]
509			<< endl;
510
511
512
513			//***************************************************************
514			// finish
515			//***************************************************************
516
517			if( !evtfail ) {
518			ret.status.setStatus(SelectionStatus::EVTPASS);
519			ret.Z1leptons.push_back(lepvec[Z1LeptonMinusIndex]);
520			ret.Z1leptons.push_back(lepvec[Z1LeptonPlusIndex]);
521			ret.Z2leptons.push_back(lepvec[Z2LeptonMinusIndex]);
522			ret.Z2leptons.push_back(lepvec[Z2LeptonPlusIndex]);
523			}
524	anlevin	1.4
525
526	khahn	1.1
527			return ret;
528			}
529
530