/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/Calibration/IsolatedParticles/plugins/IsoTrig.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    IsoTrig
// Class:      IsoTrig
// 
//
// Original Author:  Ruchi Gupta
//         Created:  Fri May 25 12:02:48 CDT 2012
// $Id: IsoTrig.cc,v 1.1 2012/10/10 15:04:59 sunanda Exp $
//
//
#include "IsoTrig.h"

//Tracks
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/TrackReco/interface/HitPattern.h"
#include "DataFormats/TrackReco/interface/TrackBase.h"
// Vertices
#include "DataFormats/VertexReco/interface/VertexFwd.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/BeamSpot/interface/BeamSpot.h"

//Triggers
#include "DataFormats/Common/interface/TriggerResults.h"
#include "DataFormats/HLTReco/interface/TriggerEvent.h"
#include "DataFormats/HLTReco/interface/TriggerObject.h"
#include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
#include "DataFormats/Luminosity/interface/LumiDetails.h"

#include "FWCore/Common/interface/TriggerNames.h"

#include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"
#include "Calibration/IsolatedParticles/interface/ChargeIsolation.h"
#include "Calibration/IsolatedParticles/interface/eCone.h"
#include "Calibration/IsolatedParticles/interface/eECALMatrix.h"

#include "MagneticField/Engine/interface/MagneticField.h"
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
#include "Geometry/CaloEventSetup/interface/CaloTopologyRecord.h"
#include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
#include "Geometry/CaloTopology/interface/HcalTopology.h"
#include "Geometry/CaloTopology/interface/CaloTopology.h"
#include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
#include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h"
#include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
#include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"

IsoTrig::IsoTrig(const edm::ParameterSet& iConfig) : changed(false) {
   //now do whatever initialization is needed
  Det                                 = iConfig.getParameter<std::string>("Det");
  verbosity                           = iConfig.getUntrackedParameter<int>("Verbosity",0);
  theTrackQuality                     = iConfig.getUntrackedParameter<std::string>("TrackQuality","highPurity");
  reco::TrackBase::TrackQuality trackQuality_=reco::TrackBase::qualityByName(theTrackQuality);
  selectionParameters.minPt           = iConfig.getUntrackedParameter<double>("MinTrackPt", 10.0);
  selectionParameters.minQuality      = trackQuality_;
  selectionParameters.maxDxyPV        = iConfig.getUntrackedParameter<double>("MaxDxyPV", 0.2);
  selectionParameters.maxDzPV         = iConfig.getUntrackedParameter<double>("MaxDzPV",  5.0);
  selectionParameters.maxChi2         = iConfig.getUntrackedParameter<double>("MaxChi2",  5.0);
  selectionParameters.maxDpOverP      = iConfig.getUntrackedParameter<double>("MaxDpOverP",  0.1);
  selectionParameters.minOuterHit     = iConfig.getUntrackedParameter<int>("MinOuterHit", 4);
  selectionParameters.minLayerCrossed = iConfig.getUntrackedParameter<int>("MinLayerCrossed", 8);
  selectionParameters.maxInMiss       = iConfig.getUntrackedParameter<int>("MaxInMiss", 0);
  selectionParameters.maxOutMiss      = iConfig.getUntrackedParameter<int>("MaxOutMiss", 0);
  dr_L1                               = iConfig.getUntrackedParameter<double>("IsolationL1",1.0);
  a_coneR                             = iConfig.getUntrackedParameter<double>("ConeRadius",34.98);
  a_charIsoR                          = a_coneR + 28.9;
  a_neutIsoR                          = a_charIsoR*0.726;
  a_mipR                              = iConfig.getUntrackedParameter<double>("ConeRadiusMIP",14.0);
  a_neutR1                            = iConfig.getUntrackedParameter<double>("ConeRadiusNeut1",21.0);
  a_neutR2                            = iConfig.getUntrackedParameter<double>("ConeRadiusNeut2",29.0);
  cutMip                              = iConfig.getUntrackedParameter<double>("MIPCut", 1.0);
  cutCharge                           = iConfig.getUntrackedParameter<double>("ChargeIsolation",  2.0);
  cutNeutral                          = iConfig.getUntrackedParameter<double>("NeutralIsolation",  2.0);
  minRunNo                            = iConfig.getUntrackedParameter<int>("minRun");
  maxRunNo                            = iConfig.getUntrackedParameter<int>("maxRun");
  if (verbosity>=0) {
    std::cout <<"Parameters read from config file \n" 
              <<"\t minPt "           << selectionParameters.minPt   
              <<"\t theTrackQuality " << theTrackQuality
              <<"\t minQuality "      << selectionParameters.minQuality
              <<"\t maxDxyPV "        << selectionParameters.maxDxyPV          
              <<"\t maxDzPV "         << selectionParameters.maxDzPV          
              <<"\t maxChi2 "         << selectionParameters.maxChi2          
              <<"\t maxDpOverP "      << selectionParameters.maxDpOverP
              <<"\t minOuterHit "     << selectionParameters.minOuterHit
              <<"\t minLayerCrossed " << selectionParameters.minLayerCrossed
              <<"\t maxInMiss "       << selectionParameters.maxInMiss
              <<"\t maxOutMiss "      << selectionParameters.maxOutMiss
              <<"\t a_coneR "         << a_coneR          
              <<"\t a_charIsoR "      << a_charIsoR          
              <<"\t a_neutIsoR "      << a_neutIsoR          
              <<"\t a_mipR "          << a_mipR 
              <<"\t a_neutR "         << a_neutR1 << ":" << a_neutR2
              <<"\t cuts (MIP "       << cutMip << " : Charge " << cutCharge
              <<" : Neutral "         << cutNeutral << ")"
              << std::endl;
  }
}

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

}

