cutflowSignalRegion.cpp

#include "LLGAnalysis.h"
#include "TLorentzVector.h"

void LLGAnalysis::SetupSignalRegion() {

  // setup the cutflow
  _cutFlow.insert(pair<string,int>("0_NoCut", 0) );
  _cutFlow.insert(pair<string,int>("1_Trigger", 0) );
  _cutFlow.insert(pair<string,int>("2_MuonVeto", 0) );
  _cutFlow.insert(pair<string,int>("3_ElectronVeto", 0) );
  _cutFlow.insert(pair<string,int>("4_HasPVJet", 0) );
  _cutFlow.insert(pair<string,int>("5_HasSV20", 0) );
  _cutFlow.insert(pair<string,int>("6_MET", 0) );
  _cutFlow.insert(pair<string,int>("6a_SVJetPT", 0 ) );
  _cutFlow.insert(pair<string,int>("7_DiJetMass", 0 ) );
  _cutFlow.insert(pair<string,int>("8_SVPVDistance", 0) );
  _cutFlow.insert(pair<string,int>("9_SVJetsDeltaEta", 0 ) );

  // and the histograms
  makeHist( "nBjetAtSV", 5, -0.5, 4.5, "Number of b-jets associated to SV", "Number of SV" );
  makeHist( "mJJSV", 100, 0., 500., "DiJet mass at SV", "Number of Jet Pairs" );
  makeHist( "nJetsSV", 7, -0.5, 6.5, "Number of Jets associated to SV", "Number of SV" );
  makeHist( "SVJet1Pt", 50, 0., 500., "SV Leading Jet pT [GeV]", "Number of SV" );
  makeHist( "SVJet2Pt", 50, 0., 500., "SV 2^{nd} Leading Jet pT [GeV]", "Number of SV" );
  makeHist( "SVJet3Pt", 50, 0., 500., "SV 3^{rd} Leading Jet pT [GeV]", "Number of SV" );
  makeHist( "SVJet4Pt", 50, 0., 500., "SV 4^{th} Leading Jet pT [GeV]", "Number of SV" );
  makeHist( "PVJet1Pt", 50, 0., 1000., "PV Leading Jet pT [GeV]", "Number of Events" );
  makeHist( "PVJet2Pt", 50, 0., 500., "PV 2^{nd} Leading Jet pT [GeV]", "Number of Events" );
  makeHist( "PVJet3Pt", 50, 0., 500., "PV 3^{rd} Leading Jet pT [GeV]", "Number of Events" );
  makeHist( "PVJet4Pt", 50, 0., 500., "PV 4^{th} Leading Jet pT [GeV]", "Number of Events" );
  makeHist( "nJetsTotal", 26, -0.5, 25.5, "Number of Jets with p_{T} > 30 GeV", "Number of Events" );
  makeHist( "BJet1Pt", 50, 0., 500., "Leading B-tagged Jet p_{T} [GeV]", "Number of Events" );
  makeHist( "BJet2Pt", 50, 0., 500., "Subleading B-tagged Jet p_{T} [GeV]", "Number of Events" );
  makeHist( "BJet3Pt", 50, 0., 500., "3^{rd} leading B-tagged Jet p_{T} [GeV]", "Number of Events" );
  makeHist( "BJet4Pt", 50, 0., 500., "4^{th} leading B-tagged Jet p_{T} [GeV]", "Number of Events" );
  makeHist( "JetLeptonDr", 50, 0., 6., "#DeltaR(jet, electron)", "# Jet-Electron Pairs" );
  makeHist( "nPVWithJet75", 10, -0.5, 9.5, "Number of PV with >= 1 Jet > 75 GeV", "# events" );
  makeHist( "nSVWith2Jets30", 10, -0.5, 9.5, "Number of SV with >= 2 Jets > 30 GeV", "# events" );
  makeHist( "mjjvsleadingjetpt", 100, 0., 500., 100, 0., 500., "DiJet Mass at SV [GeV]", "Leading Jet p_{T} at SV [GeV]", "Number of Events", "COLZ" );
  makeHist( "distancePVSV", 40, 0., 40., "Distance between leading PV and SV [mm]", "Number of PV-SV pairs" );
  makeHist( "distancePVSVSignificance", 40, 0., 400., "Significance of distance between leading PV and SV [mm]", "Number of PV-SV pairs" );
  makeHist( "radialDistancePVSV", 40, 0., 40., "Radial distance between leading PV and SV [mm]", "Number of PV-SV pairs");
  makeHist( "zDistancePVSV", 40, 0., 40., "z distance between leading PV and SV [mm]", "Number of PV-SV pairs");
  makeHist( "SVDistanceToBeamPipe", 100, 0., 100., "Distance from Beampipe [mm]", "Number of SV", "HIST" );
  makeHist( "SVDistanceSignificance", 100, 0., 10000., "Significance of Distance to Beampipe", "Number of SV", "HIST");
  makeHist( "SVJetsDPhi", 20, 0., M_PI, "#Delta #phi between SV jets", "Number of SV Jet pairs" );
  makeHist( "SVJetsDR", 20, 0., 10., "#Delta R between SV jets", "Number of SV Jet pairs" );
  makeHist( "SVJetsDEta", 20, 0., 10., "Delta #eta between SV jets", "Number of SV Jet Pairs" );
  makeHist( "SVJetsNConsidered", 71, -0.5, 70.5, "Number of Considered Jet Constituents", "Number of SV Jets", "HIST" );
  makeHist( "SV2Jets_DeltaR", 20, 0., 100.e-3, "Secondary Vertex Radial Resolution", "Number of SV" );
  makeHist( "SV2Jets_DeltaZ", 20, 0., 100.e-3, "Secondary Vertex Z Resolution", "Number of SV" );

  makeHist( "SV2Jets_JetEta", 50, -2.5, 2.5, "SV Jet #eta", "Number of SV Jets" );
  makeHist( "SV2Jets_DeltaEta", 50, -2.5, 2.5, "SV Jet #Delta #eta", "Number of SV jet pairs" );
  makeHist( "SV2Jets_2DJetEta", 50, -2.5, 2.5, 50, -2.5, 2.5, "SV Jet #eta", "SV Jet #eta", "Number of SV Jets", "COLZ" );
  makeHist( "SV2Jets_AverageDistance", 30, 0., 30., "Average Distance at PCA for Jet Constitutents [cm]", "Number of Jets" );
  makeHist( "SV2Jets_RMSDistance", 30, 0., 150., "RMS of Distance at PCA for Jet Constitutents [cm]", "Number of Jets" );
  makeHist( "SV2Jets_AverageVsRMSDistance", 30, 0., 30., 30, 0., 150., "Average Distance at PCA for Jet Constitutents [cm]", "RMS of Distance at PCA for Jet Constitutents [cm]", "Number of Jets", "COLZ" );
  makeHist( "SV2Jets_AverageDistanceVsNConst", 71, -0.5, 70.5, 30, 0., 30., "Number of Considered Jet Constituents", "Average Distance at PCA for Jet Constitutents [cm]", "Number of Jets", "COLZ" );
  makeHist( "SV2Jets_RMSDistanceVsNConst", 71, -0.5, 70.5, 30, 0., 150., "Number of Considered Jet Constituents", "RMS Distance at PCA for Jet Constitutents [cm]", "Number of Jets", "COLZ" );

  makeHist( "SV2Jets_2DAverageDistance", 30, 0., 30., 30, 0., 30., "Average Distance at PCA for Jet Constitutents [cm]", "Average Distance at PCA for Jet Constitutents [cm]", "Number of Jets", "COLZ" );
  makeHist( "SV2Jets_2DRMSDistance", 30, 0., 150., 30, 0., 150., "RMS of Distance at PCA for Jet Constitutents [cm]", "RMS of Distance at PCA for Jet Constitutents [cm]", "Number of Jets", "COLZ" );

  makeHist( "SV2Jets_AbsVertexScore", 20, 1.e-1, 2.e5, "Secondary Vertex Score", "Number of SV" );
  makeHist( "SV2Jets_VertexScoreDivN", 20, 1.e-2, 3.e4, "Secondary Vertex Score / Number of Considered Constitutents", "Number of SV" );
  makeHist( "SV2Jets_VertexScoreDivPt", 20, 1.e-3, 5.e3, "Secondary Vertex Score / Jet pT", "Number of SV" );
  makeHist( "SV2Jets_VertexScoreDivNPt", 20, 1.e-4, 5.e2, "Normalised Secondary Vertex Score", "Number of SV" );

  makeHist( "SV2Jets2D_AbsVertexScore", 50, 100, 2.e4, 50, 1000, 2.e5, "Secondary Vertex Score", "Secondary Vertex Score", "Number of SV", "COLZ" );
  makeHist( "SV2Jets2D_VertexScoreDivN", 50, 10., 3.e3, 50, 1.e2, 3.e4, "Secondary Vertex Score / Number of Considered Constitutents",  "Secondary Vertex Score / Number of Considered Constitutents","Number of SV", "COLZ" );
  makeHist( "SV2Jets2D_VertexScoreDivPt", 50, 10., 5.e2, 50, 10., 5.e3, "Secondary Vertex Score / Jet pT",  "Secondary Vertex Score / Jet pT", "Number of SV", "COLZ" );
  makeHist( "SV2Jets2D_VertexScoreDivNPt", 50, 1., 5.e1, 50, 1., 5.e2, "Normalised Secondary Vertex Score",  "Normalised Secondary Vertex Score", "Number of SV", "COLZ" );

  for( double jptpv = 30.; jptpv <= 401.; jptpv += 5. ) {
    std::vector<double> yield;
    //for( double jptsv = 40.; jptsv <= 401.; jptsv += 10. ) {
    //for( double jptsv = 15.; jptsv <= 201.; jptsv += 5. ) {
    for( double jptsv = 0.; jptsv <= 50.; jptsv += 2. ) {
      yield.push_back( 0. );
    }
    _yields2DOptimisation.push_back( yield );
  }

  passedLogFile.open( "passedEvents.log", ios::out );
  return;
}

