GenPurposeSkimmerData.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    GenPurposeSkimmerData
// Class:      GenPurposeSkimmerData
//
//
// Original Author:  Nikolaos Rompotis
//         Created:  Thu Oct 16 17:11:55 CEST 2008
//
//

#include "ElectroWeakAnalysis/WENu/interface/GenPurposeSkimmerData.h"

//
//
//
#include "FWCore/Common/interface/TriggerNames.h"
//
//
GenPurposeSkimmerData::GenPurposeSkimmerData(const edm::ParameterSet& ps)

{
//
//   I N P U T      P A R A M E T E R S
//
  // output file name
  outputFile_ = ps.getUntrackedParameter<std::string>("outputfile");
  //
  // Electron Collection
  ElectronCollectionToken_=consumes<pat::ElectronCollection>(ps.getUntrackedParameter<edm::InputTag>("ElectronCollection"));
  //
  // MC:
  //MCCollection_ = ps.getUntrackedParameter<edm::InputTag>("MCCollection");
  //MCCollectionToken_ = consumes<reco::GenParticleCollection>(MCCollection_);
  //MCMatch_Deta_ = ps.getUntrackedParameter<double>("MCMatch_Deta",0.1);
  //MCMatch_Dphi_ = ps.getUntrackedParameter<double>("MCMatch_Dphi",0.35);
  //
  // MET Collections:
  MetCollectionTag_ = ps.getUntrackedParameter<edm::InputTag>("MetCollectionTag");
  MetCollectionToken_ = consumes<reco::CaloMETCollection>(MetCollectionTag_);
  mcMetCollectionToken_ = consumes<pat::METCollection>(ps.getUntrackedParameter<edm::InputTag>("mcMetCollectionTag"));
  t1MetCollectionToken_ = consumes<pat::METCollection>(ps.getUntrackedParameter<edm::InputTag>("t1MetCollectionTag"));
  pfMetCollectionToken_ = consumes<reco::PFMETCollection>(ps.getUntrackedParameter<edm::InputTag>("pfMetCollectionTag"));
  tcMetCollectionToken_ = consumes<reco::METCollection>(ps.getUntrackedParameter<edm::InputTag>("tcMetCollectionTag"));
  //  genMetCollectionToken_ = consumes<reco::GenMETCollection>(ps.getUntrackedParameter<edm::InputTag>("genMetCollectionTag"));
  //
  // HLT parameters:
  // allow info for 2 paths and 2 filters
  // ---------------------------------------------------------------------------
  HLTCollectionE29_= ps.getUntrackedParameter<edm::InputTag>("HLTCollectionE29");
  HLTCollectionE29Token_ = consumes<trigger::TriggerEvent>(HLTCollectionE29_);
  HLTCollectionE31_= ps.getUntrackedParameter<edm::InputTag>("HLTCollectionE31");
  HLTCollectionE31Token_ = consumes<trigger::TriggerEvent>(HLTCollectionE31_);
  HLTTriggerResultsE29_ = ps.getUntrackedParameter<edm::InputTag>("HLTTriggerResultsE29");
  HLTTriggerResultsE29Token_ = consumes<edm::TriggerResults>(HLTTriggerResultsE29_);
  HLTTriggerResultsE31_ = ps.getUntrackedParameter<edm::InputTag>("HLTTriggerResultsE31");
  HLTTriggerResultsE31Token_ = consumes<edm::TriggerResults>(HLTTriggerResultsE31_);
  //HLTPath_ = ps.getUntrackedParameter<std::string>("HLTPath","HLT_Ele15_LW_L1R");
  //HLTFilterType_ =ps.getUntrackedParameter<edm::InputTag>("HLTFilterType");
  //
  // matching HLT objects to electrons
  ProbeHLTObjMaxDR= ps.getUntrackedParameter<double>("ProbeHLTObjMaxDR",0.2);
  //
  // ----------------------------------------------------------------------------
  //
  // detector geometry
  //
  BarrelMaxEta = ps.getUntrackedParameter<double>("BarrelMaxEta");
  EndcapMinEta = ps.getUntrackedParameter<double>("EndcapMinEta");
  EndcapMaxEta = ps.getUntrackedParameter<double>("EndcapMaxEta");
  //
  ctfTracksToken_ = consumes<reco::TrackCollection>(ps.getUntrackedParameter<edm::InputTag>("ctfTracksTag"));
  corHybridscToken_  = consumes<reco::SuperClusterCollection>(ps.getUntrackedParameter<edm::InputTag>("corHybridsc"));
  multi5x5scToken_   = consumes<reco::SuperClusterCollection>(ps.getUntrackedParameter<edm::InputTag>("multi5x5sc"));
  offlineBeamSpotToken_ = consumes<reco::BeamSpot>(edm::InputTag("offlineBeamSpot"));
  pMuonsToken_ = consumes<pat::MuonCollection>(edm::InputTag("selectedLayer1Muons"));

}


GenPurposeSkimmerData::~GenPurposeSkimmerData()
{

   // do anything here that needs to be done at desctruction time
   // (e.g. close files, deallocate resources etc.)

}


//
// member functions
//