void IsoTrig::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {

  if (verbosity%10 > 0) std::cout << "Event starts====================================" << std::endl;
  int RunNo = iEvent.id().run();
  int EvtNo = iEvent.id().event();
  int Lumi  = iEvent.luminosityBlock();
  int Bunch = iEvent.bunchCrossing();

  edm::ESHandle<MagneticField> bFieldH;
  iSetup.get<IdealMagneticFieldRecord>().get(bFieldH);
  const MagneticField *bField = bFieldH.product();

  // get handles to calogeometry and calotopology
  edm::ESHandle<CaloGeometry> pG;
  iSetup.get<CaloGeometryRecord>().get(pG);
  const CaloGeometry* geo = pG.product();

  edm::Handle<reco::TrackCollection> trkCollection;
  iEvent.getByLabel("generalTracks", trkCollection);
  reco::TrackCollection::const_iterator trkItr;

  edm::InputTag lumiProducer("LumiProducer", "", "RECO");
  edm::Handle<LumiDetails> Lumid;
  iEvent.getLuminosityBlock().getByLabel("lumiProducer",Lumid); 
  float mybxlumi=-1;
  if (Lumid.isValid()) 
    mybxlumi=Lumid->lumiValue(LumiDetails::kOCC1,iEvent.bunchCrossing())*6.37;
  if (verbosity%10 > 0)
    std::cout << "RunNo " << RunNo << " EvtNo " << EvtNo << " Lumi " << Lumi 
              << " Bunch " << Bunch << " mybxlumi " << mybxlumi << std::endl;

  edm::InputTag triggerEvent_ ("hltTriggerSummaryAOD","","HLT");
  trigger::TriggerEvent triggerEvent;
  edm::Handle<trigger::TriggerEvent> triggerEventHandle;
  iEvent.getByLabel(triggerEvent_,triggerEventHandle);
  if (!triggerEventHandle.isValid())
    std::cout << "Error! Can't get the product "<< triggerEvent_.label() << std::endl;
  else {
    triggerEvent = *(triggerEventHandle.product());
  
    const trigger::TriggerObjectCollection& TOC(triggerEvent.getObjects());
    edm::InputTag theTriggerResultsLabel ("TriggerResults","","HLT");
    edm::Handle<edm::TriggerResults> triggerResults;
    iEvent.getByLabel( theTriggerResultsLabel, triggerResults);
    char TrigName[50];
    sprintf(TrigName, "HLT_IsoTrack%s", Det.c_str());
    if (triggerResults.isValid()) {
      std::vector<std::string> modules;
      h_nHLT->Fill(triggerResults->size());
      const edm::TriggerNames & triggerNames = iEvent.triggerNames(*triggerResults);

      const std::vector<std::string> & triggerNames_ = triggerNames.triggerNames();
      int hlt(-1), preL1(-1), preHLT(-1), prescale(-1);
      for (unsigned int i=0; i<triggerResults->size(); i++) {
        unsigned int triggerindx = hltConfig_.triggerIndex(triggerNames_[i]);
        const std::vector<std::string>& moduleLabels(hltConfig_.moduleLabels(triggerindx));
        if (triggerNames_[i].find(TrigName)!=std::string::npos) {
          const std::pair<int,int> prescales(hltConfig_.prescaleValues(iEvent,iSetup,triggerNames_[i]));
          hlt = triggerResults->accept(i);
          preL1  = prescales.first;
          preHLT = prescales.second;
          prescale = preL1*preHLT;
          if (verbosity%10 > 0)
            std::cout << triggerNames_[i] << " accept " << hlt << " preL1 " 
                      << preL1 << " preHLT " << preHLT << std::endl;
          if (hlt>0) {
            std::vector<math::XYZTLorentzVector> vec[3];
            if (TrigList.find(RunNo) != TrigList.end() ) {
              TrigList[RunNo] += 1;
            } else {
              TrigList.insert(std::pair<unsigned int, unsigned int>(RunNo,1));
              TrigPreList.insert(std::pair<unsigned int, std::pair<int, int>>(RunNo,prescales));
            }
            //loop over all trigger filters in event (i.e. filters passed)
            for (unsigned int ifilter=0; ifilter<triggerEvent.sizeFilters(); ++ifilter) {  
              std::vector<int> Keys;
              std::string label = triggerEvent.filterTag(ifilter).label();
              //loop over keys to objects passing this filter
              for (unsigned int imodule=0; imodule<moduleLabels.size(); imodule++) {
                if (label.find(moduleLabels[imodule]) != std::string::npos) {
                  if (verbosity%10 > 0) std::cout << "FILTERNAME " << label << std::endl;
                  for (unsigned int ifiltrKey=0; ifiltrKey<triggerEvent.filterKeys(ifilter).size(); ++ifiltrKey) {
                    Keys.push_back(triggerEvent.filterKeys(ifilter)[ifiltrKey]);
                    const trigger::TriggerObject& TO(TOC[Keys[ifiltrKey]]);
                    math::XYZTLorentzVector v4(TO.px(), TO.py(), TO.pz(), TO.energy());
                    if(label.find("L2Filter") != std::string::npos) {
                      vec[1].push_back(v4);
                    } else if (label.find("L3Filter") != std::string::npos) {
                      vec[2].push_back(v4);
                    } else {
                      vec[0].push_back(v4);
                      h_L1ObjEnergy->Fill(TO.energy());
                    }
                    if (verbosity%10 > 0)
                      std::cout << "key " << ifiltrKey << " : pt " << TO.pt() << " eta " << TO.eta() << " phi " << TO.phi() << " mass " << TO.mass() << " Id " << TO.id() << std::endl;
                  }
                }
              }
            }
            h_nL3Objs  -> Fill(vec[2].size());

            for (int j=0; j<3; j++) {
              for (unsigned int k=0; k<vec[j].size(); k++) {
                if (verbosity%10 > 0) std::cout << "vec[" << j << "][" << k << "] pt " << vec[j][k].pt() << " eta " << vec[j][k].eta() << " phi " << vec[j][k].phi() << std::endl;
                fillHist(j, vec[j][k]);
              }
            }

            double deta, dphi, dr;
            for (int lvl=1; lvl<3; lvl++) {
              for (unsigned int i=0; i<vec[lvl].size(); i++) {
                deta = dEta(vec[0][0],vec[lvl][i]);
                dphi = dPhi(vec[0][0],vec[lvl][i]);
                dr   = dR(vec[0][0],vec[lvl][i]);
                if (verbosity%10 > 0) std::cout << "lvl " <<lvl << " i " << i << " deta " << deta << " dphi " << dphi << " dR " << dr << std::endl;
                h_dEtaL1[lvl-1] -> Fill(deta);
                h_dPhiL1[lvl-1] -> Fill(dphi);
                h_dRL1[lvl-1]   -> Fill(std::sqrt(dr));
              }
            }

            math::XYZTLorentzVector mindRvec;
            double mindR;
            for (unsigned int k=0; k<vec[2].size(); ++k) {
              mindR=999.9;
              if (verbosity%10 > 0) std::cout << "L3obj: pt " << vec[2][i].pt() << " eta " << vec[2][i].eta() << " phi " << vec[2][i].phi() << std::endl;
              for (unsigned int j=0; j<vec[1].size(); j++) {
                dr   = dR(vec[2][k],vec[1][j]);
                if(dr<mindR) {
                  mindR=dr;
                  mindRvec=vec[1][j];
                }
              }
              fillDifferences(0, vec[2][k], mindRvec, (verbosity%10 >0));
              if(mindR < 0.03) {
                fillDifferences(1, vec[2][k], mindRvec, (verbosity%10 >0));
                fillHist(6, mindRvec);
                fillHist(8, vec[2][k]);
              } else {
                fillDifferences(2, vec[2][k], mindRvec, (verbosity%10 >0));
                fillHist(7, mindRvec);
                fillHist(9, vec[2][k]);
              }
                      
              mindR=999.9;
              if (verbosity%10 > 0) std::cout << "vec[2][k].eta() " << vec[2][k].eta() << " vec[k][0].phi " << vec[2][k].phi() << std::endl;
              reco::TrackCollection::const_iterator goodTk = trkCollection->end();
              if (trkCollection.isValid()) {
                double mindP = 9999.9;
                for (trkItr=trkCollection->begin(); 
                     trkItr!=trkCollection->end(); trkItr++) {
                  math::XYZTLorentzVector v4(trkItr->px(), trkItr->py(), 
                                             trkItr->pz(), trkItr->p());
                  double deltaR = dR(v4, vec[2][k]);
                  double dp     = std::abs(v4.r()/vec[2][k].r()-1.0);
                  if (deltaR<mindR) {
                    mindR    = deltaR;
                    mindP    = dp;
                    mindRvec = v4;
                    goodTk   = trkItr;
                  }
                  if ((verbosity/10)%10>1 && deltaR<1.0) {
                    std::cout << "track: pt " << v4.pt() << " eta " << v4.eta()
                              << " phi " << v4.phi() << " DR " << deltaR 
                              << std::endl;
                  }
                }
                if (mindR < 0.03 && mindP > 0.1) {
                  for (trkItr=trkCollection->begin(); 
                       trkItr!=trkCollection->end(); trkItr++) {
                    math::XYZTLorentzVector v4(trkItr->px(), trkItr->py(), 
                                               trkItr->pz(), trkItr->p());
                    double deltaR = dR(v4, vec[2][k]);
                    double dp     = std::abs(v4.r()/vec[2][k].r()-1.0);
                    if (dp<mindP && deltaR<0.03) {
                      mindR    = deltaR;
                      mindP    = dp;
                      mindRvec = v4;
                      goodTk   = trkItr;
                    }
                  }
                }
                fillDifferences(3, vec[2][k], mindRvec, (verbosity%10 >0));
                fillHist(3, mindRvec);
                if(mindR < 0.03) {
                  fillDifferences(4, vec[2][k], mindRvec, (verbosity%10 >0));
                  fillHist(4, mindRvec);
                } else {
                  fillDifferences(5, vec[2][k], mindRvec, (verbosity%10 >0));
                  fillHist(5, mindRvec);
                }

                //Define the best vertex
                edm::Handle<reco::VertexCollection> recVtxs;
                iEvent.getByLabel("offlinePrimaryVertices",recVtxs);  
                // Get the beamspot
                edm::Handle<reco::BeamSpot> beamSpotH;
                iEvent.getByLabel("offlineBeamSpot", beamSpotH);
                math::XYZPoint leadPV(0,0,0);
                if (recVtxs->size()>0 && !((*recVtxs)[0].isFake())) {
                  leadPV = math::XYZPoint( (*recVtxs)[0].x(),(*recVtxs)[0].y(), (*recVtxs)[0].z() );
                } else if (beamSpotH.isValid()) {
                  leadPV = beamSpotH->position();
                }
                if ((verbosity/100)%10>0) {
                  std::cout << "Primary Vertex " << leadPV;
                  if (beamSpotH.isValid()) std::cout << " Beam Spot " 
                                                     << beamSpotH->position();
                  std::cout << std::endl;
                }

                std::vector<spr::propagatedTrackDirection> trkCaloDirections;
                spr::propagateCALO(trkCollection, geo, bField, theTrackQuality, trkCaloDirections, ((verbosity/100)%10>2));
                std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
  
                edm::Handle<EcalRecHitCollection> barrelRecHitsHandle;
                edm::Handle<EcalRecHitCollection> endcapRecHitsHandle;
                iEvent.getByLabel("ecalRecHit","EcalRecHitsEB",barrelRecHitsHandle);
                iEvent.getByLabel("ecalRecHit","EcalRecHitsEE",endcapRecHitsHandle);

                unsigned int nTracks=0, ngoodTk=0, nselTk=0;
                int          ieta=999;
                for (trkDetItr = trkCaloDirections.begin(); trkDetItr != trkCaloDirections.end(); trkDetItr++,nTracks++){
                  bool l3Track  = (trkDetItr->trkItr == goodTk);
                  const reco::Track* pTrack = &(*(trkDetItr->trkItr));
                  math::XYZTLorentzVector v4(pTrack->px(), pTrack->py(), 
                                             pTrack->pz(), pTrack->p());
                  bool selectTk = spr::goodTrack(pTrack,leadPV,selectionParameters,((verbosity/100)%10>2));
                  double eMipDR=9999., e_inCone=0, conehmaxNearP=0, mindR=999.9;
                  if (trkDetItr->okHCAL) {
                    HcalDetId detId = (HcalDetId)(trkDetItr->detIdHCAL);
                    ieta = detId.ieta();
                  }
                  for (unsigned k=0; k<vec[0].size(); ++k) {
                    double deltaR = dR(v4, vec[0][k]);
                    if (deltaR<mindR) mindR = deltaR;
                  }
                  if (selectTk && trkDetItr->okECAL && trkDetItr->okHCAL) {
                    ngoodTk++;
                    int nRH_eMipDR=0, nNearTRKs=0;
                    double e1 =  spr::eCone_ecal(geo, barrelRecHitsHandle, endcapRecHitsHandle,
                                                 trkDetItr->pointHCAL, trkDetItr->pointECAL,
                                                 a_neutR1, trkDetItr->directionECAL, nRH_eMipDR);
                    double e2 =  spr::eCone_ecal(geo, barrelRecHitsHandle, endcapRecHitsHandle,
                                                 trkDetItr->pointHCAL, trkDetItr->pointECAL,
                                                 a_neutR2, trkDetItr->directionECAL, nRH_eMipDR);
                    eMipDR = spr::eCone_ecal(geo, barrelRecHitsHandle, endcapRecHitsHandle,
                                             trkDetItr->pointHCAL, trkDetItr->pointECAL,
                                             a_mipR, trkDetItr->directionECAL, nRH_eMipDR);
                    conehmaxNearP = spr::chargeIsolationCone(nTracks, trkCaloDirections, a_charIsoR, nNearTRKs, ((verbosity/100)%10>2));
                    e_inCone = e2 - e1;
                    if (eMipDR<1.0) nselTk++;
                  }
                  if (l3Track) {
                    fillHist(10,v4);
                    if (selectTk) {
                      fillHist(11,v4);
                      fillCuts(0, eMipDR, conehmaxNearP, e_inCone, v4, ieta, (mindR>dr_L1));
                      if (conehmaxNearP < cutCharge) {
                        fillHist(12,v4);
                        if (eMipDR < cutMip) {
                          fillHist(13,v4);
                          if (e_inCone < cutNeutral) fillHist(14, v4);
                        }
                      }
                    }
                  } else {
                    fillHist(15,v4);
                    if (selectTk) {
                      fillHist(16,v4);
                      fillCuts(1, eMipDR, conehmaxNearP, e_inCone, v4, ieta, (mindR>dr_L1));
                      if (conehmaxNearP < cutCharge) {
                        fillHist(17,v4);
                        if (eMipDR < cutMip) {
                          fillHist(18,v4);
                          if (e_inCone < cutNeutral) fillHist(19, v4);
                        }
                      }
                    }
                  }
                }
              }
            }
          }
          break;
        }
      }
      
      h_HLT      -> Fill(hlt);
      h_PreL1    -> Fill(preL1);
      h_PreHLT   -> Fill(preHLT);
      h_Pre      -> Fill(prescale);
      h_PreL1wt  -> Fill(preL1, mybxlumi);
      h_PreHLTwt -> Fill(preHLT, mybxlumi);

      // check if trigger names in (new) config                       
      if (changed) {
        changed = false;
        if ((verbosity/10)%10 > 1) {
          std::cout<<"New trigger menu found !!!" << std::endl;
          const unsigned int n(hltConfig_.size());
          for (unsigned itrig=0; itrig<triggerNames_.size(); itrig++) {
            unsigned int triggerindx = hltConfig_.triggerIndex(triggerNames_[itrig]);
            std::cout << triggerNames_[itrig] << " " << triggerindx << " ";
            if (triggerindx >= n)
              std::cout << "does not exist in the current menu" << std::endl;
            else
              std::cout << "exists" << std::endl;
          }
        }
      }
    }
  }
}