void LLGAnalysis::SignalRegionSelection() {


  int IDXFIRSTCUT = -1;
  int IDXSECONDCUT = -1;

  int leadingPV = -1;
  double leadingVertexPt = 0.;


  _cutFlow.at("0_NoCut") += 1;



  bool passTrigger = false;
  for( unsigned int iTrig = 0; iTrig < triggerNames->size(); ++iTrig ) {
    if( triggerNames->at(iTrig).find("HLT_PFMET170_NoiseCleaned") != std::string::npos && triggerBits->at(iTrig) == 1 ) passTrigger = true;
  }

  if( !passTrigger ) return;
  _cutFlow.at("1_Trigger") += 1;


  // lepton veto:
  if( vetoMuons.size() > 0 ) return;
  _cutFlow.at("2_MuonVeto") += 1;


  if( vetoElectrons.size() > 0 ) return;
  _cutFlow.at("3_ElectronVeto") += 1;


  // now assign jets to the vertices:
  vector<int> nJetsToPV( vertex_x->size(), 0 );
  vector<int> nJetsToSV( secVertex_x->size(), 0 );
  vector<vector<int> > idJetsToPV;
  vector<vector<int> > idJetsToSV;
  for( unsigned int iVtx = 0; iVtx < vertex_x->size(); ++iVtx ) {
    vector<int> idx;
    idJetsToPV.push_back( idx );
  }
  for( unsigned int iVtx = 0; iVtx < secVertex_x->size(); ++iVtx ) {
    vector<int> idx;
    idJetsToSV.push_back( idx );
  }


  for( unsigned int iselJet = 0; iselJet < selectedJets.size(); ++iselJet ) {
    int iJet = selectedJets.at(iselJet);
    //calculate jet vertex position:
    //vector<double> error(3,0.);
    vector<double> position(3,0.);
    position.at(0) = recoJet_vertex_x->at(iJet);
    position.at(1) = recoJet_vertex_y->at(iJet);
    position.at(2) = recoJet_vertex_z->at(iJet);
    //vector<double> position = CalculateVertex( recoJet_constVertex_x->at(iJet), recoJet_constVertex_y->at(iJet), recoJet_constVertex_z->at(iJet), recoJet_const_pt->at(iJet), recoJet_const_charge->at(iJet), recoJet_const_closestVertex_d->at(iJet), nCons, weightednCons, error );
    int nMatch = 0;

    //std::cout << "this is jet # " << iJet << " at position " << position.at(0) << " " << position.at(1) << " " << position.at(2) << std::endl;
    for( unsigned int iVtx = 0; iVtx < vertex_x->size(); ++iVtx ) {
      if( fabs(position.at(0) - vertex_x->at(iVtx) ) < 1.e-10 &&
      fabs(position.at(1) - vertex_y->at(iVtx) ) < 1.e-10 &&
      fabs(position.at(2) - vertex_z->at(iVtx) ) < 1.e-10 ) {
        nJetsToPV.at(iVtx) += 1;
        idJetsToPV.at(iVtx).push_back( iJet );
        if( recoJet_pt->at(iJet).at(SYSJET) > leadingVertexPt ) {
          leadingVertexPt = recoJet_pt->at(iJet).at(SYSJET);
          leadingPV = iVtx;
        }
        nMatch += 1;
      }
    }
    for( unsigned int iVtx = 0; iVtx < secVertex_x->size(); ++iVtx ) {
      //std::cout << "checking SV with " << secVertex_x->at(iVtx) << " " << secVertex_y->at(iVtx) << " " << secVertex_z->at(iVtx) << std::endl;
      if( fabs(position.at(0) - secVertex_x->at(iVtx) ) < 1.e-10 &&
      fabs(position.at(1) - secVertex_y->at(iVtx) ) < 1.e-10 &&
      fabs(position.at(2) - secVertex_z->at(iVtx) ) < 1.e-10 ) {

        nJetsToSV.at(iVtx) += 1;
        if( recoJet_pt->at(iJet).at(SYSJET) > JET_PT_CUT_SV ) {
          idJetsToSV.at(iVtx).push_back( iJet );
        }
        nMatch += 1;

      }
    }
    if( nMatch > 1 ) {
      cout << "WARNING! ASSOCIATED JET TO MORE THAN 1 VERTEX ?!" << endl;
    }
  }

  // now count the number of vertices with jets:
  vector<int> PVWithJet;
  vector<int> SVWithJets;
  vector<int> SVWith2Jets;


  double allPVLeadingJetPt = -1.;
  double allSVLeadingJetPt = -1.;
  double allSVLeadingmJJ = -1.;
  for( unsigned int iPV = 0; iPV < vertex_x -> size(); ++iPV ) {
    bool hasJetPV = false;
    double leadingJetPt = 0.;
    for( unsigned int iiJet = 0; iiJet < idJetsToPV.at(iPV).size(); ++iiJet ) {
      int iJet = idJetsToPV.at(iPV).at(iiJet);
      if( recoJet_pt->at(iJet).at(SYSJET) > JET_PT_CUT_PV ) hasJetPV = true;
      if( recoJet_pt->at(iJet).at(SYSJET) > leadingJetPt ) leadingJetPt = recoJet_pt->at(iJet).at(SYSJET);
      if( leadingJetPt > allPVLeadingJetPt ) allPVLeadingJetPt = leadingJetPt;
    }
    if( hasJetPV ) {
      PVWithJet.push_back( iPV );
    }
  }

  if( PVWithJet.size() > 0 ) {
    int counter = -1;
    for( double jptpv = 30.; jptpv <= 401.; jptpv += 5. ) {
      counter++;
      if( allPVLeadingJetPt > jptpv && allPVLeadingJetPt <= jptpv + 5. ) break;
    }
    IDXFIRSTCUT = counter;

    for( unsigned int iSV = 0; iSV < secVertex_x->size(); ++iSV ) {
      if( idJetsToSV.at(iSV).size() > 0 ) SVWithJets.push_back( iSV );
      if( idJetsToSV.at(iSV).size() >= 2 ) SVWith2Jets.push_back( iSV );
    }
  }

  if( SVWith2Jets.size() > 0 ) _outputTree->Fill();


  // do the n-1 plots here:
  // 1st: Leading Lepton pT from PV:
  if( SVWith2Jets.size() > 0 && met->at(SYSMET) > MET_CUT ) {
    for( unsigned int iPV = 0; iPV < vertex_x -> size(); ++iPV ) {
      int idxLeadingJetPV = -1;
      double ptLeadingJetPV = -1.;
      int idxSubLeadingJetPV = -1;
      double ptSubLeadingJetPV = -1;
      int idxThirdLeadingJetPV = -1;
      double ptThirdLeadingJetPV = -1;
      int idxFourthLeadingJetPV = -1;
      double ptFourthLeadingJetPV = -1;
      bool hasJetPV = false;
      for( unsigned int iiJet = 0; iiJet < idJetsToPV.at(iPV).size(); ++iiJet ) {
        int iJet = idJetsToPV.at(iPV).at(iiJet);
        if( recoJet_pt->at(iJet).at(SYSJET) > JET_PT_CUT_PV ) hasJetPV = true;

        if( recoJet_pt->at(iJet).at(SYSJET) > ptLeadingJetPV ) {
          idxFourthLeadingJetPV = idxThirdLeadingJetPV;
          ptFourthLeadingJetPV = ptThirdLeadingJetPV;
          idxThirdLeadingJetPV = idxSubLeadingJetPV;
          ptThirdLeadingJetPV = ptSubLeadingJetPV;
          idxSubLeadingJetPV = idxLeadingJetPV;
          ptSubLeadingJetPV = ptLeadingJetPV;
          idxLeadingJetPV = iJet;
          ptLeadingJetPV = recoJet_pt->at(iJet).at(SYSJET);
        }
        else if( recoJet_pt->at(iJet).at(SYSJET) > ptSubLeadingJetPV ) {
          idxFourthLeadingJetPV = idxThirdLeadingJetPV;
          ptFourthLeadingJetPV = ptThirdLeadingJetPV;
          idxThirdLeadingJetPV = idxSubLeadingJetPV;
          ptThirdLeadingJetPV = ptSubLeadingJetPV;
          idxSubLeadingJetPV = iJet;
          ptSubLeadingJetPV = recoJet_pt->at(iJet).at(SYSJET);
        }
        else if( recoJet_pt->at(iJet).at(SYSJET) > ptThirdLeadingJetPV ) {
          idxFourthLeadingJetPV = idxThirdLeadingJetPV;
          ptFourthLeadingJetPV = ptThirdLeadingJetPV;
          idxThirdLeadingJetPV = iJet;
          ptThirdLeadingJetPV = recoJet_pt->at(iJet).at(SYSJET);
        }
        else if( recoJet_pt->at(iJet).at(SYSJET) > ptThirdLeadingJetPV ) {
          idxFourthLeadingJetPV = iJet;
          ptFourthLeadingJetPV = recoJet_pt->at(iJet).at(SYSJET);
        }
      }
      _histograms1D.at("PVJet1Pt").Fill( ptLeadingJetPV, evtWeight );
      _histograms1D.at("PVJet2Pt").Fill( ptSubLeadingJetPV, evtWeight );
      _histograms1D.at("PVJet3Pt").Fill( ptThirdLeadingJetPV, evtWeight );
      _histograms1D.at("PVJet4Pt").Fill( ptFourthLeadingJetPV, evtWeight );
    }
    _histograms1D.at("nPVWithJet75").Fill( PVWithJet.size(), evtWeight );
  }

  if( PVWithJet.size() >= 1 && met->at(SYSMET) > MET_CUT ) {

    for( unsigned int iSV = 0; iSV < secVertex_x->size(); ++iSV ) {
      _histograms1D.at("nJetsSV").Fill( idJetsToSV.at(iSV).size(), evtWeight );
    }
    _histograms1D.at("nSVWith2Jets30").Fill( SVWith2Jets.size(), evtWeight );
  }


  // and run the selection:
  if( PVWithJet.size() >= 1 ) {
    _cutFlow.at("4_HasPVJet") += 1;

    if( SVWith2Jets.size() > 0 ) {
      _cutFlow.at("5_HasSV20") += 1;

      vector<double> allDistances;
      double maxDist = 0.;
      for( unsigned int iPV = 0; iPV < PVWithJet.size(); ++iPV ) {
        double thispv_x = vertex_x->at(PVWithJet.at(iPV));
        double thispv_y = vertex_y->at(PVWithJet.at(iPV));
        double thispv_z = vertex_z->at(PVWithJet.at(iPV));
        for( unsigned int iSV = 0; iSV < SVWithJets.size(); ++iSV ) {
          double thissv_x = secVertex_x->at(SVWithJets.at(iSV));
          double thissv_y = secVertex_y->at(SVWithJets.at(iSV));
          double thissv_z = secVertex_z->at(SVWithJets.at(iSV));
          double dx = thissv_x - thispv_x;
          double dy = thissv_y - thispv_y;
          double dz = thissv_z - thispv_z;
          double dist = 10.*sqrt(dx*dx + dy*dy + dz*dz );
          allDistances.push_back( dist );
          if( dist > maxDist ) maxDist = dist;
        }
      }

      if (sqrt(met_x->at(SYSMET) * met_x->at(SYSMET) + met_y->at(SYSMET) * met_y->at(SYSMET)) > MET_CUT ) {
        _cutFlow.at("6_MET") += 1;


        bool hasDiJetPair100 = false;
        bool hasSVJetCut = false;
        bool hasSVFarAway = false;

        vector<double> svDistanceToPV;
        for( unsigned int iSV = 0; iSV < secVertex_x->size(); ++iSV ) {
          double dx = secVertex_x->at(iSV) - vertex_x->at(leadingPV);
          double dy = secVertex_y->at(iSV) - vertex_y->at(leadingPV);
          double dz = secVertex_z->at(iSV) - vertex_z->at(leadingPV);
          double dr = 10.*sqrt( dx*dx + dy*dy );
          svDistanceToPV.push_back( dr );
          if( idJetsToSV.at(iSV).size() <= 1 ) continue;
          _histograms1D.at("nJetsSV").Fill( idJetsToSV.at(iSV).size(), evtWeight );

          int idxLeadingJet = -1;
          double ptLeadingJet = -1.;
          int idxSubLeadingJet = -1;
          double ptSubLeadingJet = -1;
          int idxThirdLeadingJet = -1;
          double ptThirdLeadingJet = -1;
          int idxFourthLeadingJet = -1;
          double ptFourthLeadingJet = -1;

          for( unsigned int iJToSV = 0; iJToSV < idJetsToSV.at(iSV).size(); ++iJToSV ) {
            int jIdx = idJetsToSV.at(iSV).at(iJToSV);
            if( recoJet_pt->at(jIdx).at(SYSJET) > ptLeadingJet ) {
              idxFourthLeadingJet = idxThirdLeadingJet;
              ptFourthLeadingJet = ptThirdLeadingJet;
              idxThirdLeadingJet = idxSubLeadingJet;
              ptThirdLeadingJet = ptSubLeadingJet;
              idxSubLeadingJet = idxLeadingJet;
              ptSubLeadingJet = ptLeadingJet;
              idxLeadingJet = jIdx;
              ptLeadingJet = recoJet_pt->at(jIdx).at(SYSJET);
            }
            else if ( recoJet_pt->at(jIdx).at(SYSJET) > ptSubLeadingJet ) {
              idxFourthLeadingJet = idxThirdLeadingJet;
              ptFourthLeadingJet = ptThirdLeadingJet;
              idxThirdLeadingJet = idxSubLeadingJet;
              ptThirdLeadingJet = ptSubLeadingJet;
              idxSubLeadingJet = jIdx;
              ptSubLeadingJet = recoJet_pt->at(jIdx).at(SYSJET);
            }
            else if ( recoJet_pt->at(jIdx).at(SYSJET) > ptThirdLeadingJet ) {
              idxFourthLeadingJet = idxThirdLeadingJet;
              ptFourthLeadingJet = ptThirdLeadingJet;
              idxThirdLeadingJet = jIdx;
              ptThirdLeadingJet = recoJet_pt->at(jIdx).at(SYSJET);
            }
            else if ( recoJet_pt->at(jIdx).at(SYSJET) > ptFourthLeadingJet ) {
              idxFourthLeadingJet = jIdx;
              ptFourthLeadingJet = recoJet_pt->at(jIdx).at(SYSJET);
            }
          }

          _histograms1D.at("SVJet1Pt").Fill( ptLeadingJet, evtWeight );
          _histograms1D.at("SVJet2Pt").Fill( ptSubLeadingJet, evtWeight );
          _histograms1D.at("SVJet3Pt").Fill( ptThirdLeadingJet, evtWeight );
          _histograms1D.at("SVJet4Pt").Fill( ptFourthLeadingJet, evtWeight );
          TLorentzVector p4Jet1, p4Jet2;
          p4Jet1.SetPtEtaPhiM( recoJet_pt->at(idxLeadingJet).at(SYSJET), recoJet_eta->at(idxLeadingJet), recoJet_phi->at(idxLeadingJet), 0. );
          p4Jet2.SetPtEtaPhiM( recoJet_pt->at(idxSubLeadingJet).at(SYSJET), recoJet_eta->at(idxSubLeadingJet), recoJet_phi->at(idxSubLeadingJet), 0. );
          TLorentzVector p4DiJet = p4Jet1 + p4Jet2;
          _histograms1D.at("mJJSV").Fill( p4DiJet.M(), evtWeight );
          if( p4DiJet.M() > MJJ_CUT ) hasDiJetPair100 = true;
          if( ptLeadingJet > LEADING_SV_JET_CUT ) hasSVJetCut = true;
          _histograms2D.at("mjjvsleadingjetpt").Fill( p4DiJet.M(), ptLeadingJet, evtWeight );
          if( ptLeadingJet > allSVLeadingJetPt ) allSVLeadingJetPt = ptLeadingJet;
          if( p4DiJet.M() <= MJJ_CUT && svDistanceToPV.at(iSV) > 5. ) hasSVFarAway = true;
          if( p4DiJet.M() > allSVLeadingmJJ ) allSVLeadingmJJ = p4DiJet.M();
        }
        int counter = -1;
        /*
        for( double jptsv = 40.; jptsv <= 401.; jptsv += 10. ) {
        counter++;
        if( (allSVLeadingJetPt < jptsv && allSVLeadingJetPt >= jptsv - 10.) || allSVLeadingJetPt < 40. ) break;
      }*/
      /*for( double jptsv = 15.; jptsv <= 201.; jptsv += 5. ) {
      counter ++;
      if( (allSVLeadingmJJ < jptsv && allSVLeadingmJJ >= jptsv - 5. ) || allSVLeadingmJJ < 15. ) break;
    }*/
    for( double jptsv = 0.; jptsv <= 50.; jptsv += 2. ) {
      counter++;
      if( maxDist >= jptsv && maxDist < jptsv+2. ) break;
    }
    IDXSECONDCUT = counter;
    for( unsigned int firstCut = IDXFIRSTCUT; firstCut >= 0; --firstCut ) {
      //for( unsigned int secondCut = IDXSECONDCUT; secondCut < _yields2DOptimisation.at(firstCut).size(); ++secondCut ) {
      for( unsigned int secondCut = IDXSECONDCUT; secondCut >=0; --secondCut ) {
        //std::cout << "attempting to attac at " << firstCut << " and " << secondCut << endl;
        _yields2DOptimisation.at(firstCut).at(secondCut) += evtWeight;
        if( secondCut == 0 ) break;
      }
      if( firstCut == 0  ) break;
    }
    if( hasSVJetCut )  return;
    _cutFlow.at("6a_SVJetPT") += 1;
    if( hasDiJetPair100 ) return;
    _cutFlow.at("7_DiJetMass") += 1;
    passedLogFile << RunNumber << " " << LumiBlock << " " << EventNumber << endl;
    if( !hasSVFarAway ) return;
    _cutFlow.at("8_SVPVDistance") += 1;

    bool jetsAreClose = true;

    int leadingPV = -1;
    double leadingJetPt = -1.;
    for( int iiPV = 0; iiPV < PVWithJet.size(); ++iiPV ) {
      int iPV = PVWithJet.at(iiPV);
      for( unsigned int iiJet = 0; iiJet < idJetsToPV.at(iPV).size(); ++iiJet ) {
        int iJet = idJetsToPV.at(iPV).at(iiJet);
        if( leadingJetPt < recoJet_pt->at(iJet).at(SYSJET) ) {
          leadingPV = iPV;
          leadingJetPt = recoJet_pt->at(iJet).at(SYSJET);
        }
      }
    }
    double lpv_x = vertex_x->at(leadingPV);
    double lpv_y = vertex_y->at(leadingPV);
    double lpv_z = vertex_z->at(leadingPV);

    for( unsigned int iSV = 0; iSV < secVertex_x->size(); ++iSV ) {
      if( idJetsToSV.at(iSV).size() <= 1 ) continue;
      if( svDistanceToPV.at(iSV) <= 5. ) continue;
      double sigma_r = 10.*sqrt( secVertex_dx->at(iSV)*secVertex_dx->at(iSV) + secVertex_dy->at(iSV)*secVertex_dy->at(iSV) );
      double sigma_z = 10.* secVertex_dz->at(iSV);
      _histograms1D.at("SV2Jets_DeltaR").Fill( sigma_r, evtWeight );
      _histograms1D.at("SV2Jets_DeltaZ").Fill( sigma_z, evtWeight );
      for( unsigned int iiJet1 = 0; iiJet1 < idJetsToSV.at(iSV).size(); ++iiJet1 ) {
        int iJet1 = idJetsToSV.at(iSV).at(iiJet1);
        _histograms1D.at("SV2Jets_JetEta").Fill( recoJet_eta->at(iJet1), evtWeight );
        _histograms1D.at("SVJetsNConsidered").Fill( recoJet_nConsidered->at(iJet1), evtWeight );
        //_histograms1D.at("SV2Jets_AverageDistance").Fill( recoJet_averageDistance->at(iJet1), evtWeight );
        //_histograms1D.at("SV2Jets_RMSDistance").Fill( recoJet_rmsDistance->at(iJet1), evtWeight );
        //_histograms2D.at("SV2Jets_AverageVsRMSDistance").Fill( recoJet_averageDistance->at(iJet1), recoJet_rmsDistance->at(iJet1), evtWeight );
        //_histograms2D.at("SV2Jets_AverageDistanceVsNConst").Fill( recoJet_nConsidered->at(iJet1), recoJet_averageDistance->at(iJet1), evtWeight );
        //_histograms2D.at("SV2Jets_RMSDistanceVsNConst").Fill( recoJet_nConsidered->at(iJet1), recoJet_rmsDistance->at(iJet1), evtWeight );

        //double normScore = recoJet_vertex_score->at(iJet1);// / (float)recoJet_nConsidered->at(iJet1);

        //normScore /= recoJet_pt->at(iJet1);
        //_histograms1D.at("SV2Jets_AbsVertexScore").Fill( normScore, evtWeight );
        //_histograms1D.at("SV2Jets_VertexScoreDivN").Fill( normScore/(float)recoJet_nConsidered->at(iJet1), evtWeight );
        //_histograms1D.at("SV2Jets_VertexScoreDivPt").Fill( normScore/recoJet_pt->at(iJet1), evtWeight );
        //_histograms1D.at("SV2Jets_VertexScoreDivNPt").Fill( normScore/recoJet_pt->at(iJet1)/(float)recoJet_nConsidered->at(iJet1), evtWeight );
        for( unsigned int iiJet2 = iiJet1+1; iiJet2 < idJetsToSV.at(iSV).size(); ++iiJet2 ) {
          int iJet2 = idJetsToSV.at(iSV).at(iiJet2);
          /*
          double maxDist = ( recoJet_averageDistance->at(iJet1) > recoJet_averageDistance->at(iJet2) ) ? recoJet_averageDistance->at(iJet1) : recoJet_averageDistance->at(iJet2);
          double minDist = ( recoJet_averageDistance->at(iJet1) > recoJet_averageDistance->at(iJet2) ) ? recoJet_averageDistance->at(iJet2) : recoJet_averageDistance->at(iJet1);

          double maxRMS = ( recoJet_rmsDistance->at(iJet1) > recoJet_rmsDistance->at(iJet2) ) ? recoJet_rmsDistance->at(iJet1) : recoJet_rmsDistance->at(iJet2);
          double minRMS = ( recoJet_rmsDistance->at(iJet1) > recoJet_rmsDistance->at(iJet2) ) ? recoJet_rmsDistance->at(iJet2) : recoJet_rmsDistance->at(iJet1);
          _histograms2D.at("SV2Jets_2DAverageDistance").Fill( maxDist, minDist, evtWeight );
          _histograms2D.at("SV2Jets_2DRMSDistance").Fill( maxRMS, minRMS, evtWeight );
          */
          double maxEta = ( recoJet_eta->at(iJet1) > recoJet_eta->at(iJet2) ) ? recoJet_eta->at(iJet1) : recoJet_eta->at(iJet2);
          double minEta = ( recoJet_eta->at(iJet1) > recoJet_eta->at(iJet2) ) ? recoJet_eta->at(iJet2) : recoJet_eta->at(iJet1);
          _histograms2D.at("SV2Jets_2DJetEta").Fill( maxEta, minEta, evtWeight );
          _histograms1D.at("SV2Jets_DeltaEta").Fill( fabs( maxEta - minEta), evtWeight );
          /*
          double normScore2 = recoJet_vertex_score->at(iJet2);
          std::vector<double> vNormScore;
          std::vector<double> vNormScoreDivN;
          std::vector<double> vNormScoreDivPt;
          std::vector<double> vNormScoreDivNPt;
          vNormScore.push_back( normScore );
          vNormScore.push_back( normScore2 );
          vNormScoreDivN.push_back( normScore/(float)recoJet_nConsidered->at(iJet1) );
          vNormScoreDivN.push_back( normScore2/(float)recoJet_nConsidered->at(iJet2) );
          vNormScoreDivPt.push_back( normScore/recoJet_pt->at(iJet1) );
          vNormScoreDivPt.push_back( normScore2/recoJet_pt->at(iJet2) );
          vNormScoreDivNPt.push_back( normScore/(float)recoJet_nConsidered->at(iJet1)/recoJet_pt->at(iJet1) );
          vNormScoreDivNPt.push_back( normScore2/(float)recoJet_nConsidered->at(iJet2)/recoJet_pt->at(iJet2) );

          std::sort( vNormScore.begin(), vNormScore.end() );
          std::sort( vNormScoreDivN.begin(), vNormScoreDivN.end() );
          std::sort( vNormScoreDivPt.begin(), vNormScoreDivPt.end() );
          std::sort( vNormScoreDivNPt.begin(), vNormScoreDivNPt.end() );
          _histograms2D.at("SV2Jets2D_AbsVertexScore").Fill( vNormScore.at(0), vNormScore.at(1), evtWeight );
          _histograms2D.at("SV2Jets2D_VertexScoreDivN").Fill( vNormScoreDivN.at(0), vNormScoreDivN.at(1), evtWeight );
          _histograms2D.at("SV2Jets2D_VertexScoreDivPt").Fill( vNormScoreDivPt.at(0), vNormScoreDivPt.at(1), evtWeight );
          _histograms2D.at("SV2Jets2D_VertexScoreDivNPt").Fill( vNormScoreDivNPt.at(0), vNormScoreDivNPt.at(1), evtWeight );
          */

          double dphi = fabs( recoJet_phi->at(iJet1) - recoJet_phi->at(iJet2) );
          if( dphi > M_PI ) dphi = 2*M_PI-dphi;
          double deta = fabs( recoJet_eta->at(iJet1) - recoJet_eta->at(iJet2) );
          if( deta > 0.5 ) jetsAreClose = false;
          double dr = sqrt( dphi*dphi + deta*deta );
          _histograms1D.at("SVJetsDPhi").Fill( dphi, evtWeight );
          _histograms1D.at("SVJetsDEta").Fill( deta, evtWeight );
          _histograms1D.at("SVJetsDR").Fill( dr, evtWeight );
        }
      }

      double sv_x = secVertex_x->at( iSV );
      double sv_y = secVertex_y->at( iSV );
      double sv_z = secVertex_z->at( iSV );
      double dx = sv_x - lpv_x;
      double dy = sv_y - lpv_y;
      double dz = sv_z - lpv_z;
      double distance = sqrt( dx*dx + dy*dy + dz*dz );
      _histograms1D.at("distancePVSV").Fill( 10.*distance, evtWeight );
      _histograms1D.at("radialDistancePVSV").Fill( 10.*sqrt(dx*dx+dy*dy), evtWeight );
      _histograms1D.at("zDistancePVSV").Fill( 10.*fabs(dz), evtWeight );
      double err_dx = sqrt( secVertex_dx->at(iSV)*secVertex_dx->at(iSV) + vertex_dx->at(leadingPV)*vertex_dx->at(leadingPV) );
      double err_dy = sqrt( secVertex_dy->at(iSV)*secVertex_dy->at(iSV) + vertex_dy->at(leadingPV)*vertex_dy->at(leadingPV) );
      double err_dz = sqrt( secVertex_dz->at(iSV)*secVertex_dz->at(iSV) + vertex_dz->at(leadingPV)*vertex_dz->at(leadingPV) );
      double err_dist = err_dx*err_dx*dx*dx + err_dy*err_dy*dy*dy + err_dz*err_dz*dz*dz;
      err_dist = err_dist/distance/distance;
      err_dist = sqrt( err_dist );
      _histograms1D.at("distancePVSVSignificance").Fill( distance/err_dist, evtWeight );

      _histograms1D.at("SVDistanceToBeamPipe").Fill( 10.*sqrt(sv_x*sv_x+sv_y*sv_y), evtWeight );
      double err_x = secVertex_dx->at(iSV);
      double err_y = secVertex_dy->at(iSV);
      double totErr = sv_x*sv_x*err_x*err_x + sv_y*sv_y*err_y*err_y;
      totErr = totErr/( sv_x*sv_x + sv_y*sv_y);
      totErr=  sqrt(totErr);
      _histograms1D.at("SVDistanceSignificance").Fill( sqrt(sv_x*sv_x+sv_y*sv_y)/totErr, evtWeight );


    }
    if( jetsAreClose ) return;
    _cutFlow.at("9_SVJetsDeltaEta") += 1;

  }
}
}
return;
}