// ------------ method called to for each event  ------------
void
GenPurposeSkimmerData::analyze(const edm::Event& evt, const edm::EventSetup& es)
{
  // MC Collection ------------------------------------------------

  //  edm::Handle<reco::GenParticleCollection> pGenPart;
  //  evt.getByToken(MCCollectionToken_, pGenPart);
  //  if ( not  pGenPart.isValid() ) {
  //    std::cout <<"Error! Can't get "<<MCCollection_.label() << std::endl;
  //    return;
  //  }

  //  const reco::GenParticleCollection *McCand = pGenPart.product();

  // GsF Electron Collection ---------------------------------------
  edm::Handle<pat::ElectronCollection> pElectrons;

  try{
    evt.getByToken(ElectronCollectionToken_, pElectrons);
  }
  catch (cms::Exception)
    {
      edm::LogError("")<< "Error! Can't get ElectronCollection by label. ";
    }
  // ***********************************************************************
  // check which trigger has accepted the event ****************************
  // ***********************************************************************
  //
  // path allocation: first 10 paths belong to the low lum menu, the rest
  // in the high lum one
  //
  // Low Luminosity Menu (8e29)
  //
  /*
  edm::Handle<edm::TriggerResults> HLTResultsE29;
  evt.getByToken(HLTTriggerResultsE29Token_, HLTResultsE29);
  if (not HLTResultsE29.isValid()) {
    std::cout << "HLT Results with label: " << HLTTriggerResultsE29_
          << " not found" << std::endl;
    return;
  }
  //
  edm::Handle<trigger::TriggerEvent> pHLTe29;
  evt.getByToken(HLTCollectionE29Token_, pHLTe29);
  if (not pHLTe29.isValid()) {
    std::cout << "HLT Results with label: " << HLTCollectionE29_
          << " not found" << std::endl;
    return;
  }
  //
  int sum = 0;
  //
  for (int iT=0; iT<10; ++iT) {
    event_HLTPath[iT] = 0;
    numberOfHLTFilterObjects[iT] =0;
    //
    const edm::TriggerNames & triggerNames = evt.triggerNames(*HLTResultsE29);
    unsigned int trigger_size = HLTResultsE29->size();
    unsigned int trigger_position = triggerNames.triggerIndex(HLTPath_[iT]);
    if (trigger_position < trigger_size )
      event_HLTPath[iT] = (int) HLTResultsE29->accept(trigger_position);
    //
    numberOfHLTFilterObjects[iT] = 0;
    // check explicitly that the filter is there
    const int nF(pHLTe29->sizeFilters());
    const int filterInd = pHLTe29->filterIndex(HLTFilterType_[iT]);
    if (nF != filterInd) {
      const trigger::Vids& VIDS (pHLTe29->filterIds(filterInd));
      const trigger::Keys& KEYS(pHLTe29->filterKeys(filterInd));
      const int nI(VIDS.size());
      const int nK(KEYS.size());
      numberOfHLTFilterObjects[iT] = (nI>nK)? nI:nK;
    }
    //if (iT==2) // HLT_Ele15_LW_L1R only this trigger is required
      sum += numberOfHLTFilterObjects[iT];
  }
  //
  // High Luminosity Menu (1e31) DISABLED - only low lumi level
  //
  edm::Handle<edm::TriggerResults> HLTResultsE31;
  evt.getByToken(HLTTriggerResultsE31Token_, HLTResultsE31);
  if (not HLTResultsE31.isValid()) {
      std::cout << "HLT Results with label: " << HLTTriggerResultsE31_
            << " not found" << std::endl;
    return;
  }
  edm::Handle<trigger::TriggerEvent> pHLTe31;
  evt.getByToken(HLTCollectionE31Token_, pHLTe31);
  if (not pHLTe31.isValid()) {
    std::cout << "HLT Results with label: " << HLTCollectionE31_
          << " not found" << std::endl;
    return;
  }
  for (int iT=10; iT<25; ++iT) {
    event_HLTPath[iT] = 0;
    numberOfHLTFilterObjects[iT] =0;
    //
    const edm::TriggerNames & triggerNames = evt.triggerNames(*HLTResultsE31);
    unsigned int trigger_size = HLTResultsE31->size();
    unsigned int trigger_position = triggerNames.triggerIndex(HLTPath_[iT]);
    if (trigger_position < trigger_size )
      event_HLTPath[iT] = (int) HLTResultsE31->accept(trigger_position);
    //
    numberOfHLTFilterObjects[iT] = 0;
    // check explicitly that the filter is there
    const int nF(pHLTe31->sizeFilters());
    const int filterInd = pHLTe31->filterIndex(HLTFilterType_[iT]);
    if (nF != filterInd) {
      const trigger::Vids& VIDS (pHLTe31->filterIds(filterInd));
      const trigger::Keys& KEYS(pHLTe31->filterKeys(filterInd));
      const int nI(VIDS.size());
      const int nK(KEYS.size());
      numberOfHLTFilterObjects[iT] = (nI>nK)? nI:nK;
    }
    // not needed
    sum += numberOfHLTFilterObjects[iT];
  }
  if (sum == 0) {
    //std::cout << "No trigger found in this event..." << std::endl;
    return;
  }
  */
  //std::cout << "HLT objects: #" << sum << std::endl;
  // print out the triggers that exist in this event
    // comment this out if you want to see the names of the existing triggers
  edm::Handle<trigger::TriggerEvent> pHLTe29;
  evt.getByToken(HLTCollectionE29Token_, pHLTe29);
  if (not pHLTe29.isValid()){
    std::cout << "Error!!! HLT is missing!" << std::endl;
    return;
  } /*
  else {
    // check explicitly that the filter is there
    const int nF(pHLTe29->sizeFilters());
    for (int filterInd=0; filterInd< nF; ++filterInd) {
      const trigger::Vids& VIDS (pHLTe29->filterIds(filterInd));
      const trigger::Keys& KEYS(pHLTe29->filterKeys(filterInd));
      const int nI(VIDS.size());
      const int nK(KEYS.size());
      int   nObjects = (nI>nK)? nI:nK;
      const edm::InputTag filterTag = pHLTe29->filterTag(filterInd);
      std::cout << "Found filter with name " << filterTag
        << " and #objects: #" << nObjects << std::endl;
    }
  }
    */
  // *********************************************************************
  // MET Collections:
  //
  edm::Handle<reco::CaloMETCollection> caloMET;
  evt.getByToken(MetCollectionToken_, caloMET);
  //
  edm::Handle<pat::METCollection> t1MET;
  evt.getByToken(t1MetCollectionToken_, t1MET);
  //
  edm::Handle<pat::METCollection> mcMET;
  evt.getByToken(mcMetCollectionToken_, mcMET);
  //
  edm::Handle<reco::METCollection> tcMET;
  evt.getByToken(tcMetCollectionToken_, tcMET);
  //
  edm::Handle<reco::PFMETCollection> pfMET;
  evt.getByToken(pfMetCollectionToken_, pfMET);
  //
  //  edm::Handle<reco::GenMETCollection> genMET;
  //  evt.getByToken(genMetCollectionToken_, genMET);
  //
  // initialize the MET variables ........................................
  event_MET     = -99.;   event_MET_phi = -99.;    event_MET_sig = -99.;
  event_mcMET     = -99.;   event_mcMET_phi = -99.;    event_mcMET_sig = -99.;
  event_tcMET   = -99.;   event_tcMET_phi = -99.;  event_tcMET_sig = -99.;
  event_pfMET   = -99.;   event_pfMET_phi = -99.;  event_pfMET_sig = -99.;
  event_t1MET   = -99.;   event_t1MET_phi = -99.;  event_t1MET_sig = -99.;
  //
  //event_genMET  = -99.;   event_genMET_phi= -99.;  event_genMET_sig = -99.;
  //
  // get the values, if they are available
  if ( caloMET.isValid() ) {
    const reco::CaloMETRef MET(caloMET, 0);
    event_MET = MET->et();  event_MET_phi = MET->phi();
    event_MET_sig = MET->mEtSig();
  }
  else {
    std::cout << "caloMET not valid: input Tag: " << MetCollectionTag_
          << std::endl;
  }
  if ( tcMET.isValid() ) {
    const reco::METRef MET(tcMET, 0);
    event_tcMET = MET->et();  event_tcMET_phi = MET->phi();
    event_tcMET_sig = MET->mEtSig();
  }
  if ( pfMET.isValid() ) {
    const reco::PFMETRef MET(pfMET, 0);
    event_pfMET = MET->et();  event_pfMET_phi = MET->phi();
    event_pfMET_sig = MET->mEtSig();
  }
  if ( t1MET.isValid() ) {
    const pat::METRef MET(t1MET, 0);
    event_t1MET = MET->et();  event_t1MET_phi = MET->phi();
    event_t1MET_sig = MET->mEtSig();
  }
  if ( mcMET.isValid() ) {
    const pat::METRef MET(mcMET, 0);
    event_mcMET = MET->et();  event_mcMET_phi = MET->phi();
    event_mcMET_sig = MET->mEtSig();
  }

  //  if ( genMET.isValid() ) {
  //    const reco::GenMETRef MET(genMET, 0);
  //    event_genMET = MET->et();  event_genMET_phi = MET->phi();
  //    event_genMET_sig = MET->mEtSig();
  //  }

  //  std::cout << "t1MET: " << event_t1MET  << " twikiT1MET: "
  //        << event_twikiT1MET  << ", calo="<<event_MET  << std::endl;
  //
  // some supercluster collections ...........................................
  // correcyedHybridSuperClusters
  edm::Handle<reco::SuperClusterCollection> SC1;
  evt.getByToken(corHybridscToken_,SC1);
  const reco::SuperClusterCollection *sc1 = SC1.product();
  // multi5x5SuperClustersWithPreshower
  edm::Handle<reco::SuperClusterCollection> SC2;
  evt.getByToken(multi5x5scToken_,SC2);
  const reco::SuperClusterCollection *sc2 = SC2.product();
  //
  const int n1 =  sc1->size();
  const int n2 =  sc2->size();
  //std::cout << "SC found: hybrid: " << n1 << ", multi5x5: "
  //        << n2 << std::endl;
  // keep details of the 5 highest ET superclusters
  for (int i=0; i<5; ++i) {
    sc_hybrid_et[i] = -9999.;
    sc_hybrid_eta[i] = -9999.;
    sc_hybrid_phi[i] = -9999.;
    //
    sc_multi5x5_et[i] = -9999.;
    sc_multi5x5_eta[i] = -9999.;
    sc_multi5x5_phi[i] = -9999.;
    //
  }
  // sort the energies of the first sc
  std::vector<double> ETsc1;
  std::vector<reco::SuperCluster>::const_iterator sc;
  for (sc = sc1->begin(); sc !=  sc1->end(); ++sc) {
    reco::SuperCluster mySc = *sc;
    double scEt = mySc.energy()/(cosh(mySc.eta()));
    ETsc1.push_back(scEt);

  }
  int *sorted1 = new int[n1];
  double *et1 = new double[n1];
  for (int i=0; i<n1; ++i) {
    et1[i] = ETsc1[i];
  }
  // array sorted now has the indices of the highest ET electrons
  TMath::Sort(n1, et1, sorted1, true);
  // .........................................................................
  std::vector<double> ETsc2;
  for (sc = sc2->begin(); sc !=  sc2->end(); ++sc) {
    reco::SuperCluster mySc = *sc;
    double scEt = mySc.energy()/(cosh(mySc.eta()));
    ETsc2.push_back(scEt);

  }
  int *sorted2 = new int[n2];
  double *et2 = new double[n2];
  for (int i=0; i<n2; ++i) {
    et2[i] = ETsc2[i];
  }
  // array sorted now has the indices of the highest ET electrons
  TMath::Sort(n2, et2, sorted2, true);
  //
  //
  for( int probeSc = 0; probeSc < n1; ++probeSc)
    {
      //std::cout<<"sorted["<< probeIt<< "]=" << sorted[probeIt] << std::endl;
      // break if you have more than the appropriate number of electrons
      if (probeSc >= 5) break;
      //
      int sc_index = sorted1[probeSc];
      std::vector<reco::SuperCluster>::const_iterator
    Rprobe = sc1->begin() + sc_index;
      //
      reco::SuperCluster sc0 = *Rprobe;
      // now keep the relevant stuff:
      sc_hybrid_et[probeSc] =  sc0.energy()/(cosh(sc0.eta()));
      sc_hybrid_eta[probeSc] = sc0.eta();
      sc_hybrid_phi[probeSc] = sc0.phi();
    }
  // .........................................................................
  for( int probeSc = 0; probeSc < n2; ++probeSc)
    {
      //std::cout<<"sorted["<< probeIt<< "]=" << sorted[probeIt] << std::endl;
      // break if you have more than the appropriate number of electrons
      if (probeSc >= 5) break;
      //
      int sc_index = sorted2[probeSc];
      std::vector<reco::SuperCluster>::const_iterator
    Rprobe = sc2->begin() + sc_index;
      //
      reco::SuperCluster sc0 = *Rprobe;
      // now keep the relevant stuff:
      sc_multi5x5_et[probeSc] =  sc0.energy()/(cosh(sc0.eta()));
      sc_multi5x5_eta[probeSc] = sc0.eta();
      sc_multi5x5_phi[probeSc] = sc0.phi();
    }
  delete [] sorted1;  delete [] sorted2;
  delete [] et1;     delete [] et2;
  //  edm::InputTag ctfTracksTag("generalTracks", "", InputTagEnding_);
  edm::Handle<reco::TrackCollection> ctfTracks;
  evt.getByToken(ctfTracksToken_, ctfTracks);
  const reco::TrackCollection *ctf = ctfTracks.product();
  reco::TrackCollection::const_iterator tr;
  const int ntracks =  ctf->size();
  //
  // get the beam spot for the parameter of the track
  edm::Handle<reco::BeamSpot> pBeamSpot;
  evt.getByToken(offlineBeamSpotToken_, pBeamSpot);
  const reco::BeamSpot *bspot = pBeamSpot.product();
  const math::XYZPoint bspotPosition = bspot->position();
  //
  for (int i=0; i<20; ++i) {
    ctf_track_pt[i] = -9999.;
    ctf_track_eta[i] = -9999.;
    ctf_track_phi[i] = -9999.;
    ctf_track_vx[i] = -9999.; ctf_track_vy[i]=-9999.; ctf_track_vz[i] =-9999.;
    ctf_track_tip[i] = -9999.;    ctf_track_tip_bs[i] = -9999.;
  }
  //
  std::vector<double> ETtrack;
  for (tr = ctf->begin(); tr !=  ctf->end(); ++tr) {
    reco::Track mySc = *tr;
    double scEt = mySc.pt();
    ETtrack.push_back(scEt);
  }
  int *sortedTr = new int[ntracks];
  double *etTr = new double[ntracks];
  for (int i=0; i<ntracks; ++i) {
    etTr[i] = ETtrack[i];
  }
  // array sorted now has the indices of the highest ET electrons
  TMath::Sort(ntracks, etTr, sortedTr, true);
  //
  for( int probeSc = 0; probeSc < ntracks; ++probeSc)
    {
      //std::cout<<"sorted["<< probeIt<< "]=" << sorted[probeIt] << std::endl;
      // break if you have more than the appropriate number of electrons
      if (probeSc >= 20) break;
      //
      int sc_index = sortedTr[probeSc];
      std::vector<reco::Track>::const_iterator
    Rprobe = ctf->begin() + sc_index;
      //
      reco::Track sc0 = *Rprobe;
      // now keep the relevant stuff:
      ctf_track_pt[probeSc] =  sc0.pt();
      ctf_track_eta[probeSc] = sc0.eta();
      ctf_track_phi[probeSc] = sc0.phi();
      ctf_track_vx[probeSc] = sc0.vx();
      ctf_track_vy[probeSc] = sc0.vy();
      ctf_track_vz[probeSc] = sc0.vz();
      ctf_track_tip[probeSc] = -sc0.dxy();
      ctf_track_tip_bs[probeSc] = -sc0.dxy(bspotPosition);
    }
  delete [] sortedTr; delete [] etTr;
  //
  // keep 4 of the selectedLayer1Muons for reference
  edm::Handle<pat::MuonCollection> pMuons;
  evt.getByToken(pMuonsToken_, pMuons);
  const pat::MuonCollection *pmuon = pMuons.product();
  pat::MuonCollection::const_iterator muon;
  const int nmuons =  pMuons->size();
  //
  for (int i=0; i<4; ++i) {
    muon_pt[i] = -9999.;
    muon_eta[i] = -9999.;
    muon_phi[i] = -9999.;
    muon_vx[i] = -9999.; muon_vy[i] = -9999.; muon_vz[i] = -9999.;
    muon_tip[i] = -9999.;    muon_tip_bs[i] = -9999.;
  }
  //
  std::vector<double> ETmuons;
  for (muon = pmuon->begin(); muon !=  pmuon->end(); ++muon) {
    pat::Muon mySc = *muon;
    double scEt = mySc.track()->pt();
    ETmuons.push_back(scEt);
  }
  int *sortedMu = new int[nmuons];
  double *etMu = new double[nmuons];
  for (int i=0; i<nmuons; ++i) {
    etMu[i] = ETmuons[i];
  }
  // array sorted now has the indices of the highest ET electrons
  TMath::Sort(nmuons, etMu, sortedMu, true);
  //
  for( int probeSc = 0; probeSc < nmuons; ++probeSc)
    {
      //std::cout<<"sorted["<< probeIt<< "]=" << sorted[probeIt] << std::endl;
      // break if you have more than the appropriate number of electrons
      if (probeSc >= 4) break;
      //
      int sc_index = sortedMu[probeSc];
      std::vector<pat::Muon>::const_iterator
    Rprobe = pmuon->begin() + sc_index;
      //
      pat::Muon sc0 = *Rprobe;
      // now keep the relevant stuff:
      muon_pt[probeSc] =  sc0.track()->pt();
      muon_eta[probeSc] = sc0.track()->eta();
      muon_phi[probeSc] = sc0.track()->phi();
      muon_vx[probeSc] = sc0.track()->vx();
      muon_vy[probeSc] = sc0.track()->vy();
      muon_vz[probeSc] = sc0.track()->vz();
      muon_tip[probeSc] = -sc0.track()->dxy();
      muon_tip_bs[probeSc] = -sc0.track()->dxy(bspotPosition);
    }
  delete [] sortedMu; delete [] etMu;
  //
  if (n1+n2+ntracks == 0) {
    std::cout << "Return: no sc in this event" << std::endl;
    return;
  }
  // /////////////////////////////////////////////////////////////////////////
  // electron details
  const int MAX_PROBES = 4;
  for(int i =0; i < MAX_PROBES; i++){
    probe_ele_eta_for_tree[i] = -99.0;
    probe_ele_et_for_tree[i] = -99.0;
    probe_ele_phi_for_tree[i] = -99.0;
    probe_ele_Xvertex_for_tree[i] = -99.0;
    probe_ele_Yvertex_for_tree[i] = -99.0;
    probe_ele_Zvertex_for_tree[i] = -99.0;
    probe_ele_tip[i] = -999.;

    probe_sc_eta_for_tree[i] = -99.0;
    probe_sc_et_for_tree[i] = -99.0;
    probe_sc_phi_for_tree[i] = -99.0;

    probe_charge_for_tree[i] = -99;
    probe_sc_pass_fiducial_cut[i] = 0;
    probe_classification_index_for_tree[i]=-99;
    //
    // probe isolation values ............
    probe_isolation_value[i] = 999.0;
    probe_iso_user[i] = 999.0;
    probe_ecal_isolation_value[i] = 999;
    probe_ecal_iso_user[i] = 999;
    probe_hcal_isolation_value[i] = 999;
    probe_hcal_iso_user[i] = 999;

    probe_ele_hoe[i]  = 999.;
    probe_ele_shh[i]  = 999.;
    probe_ele_sihih[i] = 999.;
    probe_ele_dhi[i]  = 999.;
    probe_ele_dfi[i]  = 999.;
    probe_ele_eop[i]  = 999.;
    probe_ele_pin[i]  = 999.;
    probe_ele_pout[i] = 999.;
    probe_ele_e5x5[i] = 999.;
    probe_ele_e2x5[i] = 999.;
    probe_ele_e1x5[i] = 999.;

    //
    //
    //for (int j=0; j<25; ++j) {
    //  probe_pass_trigger_cut[i][j]=0;
    //}
    //probe_hlt_matched_dr[i]=0;
    //probe_mc_matched[i] = 0;
    //probe_mc_matched_deta[i] = 999.;
    //probe_mc_matched_dphi[i] = 999.;
    //probe_mc_matched_denergy[i] = 999.;
    //probe_mc_matched_mother[i] = 999;
    //
    //
  }
  const pat::ElectronCollection *electrons= pElectrons.product();


  elec_number_in_event = electrons->size();
  //std::cout << "In this event " << elec_number_in_event <<
  //  " electrons were found" << std::endl;
  //  if (elec_number_in_event == 0) return;

  std::vector<pat::ElectronRef> UniqueElectrons;
  // edm::LogInfo("") << "Starting loop over electrons.";
  int index =0;
  //***********************************************************************
  // NEW METHOD by D WARDROPE implemented 26.05.08 ************************
  //************* DUPLICATE ******  REMOVAL *******************************
  // 02.06.08: due to a bug in the hybrid algorithm that affects detid ****
  //           we change detid matching to superCluster ref matching ******
  for(pat::ElectronCollection::const_iterator
    elec = electrons->begin(); elec != electrons->end();++elec) {
    const pat::ElectronRef  electronRef(pElectrons, index);
    //Remove duplicate electrons which share a supercluster
    pat::ElectronCollection::const_iterator BestDuplicate = elec;
    int index2 = 0;
    for(pat::ElectronCollection::const_iterator
      elec2 = electrons->begin();
    elec2 != electrons->end(); ++elec2)
      {
    if(elec != elec2)
      {
        if( elec->superCluster() == elec2->superCluster())
          {
        if(fabs(BestDuplicate->eSuperClusterOverP()-1.)
           >= fabs(elec2->eSuperClusterOverP()-1.))
          {
            BestDuplicate = elec2;
          }
          }
      }
    ++index2;
      }
    if(BestDuplicate == elec) UniqueElectrons.push_back(electronRef);
    ++index;
  }
  //
  // debugging: store electrons after duplicate removal
  elec_1_duplicate_removal = UniqueElectrons.size();
  //std::cout << "In this event there are " << elec_1_duplicate_removal
  //        << " electrons" << std::endl;
  //
  //
  // duplicate removal is done now:
  //           the electron collection is in UniqueElectrons
  //
  // run over probes - now probe electrons and store
  //
  // the electron collection is now
  // vector<reco::PixelMatchGsfElectronRef>   UniqueElectrons
  std::vector<double> ETs;
  std::vector<pat::ElectronRef>::const_iterator  elec;
  for (elec = UniqueElectrons.begin(); elec !=  UniqueElectrons.end(); ++elec) {
    pat::ElectronRef probeEle;
    probeEle = *elec;
    double probeEt = probeEle->caloEnergy()/(cosh(probeEle->caloPosition().eta()));
    ETs.push_back(probeEt);

  }
  int *sorted = new int[elec_1_duplicate_removal];
  double *et = new double[elec_1_duplicate_removal];
  //std::cout << "Elecs: " << elec_1_duplicate_removal << std::endl;
  for (int i=0; i<elec_1_duplicate_removal; ++i) {
    et[i] = ETs[i];
    //std::cout << "et["<< i << "]=" << et[i] << std::endl;
  }
  // array sorted now has the indices of the highest ET electrons
  TMath::Sort(elec_1_duplicate_removal, et, sorted, true);
  //
  //
  for( int probeIt = 0; probeIt < elec_1_duplicate_removal; ++probeIt)
    {
      //std::cout<<"sorted["<< probeIt<< "]=" << sorted[probeIt] << std::endl;
      // break if you have more than the appropriate number of electrons
      if (probeIt >= MAX_PROBES) break;
      //
      int elec_index = sorted[probeIt];
      std::vector<pat::ElectronRef>::const_iterator
    Rprobe = UniqueElectrons.begin() + elec_index;
      //
      pat::ElectronRef probeEle;
      probeEle = *Rprobe;
      double probeEt = probeEle->caloEnergy()/(cosh(probeEle->caloPosition().eta()));
      probe_sc_eta_for_tree[probeIt] = probeEle->caloPosition().eta();
      probe_sc_phi_for_tree[probeIt] = probeEle->caloPosition().phi();
      probe_sc_et_for_tree[probeIt] = probeEt;
      // fiducial cut ...............................
      if(fabs(probeEle->caloPosition().eta()) < BarrelMaxEta ||
     (fabs(probeEle->caloPosition().eta()) > EndcapMinEta &&
      fabs(probeEle->caloPosition().eta()) < EndcapMaxEta)){
    probe_sc_pass_fiducial_cut[probeIt] = 1;
      }
      //
      probe_charge_for_tree[probeIt] = probeEle->charge();
      probe_ele_eta_for_tree[probeIt] = probeEle->eta();
      probe_ele_et_for_tree[probeIt] = probeEle->et();
      probe_ele_phi_for_tree[probeIt] =probeEle->phi();
      probe_ele_Xvertex_for_tree[probeIt] =probeEle->vx();
      probe_ele_Yvertex_for_tree[probeIt] =probeEle->vy();
      probe_ele_Zvertex_for_tree[probeIt] =probeEle->vz();
      probe_classification_index_for_tree[probeIt] =
    probeEle->classification();
      double ProbeTIP = probeEle->gsfTrack()->d0();
      probe_ele_tip[probeIt] = ProbeTIP;
      // isolation ..................................
      // these are the default values: trk 03, ecal, hcal 04
      // I know that there is a more direct way, but in this way it
      // is clearer what you get each time :P
      probe_isolation_value[probeIt] = probeEle->dr03IsolationVariables().tkSumPt;
      probe_ecal_isolation_value[probeIt] = probeEle->dr04IsolationVariables().ecalRecHitSumEt;
      probe_hcal_isolation_value[probeIt] =
    probeEle->dr04IsolationVariables().hcalDepth1TowerSumEt +
    probeEle->dr04IsolationVariables().hcalDepth2TowerSumEt;
      // one extra isos:
      probe_iso_user[probeIt] = probeEle->dr04IsolationVariables().tkSumPt;
      probe_ecal_iso_user[probeIt] = probeEle->dr03IsolationVariables().ecalRecHitSumEt;
      probe_hcal_iso_user[probeIt] =
    probeEle->dr03IsolationVariables().hcalDepth1TowerSumEt +
    probeEle->dr03IsolationVariables().hcalDepth2TowerSumEt;
      // ele id variables
      double hOverE = probeEle->hadronicOverEm();
      double deltaPhiIn = probeEle->deltaPhiSuperClusterTrackAtVtx();
      double deltaEtaIn = probeEle->deltaEtaSuperClusterTrackAtVtx();
      double eOverP = probeEle->eSuperClusterOverP();
      double pin  = probeEle->trackMomentumAtVtx().R();
      double pout = probeEle->trackMomentumOut().R();
      double sigmaee = probeEle->scSigmaEtaEta();
      double sigma_IetaIeta = probeEle->scSigmaIEtaIEta();
      // correct if in endcaps
      if( fabs (probeEle->caloPosition().eta()) > 1.479 )  {
    sigmaee = sigmaee - 0.02*(fabs(probeEle->caloPosition().eta()) -2.3);
      }
      //
      //double e5x5, e2x5Right, e2x5Left, e2x5Top, e2x5Bottom, e1x5;
      double e5x5, e2x5, e1x5;
      e5x5 = probeEle->scE5x5();
      e1x5 = probeEle->scE1x5();
      e2x5 = probeEle->scE2x5Max();
      //
      // electron ID variables
      probe_ele_hoe[probeIt] = hOverE;
      probe_ele_shh[probeIt] = sigmaee;
      probe_ele_sihih[probeIt] = sigma_IetaIeta;
      probe_ele_dfi[probeIt] = deltaPhiIn;
      probe_ele_dhi[probeIt] = deltaEtaIn;
      probe_ele_eop[probeIt] = eOverP;
      probe_ele_pin[probeIt] = pin;
      probe_ele_pout[probeIt] = pout;
      probe_ele_e5x5[probeIt] = e5x5;
      probe_ele_e2x5[probeIt] = e2x5;
      probe_ele_e1x5[probeIt] = e1x5;

      //
      // HLT filter ------------------------------------------------------
      //
      //
      // low luminosity filters
      /*************************************************************
      for (int filterNum=0; filterNum<10; ++filterNum) {
    int trigger_int_probe = 0;

    //double hlt_matched_dr   = -1.;
    const int nF(pHLTe29->sizeFilters());
    //
    // default (tag) trigger filter
    //
    // find how many relevant
    const int iF = pHLTe29->filterIndex(HLTFilterType_[filterNum]);
    // loop over these objects to see whether they match
    const trigger::TriggerObjectCollection& TOC(pHLTe29->getObjects());
    if (nF != iF) {
      // find how many objects there are
      const trigger::Keys& KEYS(pHLTe29->filterKeys(iF));
      const int nK(KEYS.size());
      for (int iTrig = 0;iTrig <nK; ++iTrig ) {
        const trigger::TriggerObject& TO(TOC[KEYS[iTrig]]);
        //std::cout << "--> filter: "<< HLTFilterType_[filterNum]  <<" TO id: " << TO.id() << std::endl;
        // this is better to be left out: HLT matching is with an HLT object
        // and we don't care what this object is
        //if (abs(TO.id())==11 ) { // demand it to be an electron
        double dr_ele_HLT =
          reco::deltaR(probeEle->eta(), probeEle->phi(), TO.eta(), TO.phi());
        if (fabs(dr_ele_HLT) < ProbeHLTObjMaxDR) {++trigger_int_probe;
        //hlt_matched_dr = dr_ele_HLT;
        }
        //}
      }
    }
    //
    if(trigger_int_probe>0) probe_pass_trigger_cut[probeIt][filterNum] = 1;
    //probe_hlt_matched_dr[probeIt] = hlt_matched_dr;
      }
      // high lumi filters
      for (int filterNum=10; filterNum<25; ++filterNum) {
        int trigger_int_probe = 0;

        //double hlt_matched_dr   = -1.;
        const int nF(pHLTe31->sizeFilters());
        //
        // default (tag) trigger filter
        //
        // find how many relevant
        const int iF = pHLTe31->filterIndex(HLTFilterType_[filterNum]);
        // loop over these objects to see whether they match
        const trigger::TriggerObjectCollection& TOC(pHLTe31->getObjects());
    if (nF != iF) {
      // find how many objects there are
      const trigger::Keys& KEYS(pHLTe31->filterKeys(iF));
      const int nK(KEYS.size());
      for (int iTrig = 0;iTrig <nK; ++iTrig ) {
        const trigger::TriggerObject& TO(TOC[KEYS[iTrig]]);
        //if (abs(TO.id())==11 ) { // demand it to be an electron
        double dr_ele_HLT =
          reco::deltaR(probeEle->eta(), probeEle->phi(), TO.eta(), TO.phi());
        if (fabs(dr_ele_HLT) < ProbeHLTObjMaxDR) {++trigger_int_probe;
        //hlt_matched_dr = dr_ele_HLT;
        }
      }
    }

    //
    if(trigger_int_probe>0) probe_pass_trigger_cut[probeIt][filterNum] = 1;
    //probe_hlt_matched_dr[probeIt] = hlt_matched_dr;
      }
      ******************************************/
      // ------------------------------------------------------------------
      //
      // MC Matching ......................................................
      // check whether these electrons are matched to a MC electron
      /*
      int mc_index = 0;
      int matched = 0; int mother_id = 999;
      double deta_matched = 999.;      double dphi_matched = 999.;
      double denergy_matched = 999.;
      for(reco::GenParticleCollection::const_iterator   McParticle =
        McCand->begin(); McParticle != McCand->end();  ++McParticle)
    {
      // check only for electrons
      if(abs(McParticle->pdgId())==11 && McParticle->status()==1) {
        mc_index++;
        // check whether it matches a gsf electron
        double deta = McParticle->eta() - probeEle->eta();
        double dphi = McParticle->phi() - probeEle->phi();
        if ( fabs(deta) < MCMatch_Deta_  && fabs(dphi) < MCMatch_Dphi_){
          ++matched;
          deta_matched = deta; dphi_matched = dphi;
          denergy_matched = McParticle->energy() - probeEle->caloEnergy();
          // find the mother of the MC electron
          const reco::Candidate *mum;
          bool mother_finder = true;
          if (abs(McParticle->mother()->pdgId()) != 11)
        mum = McParticle->mother();
          else if (abs(McParticle->mother()->mother()->pdgId())!= 11)
        mum = McParticle->mother()->mother();
          else {
        edm::LogInfo("info") << "Going too far to find the mum";
        mother_finder = false;
          }
          if (mother_finder) {
        mother_id = mum->pdgId();
          }
        }
      }
    }
      probe_mc_matched[probeIt] = matched;
      probe_mc_matched_deta[probeIt] = deta_matched;
      probe_mc_matched_dphi[probeIt] = dphi_matched;
      probe_mc_matched_denergy[probeIt] = denergy_matched;
      probe_mc_matched_mother[probeIt] = mother_id;
      */
    }

  probe_tree->Fill();
  ++ tree_fills_;
  delete []  sorted;
  delete []  et;
}