void IsoTrig::beginJob() {
  char hname[100], htit[100];
  std::string levels[20] = {"L1", "L2", "L3", 
                            "Reco", "RecoMatch", "RecoNoMatch", 
                            "L2Match", "L2NoMatch", "L3Match", "L3NoMatch", 
                            "HLTTrk", "HLTGoodTrk", "HLTIsoTrk", "HLTMip", "HLTSelect",
                            "nonHLTTrk", "nonHLTGoodTrk", "nonHLTIsoTrk", "nonHLTMip", "nonHLTSelect"};
  std::string pairs[6] = {"L2L3", "L2L3Match", "L2L3NoMatch", "L3Reco", "L3RecoMatch", "L3RecoNoMatch"};

  std::string cuts[2] = {"HLTMatched", "HLTNotMatched"};
  std::string cuts2[2] = {"All", "Away from L1"};
  h_nHLT        = fs->make<TH1I>("h_nHLT" , "size of rigger Names", 1000, 1, 1000);
  h_HLT         = fs->make<TH1I>("h_HLT"  , "HLT accept", 3, -1, 2);
  h_PreL1       = fs->make<TH1I>("h_PreL1", "L1 Prescale", 500, 0, 500);
  h_PreHLT      = fs->make<TH1I>("h_PreHLT", "HLT Prescale", 50, 0, 50);
  h_Pre         = fs->make<TH1I>("h_Pre", "Prescale", 3000, 0, 3000);
  h_nL3Objs     = fs->make<TH1I>("h_nL3Objs", "Number of L3 objects", 10, 0, 10);

  h_PreL1wt     = fs->make<TH1D>("h_PreL1wt", "Weighted L1 Prescale", 500, 0, 500);
  h_PreHLTwt    = fs->make<TH1D>("h_PreHLTwt", "Weighted HLT Prescale", 500, 0, 100);
  h_L1ObjEnergy = fs->make<TH1D>("h_L1ObjEnergy", "Energy of L1Object", 500, 0.0, 500.0);

  for(int ipair=0; ipair<6; ipair++) {
    sprintf(hname, "h_dEta%s", pairs[ipair].c_str());
    sprintf(htit, "dEta for %s", pairs[ipair].c_str());
    h_dEta[ipair]        = fs->make<TH1D>(hname, htit, 200, -10.0, 10.0);
    h_dEta[ipair]->GetXaxis()->SetTitle("d#eta");

    sprintf(hname, "h_dPhi%s", pairs[ipair].c_str());
    sprintf(htit, "dPhi for %s", pairs[ipair].c_str());
    h_dPhi[ipair]        = fs->make<TH1D>(hname, htit, 140, -7.0, 7.0);
    h_dPhi[ipair]->GetXaxis()->SetTitle("d#phi");

    sprintf(hname, "h_dPt%s", pairs[ipair].c_str());
    sprintf(htit, "dPt for %s objects", pairs[ipair].c_str());
    h_dPt[ipair]         = fs->make<TH1D>(hname, htit, 400, -200.0, 200.0);
    h_dPt[ipair]->GetXaxis()->SetTitle("dp_{T} (GeV)");

    sprintf(hname, "h_dP%s", pairs[ipair].c_str());
    sprintf(htit, "dP for %s objects", pairs[ipair].c_str());
    h_dP[ipair]         = fs->make<TH1D>(hname, htit, 400, -200.0, 200.0);
    h_dP[ipair]->GetXaxis()->SetTitle("dP (GeV)");

    sprintf(hname, "h_dinvPt%s", pairs[ipair].c_str());
    sprintf(htit, "dinvPt for %s objects", pairs[ipair].c_str());
    h_dinvPt[ipair]      = fs->make<TH1D>(hname, htit, 500, -0.4, 0.1);
    h_dinvPt[ipair]->GetXaxis()->SetTitle("d(1/p_{T})");

    sprintf(hname, "h_mindR%s", pairs[ipair].c_str());
    sprintf(htit, "mindR for %s objects", pairs[ipair].c_str());
    h_mindR[ipair]       = fs->make<TH1D>(hname, htit, 500, 0.0, 1.0);
    h_mindR[ipair]->GetXaxis()->SetTitle("dR");
  }

  for(int ilevel=0; ilevel<20; ilevel++) {
    sprintf(hname, "h_p%s", levels[ilevel].c_str());
    sprintf(htit, "p for %s objects", levels[ilevel].c_str());
    h_p[ilevel] = fs->make<TH1D>(hname, htit, 100, 0.0, 500.0);
    h_p[ilevel]->GetXaxis()->SetTitle("p (GeV)");

    sprintf(hname, "h_pt%s", levels[ilevel].c_str());
    sprintf(htit, "pt for %s objects", levels[ilevel].c_str());
    h_pt[ilevel] = fs->make<TH1D>(hname, htit, 100, 0.0, 500.0);
    h_pt[ilevel]->GetXaxis()->SetTitle("p_{T} (GeV)");

    sprintf(hname, "h_eta%s", levels[ilevel].c_str());
    sprintf(htit, "eta for %s objects", levels[ilevel].c_str());
    h_eta[ilevel] = fs->make<TH1D>(hname, htit, 100, -5.0, 5.0);
    h_eta[ilevel]->GetXaxis()->SetTitle("#eta");

    sprintf(hname, "h_phi%s", levels[ilevel].c_str());
    sprintf(htit, "phi for %s objects", levels[ilevel].c_str());
    h_phi[ilevel] = fs->make<TH1D>(hname, htit, 70, -3.5, 3.50);
    h_phi[ilevel]->GetXaxis()->SetTitle("#phi");
  }
  for(int lvl=0; lvl<2; lvl++) {
    sprintf(hname, "h_dEtaL1%s", levels[lvl+1].c_str());
    sprintf(htit, "dEta for L1 and %s objects", levels[lvl+1].c_str());
    h_dEtaL1[lvl] = fs->make<TH1D>(hname, htit, 400, -10.0, 10.0);

    sprintf(hname, "h_dPhiL1%s", levels[lvl+1].c_str());
    sprintf(htit, "dPhi for L1 and %s objects", levels[lvl+1].c_str());
    h_dPhiL1[lvl] = fs->make<TH1D>(hname, htit, 280, -7.0, 7.0);

    sprintf(hname, "h_dRL1%s", levels[lvl+1].c_str());
    sprintf(htit, "dR for L1 and %s objects", levels[lvl+1].c_str());
    h_dRL1[lvl] = fs->make<TH1D>(hname, htit, 100, 0.0, 10.0);
  }

  for(int icut=0; icut<2; icut++) {
    sprintf(hname, "h_eMip%s", cuts[icut].c_str());
    sprintf(htit, "eMip for %s tracks", cuts[icut].c_str());
    h_eMip[icut]     =fs->make<TH1D>(hname, htit, 200, 0.0, 10.0);
    h_eMip[icut]->GetXaxis()->SetTitle("E_{Mip} (GeV)");

    sprintf(hname, "h_eMaxNearP%s", cuts[icut].c_str());
    sprintf(htit, "eMaxNearP for %s tracks", cuts[icut].c_str());
    h_eMaxNearP[icut]=fs->make<TH1D>(hname, htit, 240, -2.0, 10.0);
    h_eMaxNearP[icut]->GetXaxis()->SetTitle("E_{MaxNearP} (GeV)");

    sprintf(hname, "h_eNeutIso%s", cuts[icut].c_str());
    sprintf(htit, "eNeutIso for %s ", cuts[icut].c_str());
    h_eNeutIso[icut] =fs->make<TH1D>(hname, htit, 200, 0.0, 10.0);
    h_eNeutIso[icut]->GetXaxis()->SetTitle("E_{NeutIso} (GeV)");

    for (int kcut=0; kcut<2; ++kcut) {
      sprintf(hname, "h_etaCalibTracks%s%d", cuts[icut].c_str(), kcut);
      sprintf(htit, "etaCalibTracks for %s (%s)", cuts[icut].c_str(), cuts2[kcut].c_str());
      h_etaCalibTracks[icut][kcut]=fs->make<TH1D>(hname, htit, 60, -30.0, 30.0);
      h_etaCalibTracks[icut][kcut]->GetXaxis()->SetTitle("i#eta");

      sprintf(hname, "h_etaMipTracks%s%d", cuts[icut].c_str(), kcut);
      sprintf(htit, "etaMipTracks for %s (%s)", cuts[icut].c_str(), cuts2[kcut].c_str());
      h_etaMipTracks[icut][kcut]=fs->make<TH1D>(hname, htit, 60, -30.0, 30.0);
      h_etaMipTracks[icut][kcut]->GetXaxis()->SetTitle("i#eta");
    }
  }
}
// ------------ method called once each job just after ending the event loop  ------------
void IsoTrig::endJob() {
  unsigned int preL1, preHLT;
  std::map<unsigned int, unsigned int>::iterator itr;
  std::map<unsigned int, const std::pair<int, int>>::iterator itrPre;
  std::cout << "RunNo vs HLT accepts for " << Det << std::endl;
  unsigned int n = maxRunNo - minRunNo +1;
  g_Pre = fs->make<TH1I>("h_PrevsRN", "PreScale Vs Run Number", n, minRunNo, maxRunNo);
  g_PreL1 = fs->make<TH1I>("h_PreL1vsRN", "L1 PreScale Vs Run Number", n, minRunNo, maxRunNo);
  g_PreHLT = fs->make<TH1I>("h_PreHLTvsRN", "HLT PreScale Vs Run Number", n, minRunNo, maxRunNo);
  g_Accepts = fs->make<TH1I>("h_HLTAcceptsvsRN", "HLT Accepts Vs Run Number", n, minRunNo, maxRunNo); 

  for (itr=TrigList.begin(), itrPre=TrigPreList.begin(); itr!=TrigList.end(); itr++, itrPre++) {
    preL1 = (itrPre->second).first;
    preHLT = (itrPre->second).second;
    std::cout << itr->first << " " << itr->second << " " <<  itrPre->first << " " << preL1 << " " << preHLT << std::endl;
    g_Accepts->Fill(itr->first, itr->second);
    g_PreL1->Fill(itr->first, preL1);
    g_PreHLT->Fill(itr->first, preHLT);
    g_Pre->Fill(itr->first, preL1*preHLT);
  }
}