// ------------ method called once each job just before starting event loop  --
void
GenPurposeSkimmerData::beginJob()
{
  //std::cout << "In beginJob()" << std::endl;
  TString filename_histo = outputFile_;
  histofile = new TFile(filename_histo,"RECREATE");
  tree_fills_ = 0;

  probe_tree =  new TTree("probe_tree","Tree to store probe variables");

  //probe_tree->Branch("probe_ele_eta",probe_ele_eta_for_tree,"probe_ele_eta[4]/D");
  //probe_tree->Branch("probe_ele_phi",probe_ele_phi_for_tree,"probe_ele_phi[4]/D");
  //probe_tree->Branch("probe_ele_et",probe_ele_et_for_tree,"probe_ele_et[4]/D");
  probe_tree->Branch("probe_ele_tip",probe_ele_tip,"probe_ele_tip[4]/D");
  probe_tree->Branch("probe_ele_vertex_x",probe_ele_Xvertex_for_tree,
             "probe_ele_vertex_x[4]/D");
  probe_tree->Branch("probe_ele_vertex_y",probe_ele_Yvertex_for_tree,
             "probe_ele_vertex_y[4]/D");
  probe_tree->Branch("probe_ele_vertex_z",probe_ele_Zvertex_for_tree,
             "probe_ele_vertex_z[4]/D");
  probe_tree->Branch("probe_sc_eta",probe_sc_eta_for_tree,"probe_sc_eta[4]/D");
  probe_tree->Branch("probe_sc_phi",probe_sc_phi_for_tree,"probe_sc_phi[4]/D");
  probe_tree->Branch("probe_sc_et",probe_sc_et_for_tree,"probe_sc_et[4]/D");

  // trigger related variables
  //probe_tree->Branch("probe_trigger_cut",probe_pass_trigger_cut,"probe_trigger_cut[4][25]/I");
  //probe_tree->Branch("probe_hlt_matched_dr", probe_hlt_matched_dr,"probe_hlt_matched_dr[4]/D");
  // mc matching to electrons
  //  probe_tree->Branch("probe_mc_matched",probe_mc_matched,"probe_mc_matched[4]/I");
  //probe_tree->Branch("probe_mc_matched_deta",probe_mc_matched_deta,
  //         "probe_mc_matched_deta[4]/D");
  //probe_tree->Branch("probe_mc_matched_dphi",probe_mc_matched_dphi,
  //             "probe_mc_matched_dphi[4]/D");
  //probe_tree->Branch("probe_mc_matched_denergy",probe_mc_matched_denergy,
  //             "probe_mc_matched_denergy[4]/D");
  //probe_tree->Branch("probe_mc_matched_mother",probe_mc_matched_mother,
  //             "probe_mc_matched_mother[4]/I");
  //
  probe_tree->Branch("probe_charge",probe_charge_for_tree,"probe_charge[4]/I");
  //probe_tree->Branch("probe_sc_fiducial_cut",probe_sc_pass_fiducial_cut,
  //             "probe_sc_fiducial_cut[4]/I");



  //probe_tree->Branch("probe_classification",
  //        probe_classification_index_for_tree,"probe_classification[4]/I");
  //
  // Isolation related variables ........................................
  //
  probe_tree->Branch("probe_isolation_value",probe_isolation_value, "probe_isolation_value[4]/D");
  probe_tree->Branch("probe_ecal_isolation_value",probe_ecal_isolation_value, "probe_ecal_isolation_value[4]/D");
  probe_tree->Branch("probe_hcal_isolation_value",probe_hcal_isolation_value,"probe_hcal_isolation_value[4]/D");
  //
  probe_tree->Branch("probe_iso_user",     probe_iso_user,      "probe_iso_user[4]/D");
  probe_tree->Branch("probe_ecal_iso_user",probe_ecal_iso_user, "probe_ecal_iso_user[4]/D");
  probe_tree->Branch("probe_hcal_iso_user",probe_hcal_iso_user, "probe_hcal_iso_user[4]/D");

  //......................................................................
  // Electron ID Related variables .......................................
  probe_tree->Branch("probe_ele_hoe",probe_ele_hoe, "probe_ele_hoe[4]/D");
  //probe_tree->Branch("probe_ele_shh",probe_ele_shh, "probe_ele_shh[4]/D");
  probe_tree->Branch("probe_ele_sihih",probe_ele_sihih,"probe_ele_sihih[4]/D");
  probe_tree->Branch("probe_ele_dfi",probe_ele_dfi, "probe_ele_dfi[4]/D");
  probe_tree->Branch("probe_ele_dhi",probe_ele_dhi, "probe_ele_dhi[4]/D");
  probe_tree->Branch("probe_ele_eop",probe_ele_eop, "probe_ele_eop[4]/D");
  probe_tree->Branch("probe_ele_pin",probe_ele_pin, "probe_ele_pin[4]/D");
  probe_tree->Branch("probe_ele_pout",probe_ele_pout, "probe_ele_pout[4]/D");
  // probe_tree->Branch("probe_ele_e5x5",probe_ele_e5x5, "probe_ele_e5x5[4]/D");
  //probe_tree->Branch("probe_ele_e2x5",probe_ele_e2x5, "probe_ele_e2x5[4]/D");
  //probe_tree->Branch("probe_ele_e1x5",probe_ele_e1x5, "probe_ele_e1x5[4]/D");

  //.......................................................................
  //
  // each entry for each trigger path
  //probe_tree->Branch("event_HLTPath",event_HLTPath,"event_HLTPath[25]/I");
  //probe_tree->Branch("numberOfHLTFilterObjects", numberOfHLTFilterObjects,
  //             "numberOfHLTFilterObjects[25]/I");
  //
  // debugging info:
  //probe_tree->Branch("elec_number_in_event",&elec_number_in_event,"elec_number_in_event/I");
  probe_tree->Branch("elec_1_duplicate_removal",&elec_1_duplicate_removal,"elec_1_duplicate_removal/I");
  //

  // Missing ET in the event
  probe_tree->Branch("event_MET",&event_MET,"event_MET/D");
  probe_tree->Branch("event_MET_phi",&event_MET_phi,"event_MET_phi/D");
  //  probe_tree->Branch("event_MET_sig",&event_MET_sig,"event_MET_sig/D");
  probe_tree->Branch("event_mcMET",&event_mcMET,"event_mcMET/D");
  probe_tree->Branch("event_mcMET_phi",&event_mcMET_phi,"event_mcMET_phi/D");
  //
  probe_tree->Branch("event_tcMET",&event_tcMET,"event_tcMET/D");
  probe_tree->Branch("event_tcMET_phi",&event_tcMET_phi,"event_tcMET_phi/D");
  //  probe_tree->Branch("event_tcMET_sig",&event_tcMET_sig,"event_tcMET_sig/D");

  probe_tree->Branch("event_pfMET",&event_pfMET,"event_pfMET/D");
  probe_tree->Branch("event_pfMET_phi",&event_pfMET_phi,"event_pfMET_phi/D");
  //  probe_tree->Branch("event_pfMET_sig",&event_pfMET_sig,"event_pfMET_sig/D");

  //  probe_tree->Branch("event_genMET",&event_genMET,"event_genMET/D");
  //  probe_tree->Branch("event_genMET_phi",&event_genMET_phi, "event_genMET_phi/D");
  //  probe_tree->Branch("event_genMET_sig",&event_genMET_sig, "event_genMET_sig/D");
  //..... type 1 corrected MET
  probe_tree->Branch("event_t1MET", &event_t1MET, "event_t1MET/D");
  probe_tree->Branch("event_t1MET_phi", &event_t1MET_phi,"event_t1MET_phi/D");
  //probe_tree->Branch("event_t1MET_sig",&event_t1MET_sig,"event_t1MET_sig/D");

  //
  // some sc related variables
  probe_tree->Branch("sc_hybrid_et", sc_hybrid_et, "sc_hybrid_et[5]/D");
  probe_tree->Branch("sc_hybrid_eta", sc_hybrid_eta, "sc_hybrid_eta[5]/D");
  probe_tree->Branch("sc_hybrid_phi", sc_hybrid_phi, "sc_hybrid_phi[5]/D");
  //
  probe_tree->Branch("sc_multi5x5_et",sc_multi5x5_et, "sc_multi5x5_et[5]/D");
  probe_tree->Branch("sc_multi5x5_eta",sc_multi5x5_eta,"sc_multi5x5_eta[5]/D");
  probe_tree->Branch("sc_multi5x5_phi",sc_multi5x5_phi,"sc_multi5x5_phi[5]/D");
  // /////////////////////////////////////////////////////////////////////////
  // general tracks in the event: keep 20 tracks
  probe_tree->Branch("ctf_track_pt",  ctf_track_pt,  "ctf_track_pt[20]/D");
  probe_tree->Branch("ctf_track_eta", ctf_track_eta, "ctf_track_eta[20]/D");
  probe_tree->Branch("ctf_track_phi", ctf_track_phi, "ctf_track_phi[20]/D");
  probe_tree->Branch("ctf_track_vx", ctf_track_vx, "ctf_track_vx[20]/D");
  probe_tree->Branch("ctf_track_vy", ctf_track_vy, "ctf_track_vy[20]/D");
  probe_tree->Branch("ctf_track_vz", ctf_track_vz, "ctf_track_vz[20]/D");
  probe_tree->Branch("ctf_track_tip", ctf_track_tip, "ctf_track_tip[20]/D");
  probe_tree->Branch("ctf_track_tip_bs", ctf_track_tip_bs,
             "ctf_track_tip_bs[20]/D");
  //
  probe_tree->Branch("muon_pt",  muon_pt,  "muon_pt[4]/D");
  probe_tree->Branch("muon_eta", muon_eta, "muon_eta[4]/D");
  probe_tree->Branch("muon_phi", muon_phi, "muon_phi[4]/D");
  probe_tree->Branch("muon_vx", muon_vx, "muon_vx[4]/D");
  probe_tree->Branch("muon_vy", muon_vy, "muon_vy[4]/D");
  probe_tree->Branch("muon_vz", muon_vz, "muon_vz[4]/D");
  probe_tree->Branch("muon_tip", muon_tip, "muon_tip[4]/D");
  probe_tree->Branch("muon_tip_bs", muon_tip_bs, "muon_tip_bs[4]/D");

}

// ------------ method called once each job just after ending the event loop  -
void
GenPurposeSkimmerData::endJob() {
  //std::cout << "In endJob()" << std::endl;
  if (tree_fills_ == 0) {
    std::cout << "Empty tree: no output..." << std::endl;
    return;
  }
  //probe_tree->Print();
  histofile->Write();
  histofile->Close();

}


//define this as a plug-in
DEFINE_FWK_MODULE(GenPurposeSkimmerData);