// ------------ method called when starting to processes a run  ------------
void IsoTrig::beginRun(edm::Run const& iRun, edm::EventSetup const& iSetup) {
  std::cout  << "Run " << iRun.run() << " hltconfig.init " 
             << hltConfig_.init(iRun,iSetup,"HLT",changed) << std::endl;;
}

// ------------ method called when ending the processing of a run  ------------
void IsoTrig::endRun(edm::Run const&, edm::EventSetup const&) {}

// ------------ method called when starting to processes a luminosity block  ------------
void IsoTrig::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) {}
// ------------ method called when ending the processing of a luminosity block  ------------
void IsoTrig::endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) {}
// ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
void IsoTrig::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
  //The following says we do not know what parameters are allowed so do no validation
  // Please change this to state exactly what you do use, even if it is no parameters
  edm::ParameterSetDescription desc;
  desc.setUnknown();
  descriptions.addDefault(desc);
}

void IsoTrig::fillHist(int indx, math::XYZTLorentzVector& vec) {
  h_p[indx]->Fill(vec.r());
  h_pt[indx]->Fill(vec.pt());
  h_eta[indx]->Fill(vec.eta());
  h_phi[indx]->Fill(vec.phi());
}

void IsoTrig::fillDifferences(int indx, math::XYZTLorentzVector& vec1, 
                              math::XYZTLorentzVector& vec2, bool debug) {
  double dr     = dR(vec1,vec2);
  double deta   = dEta(vec1, vec2);
  double dphi   = dPhi(vec1, vec2);
  double dpt    = dPt(vec1, vec2);
  double dp     = dP(vec1, vec2);
  double dinvpt = dinvPt(vec1, vec2);
  h_dEta[indx]  ->Fill(deta);
  h_dPhi[indx]  ->Fill(dphi);
  h_dPt[indx]   ->Fill(dpt);
  h_dP[indx]    ->Fill(dp);
  h_dinvPt[indx]->Fill(dinvpt);
  h_mindR[indx] ->Fill(dr);
  if (debug) std::cout << "mindR for index " << indx << " is " << dr << " deta " << deta 
                       << " dphi " << dphi << " dpt " << dpt <<  " dinvpt " << dinvpt <<std::endl;
}

void IsoTrig::fillCuts(int indx, double eMipDR, double conehmaxNearP, double e_inCone, math::XYZTLorentzVector& vec, int ieta, bool cut) {
  h_eMip[indx]     ->Fill(eMipDR);
  h_eMaxNearP[indx]->Fill(conehmaxNearP);
  h_eNeutIso[indx] ->Fill(e_inCone);
  if (conehmaxNearP < cutCharge && eMipDR < cutMip && vec.r()<60 && vec.r()>40) {
    h_etaMipTracks[indx][0]->Fill((double)(ieta));
    if (cut) h_etaMipTracks[indx][1]->Fill((double)(ieta));
    if (e_inCone < cutNeutral) {
      h_etaCalibTracks[indx][0]->Fill((double)(ieta));
      if (cut) h_etaCalibTracks[indx][1]->Fill((double)(ieta));
    }
  }
}

double IsoTrig::dEta(math::XYZTLorentzVector& vec1, math::XYZTLorentzVector& vec2) {
  return (vec1.eta()-vec2.eta());
}

double IsoTrig::dPhi(math::XYZTLorentzVector& vec1, math::XYZTLorentzVector& vec2) {

  double phi1 = vec1.phi();
  if (phi1 < 0) phi1 += 2.0*M_PI;
  double phi2 = vec2.phi();
  if (phi2 < 0) phi2 += 2.0*M_PI;
  double dphi = phi1-phi2;
  if (dphi > M_PI)       dphi -= 2.*M_PI;
  else if (dphi < -M_PI) dphi += 2.*M_PI;
  return dphi;
}

double IsoTrig::dR(math::XYZTLorentzVector& vec1, math::XYZTLorentzVector& vec2) {
  double deta = dEta(vec1,vec2);
  double dphi = dPhi(vec1,vec2);
  return std::sqrt(deta*deta + dphi*dphi);
}

double IsoTrig::dPt(math::XYZTLorentzVector& vec1, math::XYZTLorentzVector& vec2) {
  return (vec1.pt()-vec2.pt());
}

double IsoTrig::dP(math::XYZTLorentzVector& vec1, math::XYZTLorentzVector& vec2) {
  return (std::abs(vec1.r()-vec2.r()));
}

double IsoTrig::dinvPt(math::XYZTLorentzVector& vec1, math::XYZTLorentzVector& vec2) {
  return ((1/vec1.pt())-(1/vec2.pt()));
}
//define this as a plug-in
DEFINE_FWK_MODULE(IsoTrig);