CSCTriggerPrimitivesReader.cc

Go to the documentation of this file.

//-------------------------------------------------
//
//   Class: CSCTriggerPrimitivesReader
//
//   Description: Basic analyzer class which accesses ALCTs, CLCTs, and
//                correlated LCTs and plot various quantities.
//
//   Author List: S. Valuev, UCLA.
//
//
//   Modifications:
//
//--------------------------------------------------

//-----------------------
// This Class's Header --
//-----------------------
#include "L1Trigger/CSCTriggerPrimitives/plugins/CSCTriggerPrimitivesReader.h"

//-------------------------------
// Collaborating Class Headers --
//-------------------------------

#include "FWCore/Framework/interface/ConsumesCollector.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "DataFormats/Common/interface/Handle.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "Geometry/Records/interface/MuonGeometryRecord.h"

#include "CondFormats/CSCObjects/interface/CSCBadChambers.h"
#include "CondFormats/DataRecord/interface/CSCBadChambersRcd.h"

// MC data
#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"

// MC tests
#include "L1Trigger/CSCTriggerPrimitives/plugins/CSCAnodeLCTAnalyzer.h"
#include "L1Trigger/CSCTriggerPrimitives/plugins/CSCCathodeLCTAnalyzer.h"
#include "Geometry/CSCGeometry/interface/CSCLayer.h"
#include "Geometry/CSCGeometry/interface/CSCGeometry.h"

#include "TCanvas.h"
#include "TFile.h"
#include "TH1.h"
#include "TText.h"
#include "TPaveLabel.h"
#include "TPostScript.h"
#include "TStyle.h"
#include "TROOT.h"

#include "FWCore/ServiceRegistry/interface/Service.h"
#include "CommonTools/UtilAlgos/interface/TFileService.h"

#define MAXSTUBS 100
using namespace std;

//-----------------
// Static variables
//-----------------


// Various useful constants
const double CSCTriggerPrimitivesReader::TWOPI = 2.*M_PI;
const string CSCTriggerPrimitivesReader::csc_type[CSC_TYPES] = {
  "ME1/1",  "ME1/2",  "ME1/3",  "ME1/A",  "ME2/1",  "ME2/2",
  "ME3/1",  "ME3/2",  "ME4/1",  "ME4/2"};
const string CSCTriggerPrimitivesReader::csc_type_plus[CSC_TYPES] = {
  "ME+1/1", "ME+1/2", "ME+1/3", "ME+1/A", "ME+2/1", "ME+2/2",
  "ME+3/1", "ME+3/2", "ME+4/1", "ME+4/2"};
const string CSCTriggerPrimitivesReader::csc_type_minus[CSC_TYPES] = {
  "ME-1/1", "ME-1/2", "ME-1/3", "ME-1/A", "ME-2/1", "ME-2/2",
  "ME-3/1", "ME-3/2", "ME-4/1", "ME-4/2"};
const int CSCTriggerPrimitivesReader::NCHAMBERS[CSC_TYPES] = {
  36, 36, 36, 36, 18, 36, 18, 36, 18, 36};
const int CSCTriggerPrimitivesReader::MAX_WG[CSC_TYPES] = {
  48,  64,  32,  48, 112,  64,  96,  64,  96,  64};//max. number of wiregroups
const int CSCTriggerPrimitivesReader::MAX_HS[CSC_TYPES] = {
  128, 160, 128,  96, 160, 160, 160, 160, 160, 160}; // max. # of halfstrips
const int CSCTriggerPrimitivesReader::ptype_TMB07[CSCConstants::NUM_CLCT_PATTERNS]= {
  -999,  -5,  4, -4,  3, -3,  2, -2,  1, -1,  0}; // "signed" pattern (== phiBend)

// LCT counters
int  CSCTriggerPrimitivesReader::numALCT   = 0;
int  CSCTriggerPrimitivesReader::numCLCT   = 0;
int  CSCTriggerPrimitivesReader::numLCTTMB = 0;
int  CSCTriggerPrimitivesReader::numLCTMPC = 0;

bool CSCTriggerPrimitivesReader::bookedHotWireHistos = false;
bool CSCTriggerPrimitivesReader::bookedALCTHistos    = false;
bool CSCTriggerPrimitivesReader::bookedCLCTHistos    = false;
bool CSCTriggerPrimitivesReader::bookedLCTTMBHistos  = false;
bool CSCTriggerPrimitivesReader::bookedLCTMPCHistos  = false;

bool CSCTriggerPrimitivesReader::bookedCompHistos    = false;

bool CSCTriggerPrimitivesReader::bookedResolHistos   = false;
bool CSCTriggerPrimitivesReader::bookedEfficHistos   = false;

bool CSCTriggerPrimitivesReader::printps = false;

// LUT for which ME1/1 wire group can cross which ME1/a halfstrip
// 1st index: WG number
// 2nd index: inclusive HS range
const int CSCTriggerPrimitivesReader::lut_wg_vs_hs_me1a[48][2] = {
  {0, 95},{0, 95},{0, 95},{0, 95},{0, 95},
  {0, 95},{0, 95},{0, 95},{0, 95},{0, 95},
  {0, 95},{0, 95},{0, 77},{0, 61},{0, 39},
  {0, 22},{-1,-1},{-1,-1},{-1,-1},{-1,-1},
  {-1,-1},{-1,-1},{-1,-1},{-1,-1},{-1,-1},
  {-1,-1},{-1,-1},{-1,-1},{-1,-1},{-1,-1},
  {-1,-1},{-1,-1},{-1,-1},{-1,-1},{-1,-1},
  {-1,-1},{-1,-1},{-1,-1},{-1,-1},{-1,-1},
  {-1,-1},{-1,-1},{-1,-1},{-1,-1},{-1,-1},
  {-1,-1},{-1,-1},{-1,-1} };
// a modified LUT for ganged ME1a
const int CSCTriggerPrimitivesReader::lut_wg_vs_hs_me1ag[48][2] = {
  {0, 31},{0, 31},{0, 31},{0, 31},{0, 31},
  {0, 31},{0, 31},{0, 31},{0, 31},{0, 31},
  {0, 31},{0, 31},{0, 31},{0, 31},{0, 31},
  {0, 22},{-1,-1},{-1,-1},{-1,-1},{-1,-1},
  {-1,-1},{-1,-1},{-1,-1},{-1,-1},{-1,-1},
  {-1,-1},{-1,-1},{-1,-1},{-1,-1},{-1,-1},
  {-1,-1},{-1,-1},{-1,-1},{-1,-1},{-1,-1},
  {-1,-1},{-1,-1},{-1,-1},{-1,-1},{-1,-1},
  {-1,-1},{-1,-1},{-1,-1},{-1,-1},{-1,-1},
  {-1,-1},{-1,-1},{-1,-1} };

// LUT for which ME1/1 wire group can cross which ME1/b halfstrip
// 1st index: WG number
// 2nd index: inclusive HS range
const int CSCTriggerPrimitivesReader::lut_wg_vs_hs_me1b[48][2] = {
  {-1,-1},{-1,-1},{-1,-1},{-1,-1},{-1,-1},
  {-1,-1},{-1,-1},{-1,-1},{-1,-1},{-1,-1},
  {100, 127},{73, 127},{47, 127},{22, 127},{0, 127},
  {0, 127},{0, 127},{0, 127},{0, 127},{0, 127},
  {0, 127},{0, 127},{0, 127},{0, 127},{0, 127},
  {0, 127},{0, 127},{0, 127},{0, 127},{0, 127},
  {0, 127},{0, 127},{0, 127},{0, 127},{0, 127},
  {0, 127},{0, 127},{0, 127},{0, 127},{0, 127},
  {0, 127},{0, 127},{0, 127},{0, 127},{0, 105},
  {0, 93},{0, 78},{0, 63} };

//init My ttree
void TreePerStub::init(int run, int event){

  t_EventNumberAnalyzed = -1;
  t_RUN = run;
  t_Event = event;
  t_nStubs = 0;
  t_nStubs_readout = 0;

  t_chamber = -1;
  t_ring = -1;
  t_station = -1;
  t_chambertype = -1;
  t_endcap = -2;
}

TTree *TreePerStub::bookTree(TTree *t, const std::string & name)
{
  edm::Service< TFileService > fs;
  t = fs->make<TTree>(name.c_str(), name.c_str());

  t->Branch("t_EventNumberAnalyzed",&t_EventNumberAnalyzed,"t_EventNumberAnalyzed/I");
  t->Branch("t_RUN",&t_RUN,"t_RUN/I");
  t->Branch("t_Event",&t_Event,"t_Event/I");
  t->Branch("t_nStubs",&t_nStubs,"t_nStubs/I");
  t->Branch("t_nStubs_readout",&t_nStubs_readout,"t_nStubs_readout/I");
  t->Branch("t_nStubs_ME119",&t_nStubs_ME119,"t_nStubs_ME119/I");
  t->Branch("t_nStubs_ME11no911",&t_nStubs_ME11no911,"t_nStubs_ME11no911/I");
  t->Branch("t_nStubs_noMEpm11",&t_nStubs_noMEpm11,"t_nStubs_noMEpm11/I");

  t->Branch("t_chamber", &t_chamber, "t_chamber/I");
  t->Branch("t_ring", &t_ring, "t_ring/I");
  t->Branch("t_endcap", &t_endcap, "t_endcap/I");
  t->Branch("t_station", &t_station, "t_station/I");
  t->Branch("t_chambertype", &t_chambertype, "t_chambertype/I");

  return t;
}


//init stub comparison struct
void MyStubComparison::init(int run, int event){

  nEvents = -1;
  nRUN = run;
  nEvent = event;
  firstfill =  false;
  totStubs_data = -1;
  totStubs_emul = -1;
  totStubs_emul_readout = -1;
  nStub_data = 0;
  nStub_emul = 0;
  chamber = -1;
  ring =-1;
  endcap =-1;
  station = -1;
  chambertype = -1;
  has_data = false;
  has_emul = false;

  npretrig = 0;
  quality_pretrig = -1;
  maxquality_pretrig = -1;
  bend_pretrig = -1;
  bx_pretrig = -1;
  key_hs_pretrig = -1;
  pattern_pretrig = -1;
  maxpattern_pretrig = -1;
  quality_data = -1;
  bend_data = -1;
  bx_data = -1;
  quality_emul = -1;
  bend_emul = -1;
  pattern_data = -1;
  pattern_emul = -1;
  bx_emul = -1;
  bx_corr_emul = -1;//corrected
  key_WG_data = -1;
  key_WG_emul = -1;
  key_hs_data  = -1;
  key_hs_emul = -1;
  WGcrossHS_data = false;
  WGcrossHS_emul = false;
  trknmb_data = -1;
  trknmb_emul = -1;
  dphi_data = -1;
  dphi_emul = -1;
  eta_data = -1;
  eta_emul = -1;
  phi_data = -1;
  phi_emul = -1;
}


TTree *MyStubComparison::bookTree(TTree *t, const std::string & name)
{
  edm::Service< TFileService > fs;
  t = fs->make<TTree>(name.c_str(), name.c_str());

  t->Branch("nEvents",&nEvents);
  t->Branch("nRUN",&nRUN);
  t->Branch("nEvent",&nEvent);
  t->Branch("firstfill",&firstfill);
  t->Branch("totStubs_data",&totStubs_data);
  t->Branch("totStubs_emul",&totStubs_emul);
  t->Branch("totStubs_emul_readout",&totStubs_emul_readout);
  t->Branch("nStub_data",&nStub_data);
  t->Branch("nStub_emul",&nStub_emul);

  t->Branch("chamber",&chamber);
  t->Branch("ring",&ring);
  t->Branch("endcap",&endcap);
  t->Branch("station",&station);
  t->Branch("chambertype",&chambertype);
  t->Branch("has_data",&has_data);
  t->Branch("has_emul",&has_emul);
  t->Branch("quality_data",&quality_data);
  t->Branch("quality_emul",&quality_emul);
  t->Branch("npretrig",&npretrig);
  t->Branch("quality_pretrig",&quality_pretrig);
  t->Branch("maxquality_pretrig",&maxquality_pretrig);
  t->Branch("pattern_data",&pattern_data);
  t->Branch("pattern_emul",&pattern_emul);
  t->Branch("pattern_pretrig",&pattern_pretrig);
  t->Branch("maxpattern_pretrig",&maxpattern_pretrig);
  t->Branch("bend_data",&bend_data);
  t->Branch("bx_data",&bx_data);
  t->Branch("fullbx_data",&fullbx_data);
  t->Branch("bend_emul",&bend_emul);
  t->Branch("bx_emul",&bx_emul);
  t->Branch("fullbx_emul",&fullbx_emul);
  t->Branch("bend_pretrig",&bend_pretrig);
  t->Branch("bx_pretrig",&bx_pretrig);
  t->Branch("bx_corr_emul",&bx_corr_emul);
  t->Branch("WGcrossHS_data",&WGcrossHS_data);
  t->Branch("WGcrossHS_emul",&WGcrossHS_emul);
  t->Branch("key_WG_data",&key_WG_data);
  t->Branch("key_WG_emul",&key_WG_emul);
  t->Branch("key_hs_data",&key_hs_data);
  t->Branch("key_hs_emul",&key_hs_emul);
  t->Branch("key_hs_pretrig",&key_hs_pretrig);
  t->Branch("trknmb_data",&trknmb_data);
  t->Branch("trknmb_emul",&trknmb_emul);
  t->Branch("dphi_data",&dphi_data);
  t->Branch("dphi_emul",&dphi_emul);
  t->Branch("eta_data",&eta_data);
  t->Branch("eta_emul",&eta_emul);
  t->Branch("phi_data", &phi_data);
  t->Branch("phi_emul", &phi_emul);

  return t;

}


//----------------
// Constructor  --
//----------------
CSCTriggerPrimitivesReader::CSCTriggerPrimitivesReader(const edm::ParameterSet& conf) : eventsAnalyzed(0) {
  edm::Service<TFileService> fs;
  //  rootFileName = conf.getUntrackedParameter<string>("rootFileName","TPEHists.root");
  // Create the root file for the histograms
  //  theFile = new TFile(rootFileName.c_str(), "RECREATE");
  //  theFile->cd();

  // Various input parameters.

  printps = conf.getParameter<bool>("printps");
  dataLctsIn_ = conf.getParameter<bool>("dataLctsIn");
  emulLctsIn_ = conf.getParameter<bool>("emulLctsIn");
  edm::ParameterSet commonParams =
    conf.getParameter<edm::ParameterSet>("commonParam");

  // Switch for a new (2007) version of the TMB firmware.
  gangedME1a = commonParams.getParameter<bool>("gangedME1a");

  // is it (non-upgrade algorithm) run along with upgrade one?
  plotME1A = true;
  plotME42 = true;
  lctProducerData_ = conf.getUntrackedParameter<string>("CSCLCTProducerData",
                                                        "cscunpacker");
  mpclctProducerData_ = conf.getUntrackedParameter<string>("CSCMPCLCTProducerData",
                                                           "csctfDigis");
  lctProducerEmul_ = conf.getUntrackedParameter<string>("CSCLCTProducerEmul",
                                                        "cscTriggerPrimitiveDigis");

  simHitProducer_ = conf.getParameter<edm::InputTag>("CSCSimHitProducer");
  wireDigiProducer_ = conf.getParameter<edm::InputTag>("CSCWireDigiProducer");
  compDigiProducer_ = conf.getParameter<edm::InputTag>("CSCComparatorDigiProducer");

  simHit_token_     = consumes<edm::PSimHitContainer>(simHitProducer_);
  wireDigi_token_   = consumes<CSCWireDigiCollection>(wireDigiProducer_);
  compDigi_token_   = consumes<CSCComparatorDigiCollection>(compDigiProducer_);

  alcts_d_token_    = consumes<CSCALCTDigiCollection>(edm::InputTag(lctProducerData_, "MuonCSCALCTDigi"));
  clcts_d_token_    = consumes<CSCCLCTDigiCollection>(edm::InputTag(lctProducerData_, "MuonCSCCLCTDigi"));
  lcts_tmb_d_token_ = consumes<CSCCorrelatedLCTDigiCollection>(edm::InputTag(lctProducerData_, "MuonCSCCorrelatedLCTDigi"));
  lcts_mpc_d_token_ = consumes<CSCCorrelatedLCTDigiCollection>(edm::InputTag(mpclctProducerData_));

  alcts_e_token_    = consumes<CSCALCTDigiCollection>(edm::InputTag(lctProducerEmul_));
  clcts_e_token_    = consumes<CSCCLCTDigiCollection>(edm::InputTag(lctProducerEmul_));
  pretrigs_e_token_ = consumes<CSCCLCTPreTriggerDigiCollection>(edm::InputTag(lctProducerEmul_));
  lcts_tmb_e_token_ = consumes<CSCCorrelatedLCTDigiCollection>(edm::InputTag(lctProducerEmul_));
  lcts_mpc_e_token_ = consumes<CSCCorrelatedLCTDigiCollection>(edm::InputTag(lctProducerEmul_, "MPCSORTED"));

  consumesMany<edm::HepMCProduct>();
  resultsFileNamesPrefix_ = conf.getUntrackedParameter<string>("resultsFileNamesPrefix","");
  checkBadChambers_ = conf.getUntrackedParameter<bool>("checkBadChambers",true);
  debug = conf.getUntrackedParameter<bool>("debug", false);
  dataIsAnotherMC_ = conf.getUntrackedParameter<bool>("dataIsAnotherMC", false);

  //rootFileName = conf.getUntrackedParameter<string>("rootFileName");
  // Create the root file.
  // Not sure we really need it - comment out for now. -Slava.
  //theFile = new TFile(rootFileName.c_str(), "RECREATE");
  //theFile->cd();
  stub_tree[0] = stubs_comparison[0].bookTree(stub_tree[0],"alcttree");
  stub_tree[1] = stubs_comparison[1].bookTree(stub_tree[1],"clcttree");
  stub_tree[2] = stubs_comparison[2].bookTree(stub_tree[2],"lcttree");
  stub_tree[3] = stubs_comparison[3].bookTree(stub_tree[3],"mpclcttree");
  // Per event TTree
  event_tree[0] = perStub[0].bookTree(event_tree[0],"Ev_alcttree");
  event_tree[1] = perStub[1].bookTree(event_tree[1],"Ev_emul_alcttree");
  event_tree[2] = perStub[2].bookTree(event_tree[2],"Ev_clcttree");
  event_tree[3] = perStub[3].bookTree(event_tree[3],"Ev_emul_clcttree");
  event_tree[4] = perStub[4].bookTree(event_tree[4],"Ev_lcttree");
  event_tree[5] = perStub[5].bookTree(event_tree[5],"Ev_emul_lcttree");
  // My favourite ROOT settings.
  setRootStyle();
}

//----------------
// Destructor   --
//----------------
CSCTriggerPrimitivesReader::~CSCTriggerPrimitivesReader() {
  //  histos->writeHists(theFile);
  //  theFile->Close();
  //delete theFile;
}


int CSCTriggerPrimitivesReader::maxRing(int station)
{
  if (station == 1) {
    if (plotME1A) return 4;
    else  return 3;
  }
  return 2;
}


void CSCTriggerPrimitivesReader::analyze(const edm::Event& ev,
                                         const edm::EventSetup& setup) {
  ++eventsAnalyzed;
  //if (ev.id().event()%10 == 0)
  LogTrace("CSCTriggerPrimitivesReader")
    << "\n** CSCTriggerPrimitivesReader: processing run #"
    << ev.id().run() << " event #" << ev.id().event()
    << "; events so far: " << eventsAnalyzed << " **";
  RUN_ = ev.id().run();
  Event_ = ev.id().event();

  // Find the geometry for this event & cache it.  Needed in LCTAnalyzer
  // modules.
  edm::ESHandle<CSCGeometry> cscGeom;
  setup.get<MuonGeometryRecord>().get(cscGeom);
  geom_ = &*cscGeom;

  // Find conditions data for bad chambers & cache it.  Needed for efficiency
  // calculations.
  if (checkBadChambers_) {
    edm::ESHandle<CSCBadChambers> pBad;
    setup.get<CSCBadChambersRcd>().get(pBad);
    badChambers_ = pBad.product();
  }

  // Get the collections of ALCTs, CLCTs, and correlated LCTs from event.
  edm::Handle<CSCALCTDigiCollection> alcts_data;
  edm::Handle<CSCCLCTDigiCollection> clcts_data;
  edm::Handle<CSCCorrelatedLCTDigiCollection> lcts_tmb_data;
  edm::Handle<CSCCorrelatedLCTDigiCollection> lcts_mpc_data;
  edm::Handle<CSCALCTDigiCollection> alcts_emul;
  edm::Handle<CSCCLCTDigiCollection> clcts_emul;
  edm::Handle<CSCCLCTPreTriggerDigiCollection> pretrigs_emul;
  edm::Handle<CSCCorrelatedLCTDigiCollection> lcts_tmb_emul;
  edm::Handle<CSCCorrelatedLCTDigiCollection> lcts_mpc_emul;

  // Data
  if (dataLctsIn_) {
    HotWires(ev);
    //    ev.getByLabel(lctProducerData_,  alcts_data);
    //    ev.getByLabel(lctProducerData_,  clcts_data);
    //    ev.getByLabel(lctProducerData_,  lcts_tmb_data);
    ev.getByToken(alcts_d_token_, alcts_data);
    ev.getByToken(clcts_d_token_, clcts_data);
    ev.getByToken(lcts_tmb_d_token_, lcts_tmb_data);
    ev.getByToken(lcts_mpc_d_token_, lcts_mpc_data);

    if (!alcts_data.isValid()) {
      edm::LogWarning("L1CSCTPEmulatorWrongInput")
        << "+++ Warning: Collection of ALCTs with label MuonCSCALCTDigi"
        << " requested, but not found in the event... Skipping the rest +++\n";
      return;
    }
    if (!clcts_data.isValid()) {
      edm::LogWarning("L1CSCTPEmulatorWrongInput")
        << "+++ Warning: Collection of CLCTs with label MuonCSCCLCTDigi"
        << " requested, but not found in the event... Skipping the rest +++\n";
      return;
    }
    if (!lcts_tmb_data.isValid()) {
      edm::LogWarning("L1CSCTPEmulatorWrongInput")
        << "+++ Warning: Collection of correlated LCTs with label"
        << " MuonCSCCorrelatedLCTDigi requested, but not found in the"
        << " event... Skipping the rest +++\n";
      return;
    }
    /*
      if (!lcts_mpc_data.isValid()) {
      edm::LogWarning("L1CSCTPEmulatorWrongInput")
      << "+++ Warning: Collection of MPC correlated LCTs with label"
      << " MuonCSCCorrelatedLCTDigi + MCPSorted requested, but not found in the"
      << " event... Skipping the rest +++\n";
      //return;
      }*/
  }

  // Emulator
  if (emulLctsIn_) {
    //    ev.getByLabel(lctProducerEmul_,              alcts_emul);
    //    ev.getByLabel(lctProducerEmul_,              clcts_emul);
    //    ev.getByLabel(lctProducerEmul_,              lcts_tmb_emul);
    //    ev.getByLabel(lctProducerEmul_, "MPCSORTED", lcts_mpc_emul);
    ev.getByToken(alcts_e_token_, alcts_emul);
    ev.getByToken(clcts_e_token_, clcts_emul);
    ev.getByToken(pretrigs_e_token_, pretrigs_emul);
    ev.getByToken(lcts_tmb_e_token_, lcts_tmb_emul);
    ev.getByToken(lcts_mpc_e_token_, lcts_mpc_emul);

    if (!alcts_emul.isValid()) {
      edm::LogWarning("L1CSCTPEmulatorWrongInput")
        << "+++ Warning: Collection of emulated ALCTs"
        << " requested, but not found in the event... Skipping the rest +++\n";
      return;
    }
    if (!clcts_emul.isValid()) {
      edm::LogWarning("L1CSCTPEmulatorWrongInput")
        << "+++ Warning: Collection of emulated CLCTs"
        << " requested, but not found in the event... Skipping the rest +++\n";
      return;
    }
    if (!lcts_tmb_emul.isValid()) {
      edm::LogWarning("L1CSCTPEmulatorWrongInput")
        << "+++ Warning: Collection of emulated correlated LCTs"
        << " requested, but not found in the event... Skipping the rest +++\n";
      return;
    }
    /*
      if (!lcts_mpc_emul.isValid()) {
      edm::LogWarning("L1CSCTPEmulatorWrongInput")
      << "+++ Warning: Collection of emulated correlated LCTs (MPCs)"
      << " requested, but not found in the event... Skipping the rest +++\n";
      //return;
      }*/
  }

  // Fill histograms with reconstructed or emulated quantities.  If both are
  // present, plot LCTs in data.
  if (dataLctsIn_) {
    fillALCTHistos(alcts_data.product());
    fillCLCTHistos(clcts_data.product());
    fillLCTTMBHistos(lcts_tmb_data.product());
  }
  else if (emulLctsIn_) {
    fillALCTHistos(alcts_emul.product());
    fillCLCTHistos(clcts_emul.product());
    fillLCTTMBHistos(lcts_tmb_emul.product());
    //fillLCTMPCHistos(lcts_mpc_emul.product());
  }

  // Compare LCTs in the data with the ones produced by the emulator.
  if (dataLctsIn_ && emulLctsIn_) {
    compare(alcts_data.product(), alcts_emul.product(),
            clcts_data.product(), clcts_emul.product(),
            pretrigs_emul.product(), lcts_tmb_data.product(), lcts_tmb_emul.product()
            );
  }
  // Fill MC-based resolution/efficiency histograms, if needed.
  if (emulLctsIn_) {
    MCStudies(ev, alcts_emul.product(), clcts_emul.product());
  }
} // analyze

void CSCTriggerPrimitivesReader::endJob() {
  // Note: all operations involving ROOT should be placed here and not in the
  // destructor.
  // Plot histos if they were booked/filled.
  if(printps) {
    if (bookedALCTHistos)   drawALCTHistos();
    if (bookedCLCTHistos)   drawCLCTHistos();
    if (bookedLCTTMBHistos) drawLCTTMBHistos();
    //if (bookedLCTMPCHistos) drawLCTMPCHistos();

    if (bookedCompHistos)   drawCompHistos();

    if (bookedResolHistos)  drawResolHistos();
    if (bookedEfficHistos)  drawEfficHistos();
  }
  //drawHistosForTalks();

  //theFile->cd();
  //theFile->Write();
  //theFile->Close();

  // Job summary.
  edm::LogInfo("CSCTriggerPrimitivesReader")
    << "\n  Average number of ALCTs/event = "
    << static_cast<float>(numALCT)/eventsAnalyzed << endl;
  edm::LogInfo("CSCTriggerPrimitivesReader")
    << "  Average number of CLCTs/event = "
    << static_cast<float>(numCLCT)/eventsAnalyzed << endl;
  edm::LogInfo("CSCTriggerPrimitivesReader")
    << "  Average number of TMB LCTs/event = "
    << static_cast<float>(numLCTTMB)/eventsAnalyzed << endl;
  edm::LogInfo("CSCTriggerPrimitivesReader")
    << "  Average number of MPC LCTs/event = "
    << static_cast<float>(numLCTMPC)/eventsAnalyzed << endl;

  if (bookedEfficHistos) {
    {
      edm::LogInfo("CSCTriggerPrimitivesReader") << "\n  ALCT efficiencies:";
      double tot_simh = 0.0, tot_alct = 0.0;
      for (int idh = 0; idh < CSC_TYPES; idh++) {
        double simh = hEfficHitsEtaCsc[idh]->Integral();
        double alct = hEfficALCTEtaCsc[idh]->Integral();
        double eff  = 0.;
        if (simh > 0) eff = alct/simh;
        edm::LogInfo("CSCTriggerPrimitivesReader")
          << "    " << csc_type[idh]
          << ": alct = " << alct << ", simh = " << simh << " eff = " << eff;
        tot_simh += simh;
        tot_alct += alct;
      }
      edm::LogInfo("CSCTriggerPrimitivesReader")
        << "    overall: alct = " << tot_alct << ", simh = " << tot_simh
        << " eff = " << tot_alct/tot_simh;
    }

    {
      edm::LogInfo("CSCTriggerPrimitivesReader") << "\n  CLCT efficiencies:";
      double tot_simh = 0.0, tot_clct = 0.0;
      for (int idh = 0; idh < CSC_TYPES; idh++) {
        double simh = hEfficHitsEtaCsc[idh]->Integral();
        double clct = hEfficCLCTEtaCsc[idh]->Integral();
        double eff  = 0.;
        if (simh > 0.) eff = clct/simh;
        edm::LogInfo("CSCTriggerPrimitivesReader")
          << "    " << csc_type[idh]
          << ": clct = " << clct << ", simh = " << simh << " eff = " << eff;
        tot_simh += simh;
        tot_clct += clct;
      }
      edm::LogInfo("CSCTriggerPrimitivesReader")
        << "    overall: clct = " << tot_clct << ", simh = " << tot_simh
        << " eff = " << tot_clct/tot_simh;
    }
  }

  if (bookedResolHistos) {
    double cor = 0.0, tot = 0.0;
    cor = hResolDeltaHS->GetBinContent(hResolDeltaHS->FindBin(0.));
    tot = hResolDeltaHS->GetEntries();
    edm::LogInfo("CSCTriggerPrimitivesReader")
      << "\n  Correct half-strip assigned in " << cor << "/" << tot
      << " = " << cor/tot << " of half-strip CLCTs";
    //cor = hResolDeltaDS->GetBinContent(hResolDeltaDS->FindBin(0.));
    //tot = hResolDeltaDS->GetEntries();
    //edm::LogInfo("CSCTriggerPrimitivesReader")
    //  << "  Correct di-strip assigned in " << cor << "/" << tot
    //  << " = " << cor/tot << " of di-strip CLCTs";
    cor = hResolDeltaWG->GetBinContent(hResolDeltaWG->FindBin(0.));
    tot = hResolDeltaWG->GetEntries();
    edm::LogInfo("CSCTriggerPrimitivesReader")
      << "  Correct wire group assigned in " << cor << "/" << tot
      << " = " << cor/tot << " of ALCTs";
  }
}

//---------------
// ROOT settings
//---------------
void CSCTriggerPrimitivesReader::setRootStyle() {
  TH1::AddDirectory(false);

  gROOT->SetStyle("Plain");
  gStyle->SetFillColor(0);
  gStyle->SetOptDate();
  gStyle->SetOptStat(111110);
  gStyle->SetOptFit(1111);
  gStyle->SetPadTickX(1);
  gStyle->SetPadTickY(1);
  gStyle->SetMarkerSize(0.5);
  gStyle->SetMarkerStyle(8);
  gStyle->SetGridStyle(3);
  gStyle->SetPaperSize(TStyle::kA4);
  gStyle->SetStatW(0.25); // width of statistics box; default is 0.19
  gStyle->SetStatH(0.10); // height of statistics box; default is 0.1
  gStyle->SetStatFormat("6.4g");  // leave default format for now
  gStyle->SetTitleSize(0.055, "");   // size for pad title; default is 0.02
  // Really big; useful for talks.
  //gStyle->SetTitleSize(0.1, "");   // size for pad title; default is 0.02
  gStyle->SetLabelSize(0.05, "XYZ"); // size for axis labels; default is 0.04
  gStyle->SetStatFontSize(0.06);     // size for stat. box
  gStyle->SetTitleFont(32, "XYZ"); // times-bold-italic font (p. 153) for axes
  gStyle->SetTitleFont(32, "");    // same for pad title
  gStyle->SetLabelFont(32, "XYZ"); // same for axis labels
  gStyle->SetStatFont(32);         // same for stat. box
  gStyle->SetLabelOffset(0.006, "Y"); // default is 0.005
}

//---------------------
// Histograms for LCTs
//---------------------
void CSCTriggerPrimitivesReader::bookHotWireHistos() {
  edm::Service<TFileService> fs;
  hHotWire1  = fs->make<TH1F>("hHotWire1", "hHotWire1",570*6*112,0,570*6*112);
  hHotCham1  = fs->make<TH1F>("hHotCham1", "hHotCham1",570,0,570);
  bookedHotWireHistos = true;
}

void CSCTriggerPrimitivesReader::bookALCTHistos() {
  string s;

  edm::Service<TFileService> fs;
  hAlctPerEvent  = fs->make<TH1F>("ALCTs_per_event", "ALCTs per event",     31, -0.5,  30.5);
  hAlctPerChamber= fs->make<TH1F>("ALCTs_per_chamber", "ALCTs per chamber",    4, -0.5,   3.5);
  hAlctPerCSC    = fs->make<TH1F>("ALCTs_per_CSCtype", "ALCTs per CSC type",  10, -0.5,   9.5);
  for (int i = 0; i < MAX_ENDCAPS; i++) { // endcaps
    for (int j = 0; j < CSC_TYPES; j++) { // station/ring
      float csc_max = NCHAMBERS[j] + 0.5;
      char asdf[256];
      sprintf(asdf,"ALCTs_%i",i*CSC_TYPES+j);
      if (i == 0) s = "ALCTs, " + csc_type_plus[j];
      else        s = "ALCTs, " + csc_type_minus[j];
      hAlctCsc[i][j] = fs->make<TH1F>(asdf, s.c_str(), NCHAMBERS[j], 0.5, csc_max);
    }
  }

  hAlctValid     = fs->make<TH1F>("ALCT_validity", "ALCT validity",        3, -0.5,   2.5);
  hAlctQuality   = fs->make<TH1F>("ALCT_quality", "ALCT quality",         5, -0.5,   4.5);
  hAlctAccel     = fs->make<TH1F>("ALCT_accel_flag", "ALCT accel. flag",     3, -0.5,   2.5);
  hAlctCollis    = fs->make<TH1F>("ALCT_collision_flag", "ALCT collision. flag", 3, -0.5,   2.5);
  hAlctKeyGroup  = fs->make<TH1F>("ALCT_key_wiregroup", "ALCT key wiregroup", 120, -0.5, 119.5);
  hAlctBXN       = fs->make<TH1F>("ALCT_bx", "ALCT bx",             20, -0.5,  19.5);

  hAlctKeyGroupME11 = fs->make<TH1F>("hAlctKeyGroupME11", "ALCT key wiregroup ME1/1", 50, -0.5, 49.5);

  bookedALCTHistos = true;
  std::cout <<" bookedALCTHistos " << std::endl;
}

void CSCTriggerPrimitivesReader::bookCLCTHistos() {
  string s;

  edm::Service<TFileService> fs;
  hClctPerEvent  = fs->make<TH1F>("CLCTs_per_event", "CLCTs per event",    31, -0.5, 30.5);
  hClctPerChamber= fs->make<TH1F>("CLCTs_per_chamber", "CLCTs per chamber",   3, -0.5,  2.5);
  hClctPerCSC    = fs->make<TH1F>("CLCTs_per_CSCtype", "CLCTs per CSC type", 10, -0.5,  9.5);
  for (int i = 0; i < MAX_ENDCAPS; i++) { // endcaps
    for (int j = 0; j < CSC_TYPES; j++) { // station/ring
      char asdf[256];
      sprintf(asdf,"CLCTs_%i",i*CSC_TYPES+j);
      float csc_max = NCHAMBERS[j] + 0.5;
      if (i == 0) s = "CLCTs, " + csc_type_plus[j];
      else        s = "CLCTs, " + csc_type_minus[j];
      hClctCsc[i][j] = fs->make<TH1F>(asdf, s.c_str(), NCHAMBERS[j], 0.5, csc_max);
    }
  }

  hClctValid     = fs->make<TH1F>("CLCT_validity", "CLCT validity",       3, -0.5,  2.5);
  hClctQuality   = fs->make<TH1F>("CLCT_layers_hit", "CLCT layers hit",     9, -0.5,  8.5);
  hClctStripType = fs->make<TH1F>("CLCT_strip_type", "CLCT strip type",     3, -0.5,  2.5);
  hClctSign      = fs->make<TH1F>("CLCT_sing_(L/R)", "CLCT sign (L/R)",     3, -0.5,  2.5);
  hClctCFEB      = fs->make<TH1F>("CLCT_cfeb_#", "CLCT cfeb #",         6, -0.5,  5.5);
  hClctBXN       = fs->make<TH1F>("CLCT_bx", "CLCT bx",            20, -0.5, 19.5);

  hClctKeyStrip[0] = fs->make<TH1F>("CLCT_keystrip_distrips","CLCT keystrip, distrips",   40, -0.5,  39.5);
  //hClctKeyStrip[0] = fs->make<TH1F>("","CLCT keystrip, distrips",  160, -0.5, 159.5);
  hClctKeyStrip[1] = fs->make<TH1F>("CLCT_keystrip_halfstrips","CLCT keystrip, halfstrips",160, -0.5, 159.5);
  hClctPattern[0]  = fs->make<TH1F>("CLCT_pattern_distrips","CLCT pattern, distrips",    13, -0.5,  12.5);
  hClctPattern[1]  = fs->make<TH1F>("CLCT_pattern_halfstrips","CLCT pattern, halfstrips",  13, -0.5,  12.5);

  for (int i = 0; i < CSC_TYPES; i++) {
    char asdf[256];
    string s1 = "CLCT bend, " + csc_type[i];
    sprintf(asdf,"CLCT_bend0_%i",i+1);
    hClctBendCsc[i][0] = fs->make<TH1F>(asdf, s1.c_str(),  5, -0.5, 4.5);
    sprintf(asdf,"CLCT_bend1_%i",i+1);
    hClctBendCsc[i][1] = fs->make<TH1F>(asdf, s1.c_str(),  5, -0.5, 4.5);

    sprintf(asdf,"CLCT_keystrip_%i",i+1);
    string s2 = "CLCT keystrip, " + csc_type[i];
    int max_ds = MAX_HS[i]/4;
    hClctKeyStripCsc[i]   = fs->make<TH1F>(asdf, s2.c_str(), max_ds, 0., max_ds);
  }

  hClctKeyStripME11 = fs->make<TH1F>("hClctKeyStripME11","CLCT keystrip, halfstrips ME1/1",161, -0.5, 160.5);

  bookedCLCTHistos = true;
  std::cout <<" bookedCLCTHistos " << std::endl;
}

void CSCTriggerPrimitivesReader::bookLCTTMBHistos() {
  string s;

  edm::Service<TFileService> fs;
  hLctTMBPerEvent  = fs->make<TH1F>("LCTs_per_event", "LCTs per event",    31, -0.5, 30.5);
  hLctTMBPerChamber= fs->make<TH1F>("LCTs_per_chamber", "LCTs per chamber",   3, -0.5,  2.5);
  hLctTMBPerCSC    = fs->make<TH1F>("LCTs_per_CSCtype", "LCTs per CSC type", 10, -0.5,  9.5);
  hCorrLctTMBPerCSC= fs->make<TH1F>("CorrLCTs_per_CSCtype", "Corr. LCTs per CSC type", 10, -0.5, 9.5);
  hLctTMBEndcap    = fs->make<TH1F>("LCTS_endcap", "Endcap",             4, -0.5,  3.5);
  hLctTMBStation   = fs->make<TH1F>("LCTS_station", "Station",            6, -0.5,  5.5);
  hLctTMBSector    = fs->make<TH1F>("LCTS_sector", "Sector",             8, -0.5,  7.5);
  hLctTMBRing      = fs->make<TH1F>("LCTS_ring", "Ring",               5, -0.5,  4.5);
  for (int i = 0; i < MAX_ENDCAPS; i++) { // endcaps
    for (int j = 0; j < CSC_TYPES; j++) { // station/ring
      char asdf[256];
      float csc_max = NCHAMBERS[j] + 0.5;
      if (i == 0) s = "LCTs, " + csc_type_plus[j];
      else        s = "LCTs, " + csc_type_minus[j];
      sprintf(asdf,"LCTs_%i",i*CSC_TYPES+j);
      hLctTMBCsc[i][j] = fs->make<TH1F>(s.c_str(), s.c_str(), NCHAMBERS[j], 0.5, csc_max);
    }
  }

  hLctTMBValid     = fs->make<TH1F>("LCT_validity", "LCT validity",        3, -0.5,   2.5);
  hLctTMBQuality   = fs->make<TH1F>("LCT_quality", "LCT quality",        17, -0.5,  16.5);
  hLctTMBKeyGroup  = fs->make<TH1F>("LCT_key_wiregroup", "LCT key wiregroup", 120, -0.5, 119.5);
  hLctTMBKeyStrip  = fs->make<TH1F>("LCT_key_strip", "LCT key strip",     160, -0.5, 159.5);
  hLctTMBStripType = fs->make<TH1F>("LCT_strip_type", "LCT strip type",      3, -0.5,   2.5);
  hLctTMBPattern   = fs->make<TH1F>("LCT_pattern", "LCT pattern",        13, -0.5,  12.5);
  hLctTMBBend      = fs->make<TH1F>("LCT_bend", "LCT L/R bend",        3, -0.5,   2.5);
  hLctTMBBXN       = fs->make<TH1F>("LCT_bx", "LCT bx",             20, -0.5,  19.5);

  // LCT quantities per station
  char histname[60];
  for (int istat = 0; istat < MAX_STATIONS; istat++) {
    sprintf(histname, "LCT_CSCId, station %d", istat+1);
    hLctTMBChamber[istat] = fs->make<TH1F>("", histname,  10, -0.5, 9.5);
  }

  hLctTMBKeyGroupME11  = fs->make<TH1F>("hLctTMBKeyGroupME11", "LCT key wiregroup ME1/1", 50, -0.5, 49.5);
  hLctTMBKeyStripME11  = fs->make<TH1F>("hLctTMBKeyStripME11", "LCT key strip ME1/1",     161, -0.5, 160.5);

  bookedLCTTMBHistos = true;
  std::cout <<" bookedLCTTMBHistos " << std::endl;
}

int CSCTriggerPrimitivesReader::chamberIXi(CSCDetId id) {
  //    1/1 1/2 1/3 2/1 2/2 3/1 3/2 4/1 4/2 -1/1 -1/2 -1/3 -2/1 -2/2 -3/1 -3/2 -4/1 -4/2
  //ix= 0   1   2   3   4   5   6   7   8    9    10   11   12   13   14   15   16   17
  int ix=0;
  if(id.station()!=1) {

    ix=(id.station()-2)*2+3;
  }
  ix+=id.ring()-1;
  if(id.endcap()==2) {
    ix+=9;
  }
  return ix;
}

int CSCTriggerPrimitivesReader::chamberIX(CSCDetId id) {
  int ix=1;
  if(id.station()!=1) {
    ix=(id.station()-2)*2+3+1;
  }
  ix+=id.ring()-1;
  if(id.endcap()==2) {
    ix*=-1;
  }
  return ix;
}

int CSCTriggerPrimitivesReader::chamberSerial( CSCDetId id ) {
  int st = id.station();
  int ri = id.ring();
  int ch = id.chamber();
  int ec = id.endcap();
  int kSerial = ch;
  if (st == 1 && ri == 1) kSerial = ch;
  if (st == 1 && ri == 2) kSerial = ch + 36;
  if (st == 1 && ri == 3) kSerial = ch + 72;
  if (st == 1 && ri == 4) kSerial = ch;
  if (st == 2 && ri == 1) kSerial = ch + 108;
  if (st == 2 && ri == 2) kSerial = ch + 126;
  if (st == 3 && ri == 1) kSerial = ch + 162;
  if (st == 3 && ri == 2) kSerial = ch + 180;
  if (st == 4 && ri == 1) kSerial = ch + 216;
  if (st == 4 && ri == 2) kSerial = ch + 234;  // one day...
  if (ec == 2) kSerial = kSerial + 300;
  return kSerial;
}


void CSCTriggerPrimitivesReader::bookLCTMPCHistos() {
  edm::Service<TFileService> fs;
  hLctMPCPerEvent  = fs->make<TH1F>("MPC_per_event", "LCTs per event",    31, -0.5, 30.5);
  hLctMPCPerCSC    = fs->make<TH1F>("MPC_per_CSCtype", "LCTs per CSC type", 10, -0.5,  9.5);
  hCorrLctMPCPerCSC= fs->make<TH1F>("CorrMPC_per_CSCtype", "Corr. LCTs per CSC type", 10, -0.5,9.5);
  hLctMPCEndcap    = fs->make<TH1F>("MPC_Endcap", "Endcap",             4, -0.5,  3.5);
  hLctMPCStation   = fs->make<TH1F>("MPC_Station", "Station",            6, -0.5,  5.5);
  hLctMPCSector    = fs->make<TH1F>("MPC_Sector", "Sector",             8, -0.5,  7.5);
  hLctMPCRing      = fs->make<TH1F>("MPC_Ring", "Ring",               5, -0.5,  4.5);

  hLctMPCValid     = fs->make<TH1F>("MPC_validity", "LCT validity",        3, -0.5,   2.5);
  hLctMPCQuality   = fs->make<TH1F>("MPC_quality", "LCT quality",        17, -0.5,  16.5);
  hLctMPCKeyGroup  = fs->make<TH1F>("MPC_key_wiregroup", "LCT key wiregroup", 120, -0.5, 119.5);
  hLctMPCKeyStrip  = fs->make<TH1F>("MPC_key_strip", "LCT key strip",     160, -0.5, 159.5);
  hLctMPCStripType = fs->make<TH1F>("MPC_strip_type", "LCT strip type",      3, -0.5,   2.5);
  hLctMPCPattern   = fs->make<TH1F>("MPC_pattern", "LCT pattern",        13, -0.5,  12.5);
  hLctMPCBend      = fs->make<TH1F>("MPC_bend", "LCT L/R bend",        3, -0.5,   2.5);
  hLctMPCBXN       = fs->make<TH1F>("MPC_bx", "LCT bx",             20, -0.5,  19.5);

  // LCT quantities per station
  char histname[60];
  for (int istat = 0; istat < MAX_STATIONS; istat++) {
    sprintf(histname, "MPC_CSCId, station %d", istat+1);
    hLctMPCChamber[istat] = fs->make<TH1F>("", histname,  10, -0.5, 9.5);
  }

  hLctMPCKeyGroupME11  = fs->make<TH1F>("hLctMPCKeyGroupME11", "MPC LCT key wiregroup ME1/1", 50, -0.5, 49.5);
  hLctMPCKeyStripME11  = fs->make<TH1F>("hLctMPCKeyStripME11", "MPC LCT key strip ME1/1",     161, -0.5, 160.5);

  bookedLCTMPCHistos = true;
  if (debug) LogTrace("CSCTriggerPrimitivesReader")
               <<"bookedLCTMPCHistos  ";
  std::cout <<" bookedLCTMPCHistos " << std::endl;

}

void CSCTriggerPrimitivesReader::bookCompHistos() {
  string s;



  edm::Service<TFileService> fs;
  //  hAlctCompMatch;

  hAlctCompFound = fs->make<TH1F>("h_ALCT_found","h_ALCT_found",600,0.5,600.5);
  hAlctCompFound2 = fs->make<TH2F>("h_ALCT_found2","h_ALCT_found2",19,-9.5,9.5,36,0.5,36.5);
  hAlctCompFound2x = fs->make<TH2F>("h_ALCT_found2x","h_ALCT_found2x",19,-9.5,9.5,18,0.5,36.5);


  hAlctCompSameN = fs->make<TH1F>("h_ALCT_SameN","h_ALCT_SameN",600,0.5,600.5);
  hAlctCompSameN2 = fs->make<TH2F>("h_ALCT_SameN2","h_ALCT_SameN2",19,-9.5,9.5,36,0.5,36.5);
  hAlctCompSameN2x = fs->make<TH2F>("h_ALCT_SameN2x","h_ALCT_SameN2x",19,-9.5,9.5,18,0.5,36.5);


  hAlctCompMatch = fs->make<TH1F>("h_ALCT_match","h_ALCT_match",600,0.5,600.5);
  hAlctCompMatch2 = fs->make<TH2F>("h_ALCT_match2","h_ALCT_match2",19,-9.5,9.5,36,0.5,36.5);
  hAlctCompMatch2x = fs->make<TH2F>("h_ALCT_match2x","h_ALCT_match2x",19,-9.5,9.5,18,0.5,36.5);

  hAlctCompTotal = fs->make<TH1F>("h_ALCT_total","h_ALCT_total",600,0.5,600.5);
  hAlctCompTotal2 = fs->make<TH2F>("h_ALCT_total2","h_ALCT_total2",19,-9.5,9.5,36,0.5,36.5);
  hAlctCompTotal2x = fs->make<TH2F>("h_ALCT_total2x","h_ALCT_total2x",19,-9.5,9.5,18,0.5,36.5);

  hClctCompFound = fs->make<TH1F>("h_CLCT_found","h_CLCT_found",600,0.5,600.5);
  hClctCompFound2 = fs->make<TH2F>("h_CLCT_found2","h_CLCT_found2",19,-9.5,9.5,36,0.5,36.5);
  hClctCompFound2x = fs->make<TH2F>("h_CLCT_found2x","h_CLCT_found2x",19,-9.5,9.5,18,0.5,36.5);

  hClctCompSameN = fs->make<TH1F>("h_CLCT_SameN","h_CLCT_SameN",600,0.5,600.5);
  hClctCompSameN2 = fs->make<TH2F>("h_CLCT_SameN2","h_CLCT_SameN2",19,-9.5,9.5,36,0.5,36.5);
  hClctCompSameN2x = fs->make<TH2F>("h_CLCT_SameN2x","h_CLCT_SameN2x",19,-9.5,9.5,18,0.5,36.5);

  hClctCompMatch = fs->make<TH1F>("h_CLCT_match","h_CLCT_match",600,0.5,600.5);
  hClctCompMatch2 = fs->make<TH2F>("h_CLCT_match2","h_CLCT_match2",19,-9.5,9.5,36,0.5,36.5);
  hClctCompMatch2x = fs->make<TH2F>("h_CLCT_match2x","h_CLCT_match2x",19,-9.5,9.5,18,0.5,36.5);

  hClctCompTotal = fs->make<TH1F>("h_CLCT_total","h_CLCT_total",600,0.5,600.5);
  hClctCompTotal2 = fs->make<TH2F>("h_CLCT_total2","h_CLCT_total2",19,-9.5,9.5,36,0.5,36.5);
  hClctCompTotal2x = fs->make<TH2F>("h_CLCT_total2x","h_CLCT_total2x",19,-9.5,9.5,18,0.5,36.5);

  hLCTCompFound = fs->make<TH1F>("h_LCT_found","h_LCT_found",600,0.5,600.5);
  hLCTCompFound2 = fs->make<TH2F>("h_LCT_found2","h_LCT_found2",19,-9.5,9.5,36,0.5,36.5);
  hLCTCompFound2x = fs->make<TH2F>("h_LCT_found2x","h_LCT_found2x",19,-9.5,9.5,18,0.5,36.5);

  hLCTCompSameN = fs->make<TH1F>("h_LCT_SameN","h_LCT_SameN",600,0.5,600.5);
  hLCTCompSameN2 = fs->make<TH2F>("h_LCT_SameN2","h_LCT_SameN2",19,-9.5,9.5,36,0.5,36.5);
  hLCTCompSameN2x = fs->make<TH2F>("h_LCT_SameN2x","h_LCT_SameN2x",19,-9.5,9.5,18,0.5,36.5);

  hLCTCompMatch = fs->make<TH1F>("h_LCT_match","h_LCT_match",600,0.5,600.5);
  hLCTCompMatch2 = fs->make<TH2F>("h_LCT_match2","h_LCT_match2",19,-9.5,9.5,36,0.5,36.5);
  hLCTCompMatch2x = fs->make<TH2F>("h_LCT_match2x","h_LCT_match2x",19,-9.5,9.5,18,0.5,36.5);

  hLCTCompTotal = fs->make<TH1F>("h_LCT_total","h_LCT_total",600,0.5,600.5);
  hLCTCompTotal2 = fs->make<TH2F>("h_LCT_total2","h_LCT_total2",19,-9.5,9.5,36,0.5,36.5);
  hLCTCompTotal2x = fs->make<TH2F>("h_LCT_total2x","h_LCT_total2x",19,-9.5,9.5,18,0.5,36.5);

  //Chad's improved historgrams
  hAlctCompFound2i = fs->make<TH2F>("h_ALCT_found2i","h_ALCT_found2i",18,0,18,36,0.5,36.5);
  hAlctCompSameN2i = fs->make<TH2F>("h_ALCT_SameN2i","h_ALCT_SameN2i",18,0,18,36,0.5,36.5);
  hAlctCompMatch2i = fs->make<TH2F>("h_ALCT_match2i","h_ALCT_match2i",18,0,18,36,0.5,36.5);
  hAlctCompTotal2i = fs->make<TH2F>("h_ALCT_total2i","h_ALCT_total2i",18,0,18,36,0.5,36.5);
  hClctCompFound2i = fs->make<TH2F>("h_CLCT_found2i","h_CLCT_found2i",18,0,18,36,0.5,36.5);
  hClctCompSameN2i = fs->make<TH2F>("h_CLCT_SameN2i","h_CLCT_SameN2i",18,0,18,36,0.5,36.5);
  hClctCompMatch2i = fs->make<TH2F>("h_CLCT_match2i","h_CLCT_match2i",18,0,18,36,0.5,36.5);
  hClctCompTotal2i = fs->make<TH2F>("h_CLCT_total2i","h_CLCT_total2i",18,0,18,36,0.5,36.5);
  hLCTCompFound2i = fs->make<TH2F>("h_LCT_found2i","h_LCT_found2i",18,0,18,36,0.5,36.5);
  hLCTCompSameN2i = fs->make<TH2F>("h_LCT_SameN2i","h_LCT_SameN2i",18,0,18,36,0.5,36.5);
  hLCTCompMatch2i = fs->make<TH2F>("h_LCT_match2i","h_LCT_match2i",18,0,18,36,0.5,36.5);
  hLCTCompTotal2i = fs->make<TH2F>("h_LCT_total2i","h_LCT_total2i",18,0,18,36,0.5,36.5);

  //  hAlctCompFound = fs->make<TH1F>("h_ALCT_found","h_ALCT_found",600,0.5,600.5);

  // ALCTs.
  for (int i = 0; i < MAX_ENDCAPS; i++) { // endcaps
    for (int j = 0; j < CSC_TYPES; j++) { // station/ring
      char asdf[256];

      float csc_max = NCHAMBERS[j] + 0.5;
      if (i == 0) s = "Comp_ALCTs, " + csc_type_plus[j];
      else        s = "Comp_ALCTs, " + csc_type_minus[j];
      sprintf(asdf,"Comp_ALCTsFound_%i",i*CSC_TYPES+j);
      hAlctCompFoundCsc[i][j] =
        fs->make<TH1F>(asdf, s.c_str(), NCHAMBERS[j], 0.5, csc_max);
      sprintf(asdf,"Comp_ALCTsSame_%i",i*CSC_TYPES+j);
      hAlctCompSameNCsc[i][j] =
        fs->make<TH1F>(asdf, s.c_str(), NCHAMBERS[j], 0.5, csc_max);
      sprintf(asdf,"Comp_ALCTsTotal_%i",i*CSC_TYPES+j);
      hAlctCompTotalCsc[i][j] =
        fs->make<TH1F>(asdf, s.c_str(), NCHAMBERS[j], 0.5, csc_max);
      sprintf(asdf,"Comp_ALCTsMatch_%i",i*CSC_TYPES+j);
      hAlctCompMatchCsc[i][j] =
        fs->make<TH1F>(asdf, s.c_str(), NCHAMBERS[j], 0.5, csc_max);
      hAlctCompFoundCsc[i][j]->Sumw2();
      hAlctCompSameNCsc[i][j]->Sumw2();
      hAlctCompTotalCsc[i][j]->Sumw2();
      hAlctCompMatchCsc[i][j]->Sumw2();
    }
  }

  // CLCTs.
  for (int i = 0; i < MAX_ENDCAPS; i++) { // endcaps
    for (int j = 0; j < CSC_TYPES; j++) { // station/ring
      float csc_max = NCHAMBERS[j] + 0.5;
      char asdf[256];
      if (i == 0) s = "Comp_CLCTs, " + csc_type_plus[j];
      else        s = "Comp_CLCTs, " + csc_type_minus[j];
      sprintf(asdf,"Comp_CLCTsFound_%i",i*CSC_TYPES+j);
      hClctCompFoundCsc[i][j] =
        fs->make<TH1F>(asdf, s.c_str(), NCHAMBERS[j], 0.5, csc_max);
      sprintf(asdf,"Comp_CLCTsSame_%i",i*CSC_TYPES+j);
      hClctCompSameNCsc[i][j] =
        fs->make<TH1F>(asdf, s.c_str(), NCHAMBERS[j], 0.5, csc_max);
      sprintf(asdf,"Comp_CLCTsTotal_%i",i*CSC_TYPES+j);
      hClctCompTotalCsc[i][j] =
        fs->make<TH1F>(asdf, s.c_str(), NCHAMBERS[j], 0.5, csc_max);
      sprintf(asdf,"Comp_CLCTsMatch_%i",i*CSC_TYPES+j);
      hClctCompMatchCsc[i][j] =
        fs->make<TH1F>(asdf, s.c_str(), NCHAMBERS[j], 0.5, csc_max);
      hClctCompFoundCsc[i][j]->Sumw2();
      hClctCompSameNCsc[i][j]->Sumw2();
      hClctCompTotalCsc[i][j]->Sumw2();
      hClctCompMatchCsc[i][j]->Sumw2();
    }
  }

  // Correlated LCTs.
  for (int i = 0; i < MAX_ENDCAPS; i++) { // endcaps
    for (int j = 0; j < CSC_TYPES; j++) { // station/ring
      float csc_max = NCHAMBERS[j] + 0.5;
      char asdf[256];
      if (i == 0) s = "LCTs, " + csc_type_plus[j];
      else        s = "LCTs, " + csc_type_minus[j];
      sprintf(asdf,"LCTs_CompFound_%i",i*CSC_TYPES+j);
      hLctCompFoundCsc[i][j] =
        fs->make<TH1F>(asdf, s.c_str(), NCHAMBERS[j], 0.5, csc_max);
      sprintf(asdf,"LCTs_CompSame_%i",i*CSC_TYPES+j);
      hLctCompSameNCsc[i][j] =
        fs->make<TH1F>(asdf, s.c_str(), NCHAMBERS[j], 0.5, csc_max);
      sprintf(asdf,"LCTs_CompTotal_%i",i*CSC_TYPES+j);
      hLctCompTotalCsc[i][j] =
        fs->make<TH1F>(asdf, s.c_str(), NCHAMBERS[j], 0.5, csc_max);
      sprintf(asdf,"LCTs_CompMatch_%i",i*CSC_TYPES+j);
      hLctCompMatchCsc[i][j] =
        fs->make<TH1F>(asdf, s.c_str(), NCHAMBERS[j], 0.5, csc_max);
      hLctCompFoundCsc[i][j]->Sumw2();
      hLctCompSameNCsc[i][j]->Sumw2();
      hLctCompTotalCsc[i][j]->Sumw2();
      hLctCompMatchCsc[i][j]->Sumw2();
    }
  }

  bookedCompHistos = true;
  if (debug) LogTrace("CSCTriggerPrimitivesReader")
               <<"bookedCompHistos  ";
  std::cout <<"bookCompHistos "<< std::endl;
}

void CSCTriggerPrimitivesReader::bookResolHistos() {
  edm::Service<TFileService> fs;

  // Limits for resolution histograms
  const double EDMIN = -0.05; // eta min
  const double EDMAX =  0.05; // eta max
  const double PDMIN = -5.0;  // phi min (mrad)
  const double PDMAX =  5.0;  // phi max (mrad)

  hResolDeltaWG = fs->make<TH1F>("", "Delta key wiregroup", 10, -5., 5.);

  hResolDeltaHS = fs->make<TH1F>("", "Delta key halfstrip", 10, -5., 5.);
  hResolDeltaDS = fs->make<TH1F>("", "Delta key distrip",   10, -5., 5.);

  hResolDeltaEta   = fs->make<TH1F>("", "#eta_rec-#eta_sim", 100, EDMIN, EDMAX);
  hResolDeltaPhi   = fs->make<TH1F>("", "#phi_rec-#phi_sim (mrad)", 100, -10., 10.);
  hResolDeltaPhiHS = fs->make<TH1F>("", "#phi_rec-#phi_sim (mrad), halfstrips", 100, -10., 10.);
  hResolDeltaPhiDS = fs->make<TH1F>("", "#phi_rec-#phi_sim (mrad), distrips", 100, -10., 10.);

  hEtaRecVsSim = fs->make<TH2F>("", "#eta_rec vs #eta_sim", 64, 0.9,  2.5,  64, 0.9,  2.5);
  hPhiRecVsSim = fs->make<TH2F>("", "#phi_rec vs #phi_sim", 100, 0., TWOPI, 100, 0., TWOPI);

  // LCT quantities per station
  char histname[60];
  for (int i = 0; i < MAX_STATIONS; i++) {
    sprintf(histname, "ALCTs vs eta, station %d", i+1);
    hAlctVsEta[i]    = fs->make<TH1F>("", histname, 66, 0.875, 2.525);

    sprintf(histname, "CLCTs vs phi, station %d", i+1);
    hClctVsPhi[i]    = fs->make<TH1F>("", histname, 100, 0.,   TWOPI);

    sprintf(histname, "#LT#eta_rec-#eta_sim#GT, station %d", i+1);
    hEtaDiffVsEta[i] = fs->make<TH1F>("", histname, 66, 0.875, 2.525);

    sprintf(histname, "#LT#phi_rec-#phi_sim#GT, station %d", i+1);
    hPhiDiffVsPhi[i] = fs->make<TH1F>("", histname, 100, 0.,   TWOPI);
  }

  for (int i = 0; i < CSC_TYPES; i++) {
    string t0 = "#eta_rec-#eta_sim, " + csc_type[i];
    hEtaDiffCsc[i][0] = fs->make<TH1F>("", t0.c_str(), 100, EDMIN, EDMAX);
    string t1 = t0 + ", endcap1";
    hEtaDiffCsc[i][1] = fs->make<TH1F>("", t1.c_str(), 100, EDMIN, EDMAX);
    string t2 = t0 + ", endcap2";
    hEtaDiffCsc[i][2] = fs->make<TH1F>("", t2.c_str(), 100, EDMIN, EDMAX);

    string t4 = "#eta_rec-#eta_sim vs wiregroup, " + csc_type[i];
    hEtaDiffVsWireCsc[i] =
      fs->make<TH2F>("", t4.c_str(), MAX_WG[i], 0., MAX_WG[i], 100, EDMIN, EDMAX);

    string u0 = "#phi_rec-#phi_sim, " + csc_type[i];
    hPhiDiffCsc[i][0] = fs->make<TH1F>("", u0.c_str(), 100, PDMIN, PDMAX);
    string u1 = u0 + ", endcap1";
    hPhiDiffCsc[i][1] = fs->make<TH1F>("", u1.c_str(), 100, PDMIN, PDMAX);
    string u2 = u0 + ", endcap2";
    hPhiDiffCsc[i][2] = fs->make<TH1F>("", u2.c_str(), 100, PDMIN, PDMAX);
    hPhiDiffCsc[i][3] = fs->make<TH1F>("", u0.c_str(), 100, PDMIN, PDMAX);
    hPhiDiffCsc[i][4] = fs->make<TH1F>("", u0.c_str(), 100, PDMIN, PDMAX);

    int MAX_DS = MAX_HS[i]/4;
    string u5 = "#phi_rec-#phi_sim (mrad) vs distrip, " + csc_type[i];
    hPhiDiffVsStripCsc[i][0] =
      fs->make<TH2F>("", u5.c_str(), MAX_DS,    0., MAX_DS,    100, PDMIN, PDMAX);
    string u6 = "#phi_rec-#phi_sim (mrad) vs halfstrip, " + csc_type[i];
    hPhiDiffVsStripCsc[i][1] =
      fs->make<TH2F>("", u6.c_str(), MAX_HS[i], 0., MAX_HS[i], 100, PDMIN, PDMAX);

    string u7 = "#phi(layer 1)-#phi(layer 6), mrad, " + csc_type[i];
    hTrueBendCsc[i] =fs->make<TH1F>("", u7.c_str(), 100, -10., 10.);
  }

  int max_patterns, phibend;
  max_patterns = CSCConstants::NUM_CLCT_PATTERNS;
  for (int i = 0; i < max_patterns; i++) {
    phibend = ptype_TMB07[i];
    sprintf(histname, "#phi_rec-#phi_sim, bend = %d", phibend);
    hPhiDiffPattern[i] = fs->make<TH1F>("", histname, 100, PDMIN, PDMAX);
  }

  bookedResolHistos = true;
  if (debug) LogTrace("CSCTriggerPrimitivesReader")
               <<"bookedResolHistos  ";
}

void CSCTriggerPrimitivesReader::bookEfficHistos() {
  edm::Service<TFileService> fs;

  // Efficiencies per station.
  char histname[60];
  for (int i = 0; i < MAX_STATIONS; i++) {
    sprintf(histname, "SimHits vs eta, station %d", i+1);
    hEfficHitsEta[i] = fs->make<TH1F>("", histname, 66, 0.875, 2.525);

    sprintf(histname, "ALCTs vs eta, station %d", i+1);
    hEfficALCTEta[i] = fs->make<TH1F>("", histname, 66, 0.875, 2.525);

    sprintf(histname, "CLCTs vs eta, station %d", i+1);
    hEfficCLCTEta[i] = fs->make<TH1F>("", histname, 66, 0.875, 2.525);
  }

  // Efficiencies per chamber type.
  for (int i = 0; i < CSC_TYPES; i++) {
    string t0 = "SimHits vs eta, " + csc_type[i];
    hEfficHitsEtaCsc[i] = fs->make<TH1F>("", t0.c_str(), 66, 0.875, 2.525);
    string t1 = "ALCTs vs eta, " + csc_type[i];
    hEfficALCTEtaCsc[i] = fs->make<TH1F>("", t1.c_str(), 66, 0.875, 2.525);
    string t2 = "CLCTs vs eta, " + csc_type[i];
    hEfficCLCTEtaCsc[i] = fs->make<TH1F>("", t1.c_str(), 66, 0.875, 2.525);
  }

  bookedEfficHistos = true;
  if (debug) LogTrace("CSCTriggerPrimitivesReader")
               <<"bookedEfficHistos  ";
  std::cout <<" bookedEfficHistos " << std::endl;
}

void CSCTriggerPrimitivesReader::fillALCTHistos(const CSCALCTDigiCollection* alcts) {
  // Book histos when called for the first time.
  if (!bookedALCTHistos) bookALCTHistos();

  int nValidALCTs = 0;
  for (auto detUnitIt = alcts->begin(); detUnitIt != alcts->end(); detUnitIt++) {
    int nValidALCTsPerCSC = 0;
    const CSCDetId& id = (*detUnitIt).first;
    if (checkBadChambers_ && badChambers_->isInBadChamber(id)) continue;
    const auto& range = (*detUnitIt).second;
    for (auto digiIt = range.first;
         digiIt != range.second; digiIt++) {

      bool alct_valid = (*digiIt).isValid();
      hAlctValid->Fill(alct_valid);
      if (alct_valid) {
        hAlctQuality->Fill((*digiIt).getQuality());
        hAlctAccel->Fill((*digiIt).getAccelerator());
        hAlctCollis->Fill((*digiIt).getCollisionB());
        hAlctKeyGroup->Fill((*digiIt).getKeyWG());
        hAlctBXN->Fill((*digiIt).getBX());

        int csctype = getCSCType(id);
        hAlctPerCSC->Fill(csctype);
        hAlctCsc[id.endcap()-1][csctype]->Fill(id.chamber());

        if (csctype==0) hAlctKeyGroupME11->Fill((*digiIt).getKeyWG());

        nValidALCTs++;
        nValidALCTsPerCSC++;

        if (debug) LogTrace("CSCTriggerPrimitivesReader")
                     << (*digiIt) << " found in ME" << ((id.endcap() == 1) ? "+" : "-")
                     << id.station() << "/" << id.ring() << "/" << id.chamber()
                     << " (sector " << id.triggerSector()
                     << " trig id. " << id.triggerCscId() << ")";
        //cout << "raw id = " << id.rawId() << endl;
      }
    }
    hAlctPerChamber->Fill(nValidALCTsPerCSC);
  }
  hAlctPerEvent->Fill(nValidALCTs);
  if (debug) LogTrace("CSCTriggerPrimitivesReader")
               << nValidALCTs << " valid ALCTs found in this event";
  numALCT += nValidALCTs;
}

void CSCTriggerPrimitivesReader::fillCLCTHistos(const CSCCLCTDigiCollection* clcts) {
  // Book histos when called for the first time.
  if (!bookedCLCTHistos) bookCLCTHistos();

  int nValidCLCTs = 0;
  for (auto detUnitIt = clcts->begin(); detUnitIt != clcts->end(); detUnitIt++) {
    int nValidCLCTsPerCSC = 0;
    const CSCDetId& id = (*detUnitIt).first;
    if (checkBadChambers_ && badChambers_->isInBadChamber(id)) continue;
    const auto& range = (*detUnitIt).second;
    for (auto digiIt = range.first;
         digiIt != range.second; digiIt++) {

      bool clct_valid = (*digiIt).isValid();
      hClctValid->Fill(clct_valid);
      if (clct_valid) {
        int striptype = (*digiIt).getStripType();
        int keystrip  = (*digiIt).getKeyStrip(); // halfstrip #
        if (striptype == 0) keystrip /= 4;       // distrip # for distrip ptns
        hClctQuality->Fill((*digiIt).getQuality());
        hClctStripType->Fill(striptype);
        hClctSign->Fill((*digiIt).getBend());
        hClctCFEB->Fill((*digiIt).getCFEB());
        hClctBXN->Fill((*digiIt).getBX());
        hClctKeyStrip[striptype]->Fill(keystrip);
        hClctPattern[striptype]->Fill((*digiIt).getPattern());

        int csctype = getCSCType(id);
        hClctPerCSC->Fill(csctype);
        hClctCsc[id.endcap()-1][csctype]->Fill(id.chamber());

        if (striptype==1 && csctype==0) hClctKeyStripME11->Fill(keystrip);

        int phibend;
        int pattern = (*digiIt).getPattern();
        phibend = ptype_TMB07[pattern];
        hClctBendCsc[csctype][striptype]->Fill(abs(phibend));

        if (striptype == 0) // distrips
          hClctKeyStripCsc[csctype]->Fill(keystrip);

        nValidCLCTs++;
        nValidCLCTsPerCSC++;

        if (debug) LogTrace("CSCTriggerPrimitivesReader")
                     << (*digiIt) << " found in ME" << ((id.endcap() == 1) ? "+" : "-")
                     << id.station() << "/" << id.ring() << "/" << id.chamber()
                     << " (sector " << id.triggerSector()
                     << " trig id. " << id.triggerCscId() << ")";
      }
    }
    hClctPerChamber->Fill(nValidCLCTsPerCSC);
  }
  hClctPerEvent->Fill(nValidCLCTs);
  if (debug) LogTrace("CSCTriggerPrimitivesReader")
               << nValidCLCTs << " valid CLCTs found in this event";
  numCLCT += nValidCLCTs;
}

void CSCTriggerPrimitivesReader::fillLCTTMBHistos(const CSCCorrelatedLCTDigiCollection* lcts) {
  // Book histos when called for the first time.
  if (!bookedLCTTMBHistos) bookLCTTMBHistos();

  int nValidLCTs = 0;
  bool alct_valid, clct_valid;

  for (auto detUnitIt = lcts->begin(); detUnitIt != lcts->end(); detUnitIt++) {
    int nValidLCTsPerCSC = 0;
    const CSCDetId& id = (*detUnitIt).first;
    if (checkBadChambers_ && badChambers_->isInBadChamber(id)) continue;
    const auto& range = (*detUnitIt).second;
    for (auto digiIt = range.first;
         digiIt != range.second; digiIt++) {

      bool lct_valid = (*digiIt).isValid();
      hLctTMBValid->Fill(lct_valid);
      if (lct_valid) {
        hLctTMBEndcap->Fill(id.endcap());
        hLctTMBStation->Fill(id.station());
        hLctTMBSector->Fill(id.triggerSector());
        hLctTMBRing->Fill(id.ring());
        hLctTMBChamber[id.station()-1]->Fill(id.triggerCscId());

        int quality = (*digiIt).getQuality();
        hLctTMBQuality->Fill(quality);
        hLctTMBBXN->Fill((*digiIt).getBX());

        alct_valid = (quality != 0 && quality != 2);
        if (alct_valid) {
          hLctTMBKeyGroup->Fill((*digiIt).getKeyWG());
        }

        clct_valid = (quality != 0 && quality != 1);
        if (clct_valid) {
          hLctTMBKeyStrip->Fill((*digiIt).getStrip());
          hLctTMBStripType->Fill(1.);
          hLctTMBPattern->Fill((*digiIt).getPattern());
          hLctTMBBend->Fill((*digiIt).getBend());
        }

        int csctype = getCSCType(id);
        hLctTMBPerCSC->Fill(csctype);
        hLctTMBCsc[id.endcap()-1][csctype]->Fill(id.chamber());
        // Truly correlated LCTs; for DAQ
        if (alct_valid && clct_valid) hCorrLctTMBPerCSC->Fill(csctype);

        if (alct_valid && csctype==0) {
          hLctTMBKeyGroupME11->Fill((*digiIt).getKeyWG());
        }
        if (clct_valid && csctype==0) {
          hLctTMBKeyStripME11->Fill((*digiIt).getStrip());
        }

        nValidLCTs++;
        nValidLCTsPerCSC++;

        if (debug) LogTrace("CSCTriggerPrimitivesReader")
                     << (*digiIt) << " found in ME" << ((id.endcap() == 1) ? "+" : "-")
                     << id.station() << "/" << id.ring() << "/" << id.chamber()
                     << " (sector " << id.triggerSector()
                     << " trig id. " << id.triggerCscId() << ")";
      }
    }
    hLctTMBPerChamber->Fill(nValidLCTsPerCSC);
  }
  hLctTMBPerEvent->Fill(nValidLCTs);
  if (debug) LogTrace("CSCTriggerPrimitivesReader")
               << nValidLCTs << " valid LCTs found in this event";
  numLCTTMB += nValidLCTs;
}

void CSCTriggerPrimitivesReader::fillLCTMPCHistos(const CSCCorrelatedLCTDigiCollection* lcts) {
  // Book histos when called for the first time.
  if (!bookedLCTMPCHistos) bookLCTMPCHistos();

  int nValidLCTs = 0;
  bool alct_valid, clct_valid;

  for (auto detUnitIt = lcts->begin(); detUnitIt != lcts->end(); detUnitIt++) {
    const CSCDetId& id = (*detUnitIt).first;
    if (checkBadChambers_ && badChambers_->isInBadChamber(id)) continue;
    const auto& range = (*detUnitIt).second;
    for (auto digiIt = range.first;
         digiIt != range.second; digiIt++) {

      bool lct_valid = (*digiIt).isValid();
      hLctMPCValid->Fill(lct_valid);
      if (lct_valid) {
        hLctMPCEndcap->Fill(id.endcap());
        hLctMPCStation->Fill(id.station());
        hLctMPCSector->Fill(id.triggerSector());
        hLctMPCRing->Fill(id.ring());
        hLctMPCChamber[id.station()-1]->Fill(id.triggerCscId());

        int quality = (*digiIt).getQuality();
        hLctMPCQuality->Fill(quality);
        hLctMPCBXN->Fill((*digiIt).getBX());

        alct_valid = (quality != 0 && quality != 2);
        if (alct_valid) {
          hLctMPCKeyGroup->Fill((*digiIt).getKeyWG());
        }

        clct_valid = (quality != 0 && quality != 1);
        if (clct_valid) {
          hLctMPCKeyStrip->Fill((*digiIt).getStrip());
          hLctMPCStripType->Fill(1.);
          hLctMPCPattern->Fill((*digiIt).getPattern());
          hLctMPCBend->Fill((*digiIt).getBend());
        }

        int csctype = getCSCType(id);
        hLctMPCPerCSC->Fill(csctype);
        // Truly correlated LCTs; for DAQ
        if (alct_valid && clct_valid) hCorrLctMPCPerCSC->Fill(csctype);


        if (alct_valid && csctype==0) {
          hLctMPCKeyGroupME11->Fill((*digiIt).getKeyWG());
        }
        if (clct_valid && csctype==0) {
          hLctMPCKeyStripME11->Fill((*digiIt).getStrip());
        }

        nValidLCTs++;

        if (debug) LogTrace("CSCTriggerPrimitivesReader")
                     << "MPC "
                     << (*digiIt) << " found in ME" << ((id.endcap() == 1) ? "+" : "-")
                     << id.station() << "/" << id.ring() << "/" << id.chamber()
                     << " (sector " << id.triggerSector()
                     << " trig id. " << id.triggerCscId() << ")";
      }
    }
  }
  hLctMPCPerEvent->Fill(nValidLCTs);
  if (debug) LogTrace("CSCTriggerPrimitivesReader")
               << nValidLCTs << " MPC LCTs found in this event";
  numLCTMPC += nValidLCTs;
}

void CSCTriggerPrimitivesReader::compare(const CSCALCTDigiCollection* alcts_data,
                                         const CSCALCTDigiCollection* alcts_emul,
                                         const CSCCLCTDigiCollection* clcts_data,
                                         const CSCCLCTDigiCollection* clcts_emul,
                                         const CSCCLCTPreTriggerDigiCollection* pretrigs_emul,
                                         const CSCCorrelatedLCTDigiCollection* lcts_data,
                                         const CSCCorrelatedLCTDigiCollection* lcts_emul){

  // Book histos when called for the first time.
  if (!bookedCompHistos) bookCompHistos();

  // Comparisons
  compareALCTs(alcts_data, alcts_emul);
  compareCLCTs(clcts_data, clcts_emul, pretrigs_emul);
  compareLCTs(lcts_data,  lcts_emul, alcts_data, clcts_data);
  //compareMPCLCTs(mpclcts_data,  mpclcts_emul, alcts_data, clcts_data);
}

void CSCTriggerPrimitivesReader::compareALCTs(const CSCALCTDigiCollection* alcts_data,
                                              const CSCALCTDigiCollection* alcts_emul) {
  int emul_corr_bx;

  // Should be taken from config. parameters.
  int fifo_pretrig     = 10;
  int fpga_latency     =  6;
  int l1a_window_width =  7;
  // Time offset of raw hits w.r.t. the full 12-bit BXN.
  int rawhit_tbin_offset =
    (fifo_pretrig - fpga_latency) + (l1a_window_width-1)/2;
  // Extra difference due to additional register stages; determined
  // empirically.
  int register_delay =  2;
  // Loop over all chambers in search for ALCTs.
  std::vector<CSCALCTDigi>::const_iterator pd, pe;
  for (int endc = 1; endc <= 2; endc++) {
    for (int stat = 1; stat <= 4; stat++) {
      for (int ring = 1; ring <= maxRing(stat); ring++) {
        for (int cham = 1; cham <= 36; cham++) {
          // Calculate DetId.  0th layer means whole chamber.
          CSCDetId detid(endc, stat, ring, cham, 0);

          // Skip chambers marked as bad.
          if (checkBadChambers_ && badChambers_->isInBadChamber(detid)) continue;
          std::vector<CSCALCTDigi> alctV_data, alctV_emul;
          const auto& drange = alcts_data->get(detid);
          for (auto digiIt = drange.first; digiIt != drange.second; digiIt++) {
            if ((*digiIt).isValid()) {
              alctV_data.push_back(*digiIt);
            }
          }

          const auto& erange = alcts_emul->get(detid);
          std::vector<bool> bookedalctV_emul;
          for (auto digiIt = erange.first; digiIt != erange.second; digiIt++) {
            if ((*digiIt).isValid()) {
              alctV_emul.push_back(*digiIt);
              bookedalctV_emul.push_back(false);
            }
          }

          int ndata = alctV_data.size();
          int nemul = alctV_emul.size();

          if (ndata == 0 && nemul == 0) continue;
          bool firstfill = true;

          if (debug or (stat == 1 and ring== 1 and (ndata != nemul)) ) {
            ostringstream strstrm;
            strstrm << "\n--- ME" << ((detid.endcap() == 1) ? "+" : "-")
                    << detid.station() << "/" << detid.ring() << "/"
                    << detid.chamber()
                    << " (sector "  << detid.triggerSector()
                    << " trig id. " << detid.triggerCscId() << "):\n";
            strstrm << "  **** " << ndata << " valid data ALCTs found:\n";
            for (pd = alctV_data.begin(); pd != alctV_data.end(); pd++) {
              strstrm << "     " << (*pd)
                      << " Full BX = " << (*pd).getFullBX() << "\n";
            }
            strstrm << "  **** " << nemul << " valid emul ALCTs found:\n";
            for (pe = alctV_emul.begin(); pe != alctV_emul.end(); pe++) {
              strstrm << "     " << (*pe);
              for (pd = alctV_data.begin(); pd != alctV_data.end(); pd++) {
                if ((*pd).getTrknmb() == (*pe).getTrknmb()) {
                  int emul_bx = (*pe).getBX();
                  emul_corr_bx =
                    emul_bx - rawhit_tbin_offset + register_delay;
                  strstrm << " Corr BX = " << emul_corr_bx;
                  break;
                }
              }
              strstrm << "\n";
            }
            LogTrace("CSCTriggerPrimitivesReader") << strstrm.str();
            //if (stat==1 && ring==1)
            //  std::cout <<"ME11  CompareALCTs "<< strstrm.str()<< std::endl;
          }

          //Data, add HS quality later
          int perEv_nStub_data = 0;
          int perEv_nStub_emul = 0;
          perStub[0].init(RUN_, Event_);
          perStub[1].init(RUN_, Event_);
          for (pd = alctV_data.begin(); pd != alctV_data.end(); pd++){
            perEv_nStub_data++;
          }
          perStub[0].t_endcap = endc;
          perStub[0].t_station = stat;
          perStub[0].t_chamber = cham;
          perStub[0].t_ring = ring;
          perStub[0].t_EventNumberAnalyzed = eventsAnalyzed;
          perStub[0].t_nStubs              = ndata;
          perStub[0].t_nStubs_readout              = ndata;
          event_tree[0]->Fill();
          //Emul
          for (pe = alctV_emul.begin(); pe != alctV_emul.end(); pe++){
            perEv_nStub_emul++;
          }
          perStub[1].t_endcap = endc;
          perStub[1].t_station = stat;
          perStub[1].t_chamber = cham;
          perStub[1].t_ring = ring;
          perStub[1].t_EventNumberAnalyzed = eventsAnalyzed;
          perStub[1].t_nStubs              = nemul;
          perStub[1].t_nStubs_readout      = nemul;
          event_tree[1]->Fill();

          int csctype = getCSCType(detid);
          hAlctCompFoundCsc[endc-1][csctype]->Fill(cham);
          int mychamber = chamberSerial(detid);
          hAlctCompFound->Fill(mychamber);
          int ix = chamberIX(detid);
          int ix2 = chamberIXi(detid);
          //      printf("station %i, ring %i, chamber %i, ix+(detid.ring()-1) %i\n",
          //         detid.station(),detid.ring(),detid.chamber(),ix);
          if(detid.station()>1 && detid.ring()==1) {
            hAlctCompFound2x->Fill(ix,detid.chamber()*2);
          }
          else {
            hAlctCompFound2->Fill(ix,detid.chamber());
          }
          hAlctCompFound2i->Fill(ix2,detid.chamber());

          if (ndata != nemul) {
            LogTrace("CSCTriggerPrimitivesReader")
              << "   +++ Different numbers of ALCTs found in ME"
              << ((endc == 1) ? "+" : "-") << stat << "/"
              << ring << "/" << cham
              << ": data = " << ndata << " emulator = " << nemul << " +++\n";
          }
          else {
            hAlctCompSameNCsc[endc-1][csctype]->Fill(cham);
            if(detid.station()>1 && detid.ring()==1) {
              hAlctCompSameN2x->Fill(ix,detid.chamber()*2);
            }
            else {
              hAlctCompSameN2->Fill(ix,detid.chamber());
            }
            hAlctCompSameN2i->Fill(ix2,detid.chamber());
          }

          for (int i = 0; i < ndata; i++) {
            if (alctV_data[i].isValid() == 0) continue;
            int data_trknmb    = alctV_data[i].getTrknmb();
            //int data_quality   = alctV_data[i].getQuality();
            //int data_accel     = alctV_data[i].getAccelerator();
            //int data_collB     = alctV_data[i].getCollisionB();
            int data_wiregroup = alctV_data[i].getKeyWG();
            int data_bx        = alctV_data[i].getBX();

            stubs_comparison[0].init(RUN_, Event_);
            stubs_comparison[0].firstfill = firstfill;
            if (firstfill) firstfill = false;
            stubs_comparison[0].nEvents = eventsAnalyzed;
            stubs_comparison[0].endcap = endc;
            stubs_comparison[0].station = stat;
            stubs_comparison[0].ring = ring;
            stubs_comparison[0].chamber = cham;
            stubs_comparison[0].chambertype = detid.iChamberType();
            stubs_comparison[0].totStubs_data = ndata;
            stubs_comparison[0].totStubs_emul = nemul;
            stubs_comparison[0].totStubs_emul_readout = nemul;
            stubs_comparison[0].nStub_data = i+1;
            stubs_comparison[0].has_data = true;
            stubs_comparison[0].quality_data = alctV_data[i].getQuality();
            stubs_comparison[0].key_WG_data = alctV_data[i].getKeyWG();
            stubs_comparison[0].bx_data = data_bx;
            stubs_comparison[0].fullbx_data = alctV_data[i].getFullBX();
            stubs_comparison[0].trknmb_data = data_trknmb;
            GlobalPoint gp_alct_data(getGlobalPosition(detid.rawId(), data_wiregroup, 60));
            stubs_comparison[0].eta_data = gp_alct_data.eta();
            for (int j=0; j<nemul; j++){
              if (alctV_emul[j].isValid() == 0) continue;
              if (bookedalctV_emul[j]) continue; //used alct

              int emul_trknmb    = alctV_emul[j].getTrknmb();
              //int emul_quality   = alctV_emul[j].getQuality();
              //int emul_accel     = alctV_emul[j].getAccelerator();
              //int emul_collB     = alctV_emul[j].getCollisionB();
              int emul_wiregroup = alctV_emul[j].getKeyWG();
              int emul_bx        = alctV_emul[j].getBX();
              // Emulator BX re-calculated for comparison with BX in the data.
              emul_corr_bx = emul_bx - rawhit_tbin_offset + register_delay;
              if (dataIsAnotherMC_)
                emul_corr_bx = emul_bx;

              if (ndata == nemul) {
                hAlctCompTotal->Fill(mychamber);
                hAlctCompTotalCsc[endc-1][csctype]->Fill(cham);
                if(detid.station()>1 && detid.ring()==1) {
                  hAlctCompTotal2x->Fill(ix,detid.chamber()*2);
                }
                else {
                  hAlctCompTotal2->Fill(ix,detid.chamber());
                }
                hAlctCompTotal2i->Fill(ix2,detid.chamber());
              }
              if (data_trknmb    == emul_trknmb )  {
                if (abs(data_wiregroup - emul_wiregroup) <= 2){
                  hAlctCompMatchCsc[endc-1][csctype]->Fill(cham);
                  hAlctCompMatch->Fill(mychamber);
                  if(detid.station()>1 && detid.ring()==1) {
                    hAlctCompMatch2x->Fill(ix,detid.chamber()*2);
                  }
                  else {
                    hAlctCompMatch2->Fill(ix,detid.chamber());
                  }
                  hAlctCompMatch2i->Fill(ix2,detid.chamber());
                }
                if (debug) LogTrace("CSCTriggerPrimitivesReader")
                             << "       Identical ALCTs #" << data_trknmb;
                stubs_comparison[0].nStub_emul = j+1;
                stubs_comparison[0].has_emul = true;
                stubs_comparison[0].quality_emul = alctV_emul[j].getQuality();
                stubs_comparison[0].key_WG_emul = alctV_emul[j].getKeyWG();
                stubs_comparison[0].bx_emul = alctV_emul[j].getBX();
                stubs_comparison[0].trknmb_emul = emul_trknmb;
                GlobalPoint gp_alct_emul(getGlobalPosition(detid.rawId(), emul_wiregroup, 60));
                stubs_comparison[0].eta_emul = gp_alct_emul.eta();
                stubs_comparison[0].bx_corr_emul = emul_corr_bx;
                bookedalctV_emul[j] = true;
                break;
              }
              else {
                LogTrace("CSCTriggerPrimitivesReader")
                  << "       Different ALCTs #" << data_trknmb << " in ME"
                  << ((endc == 1) ? "+" : "-") << stat << "/"
                  << ring << "/" << cham;
              }

            }//loop emul
            if (debug and stubs_comparison[0].key_WG_data != stubs_comparison[0].key_WG_emul)
              LogTrace("CSCTriggerPrimitivesReader")
                <<" not matched ALCT from data "<< alctV_data[i]
                <<" stubs_comparison 0 key_WG_data "<<stubs_comparison[0].key_WG_data <<" key_WG_emul "<< stubs_comparison[0].key_WG_emul;
            //if (stat==1) std::cout <<" stub_tree filled , ring "<< stubs_comparison[0].ring << std::endl;
            //cout <<"ALCT data BX "<< stubs_comparison[0].bx_data <<" WG "<< stubs_comparison[0].key_WG_data <<" emul BX "<< stubs_comparison[0].bx_emul<<" emul BX corrected "<< stubs_comparison[0].bx_corr_emul <<" WG "<< stubs_comparison[0].key_WG_emul << endl;
            stub_tree[0]->Fill();
          }//loop data
          for (int i = 0; i<nemul; i++){
            if (alctV_emul[i].isValid() == 0 or bookedalctV_emul[i]) continue;
            int emul_bx        = alctV_emul[i].getBX();
            stubs_comparison[0].init(RUN_, Event_);
            stubs_comparison[0].firstfill = firstfill;
            if (firstfill) firstfill = false;
            stubs_comparison[0].nEvents = eventsAnalyzed;
            stubs_comparison[0].endcap = endc;
            stubs_comparison[0].station = stat;
            stubs_comparison[0].ring = ring;
            stubs_comparison[0].chamber = cham;
            stubs_comparison[0].chambertype = detid.iChamberType();
            stubs_comparison[0].totStubs_data = ndata;
            stubs_comparison[0].totStubs_emul = nemul;
            stubs_comparison[0].totStubs_emul_readout = nemul;
            stubs_comparison[0].nStub_data = -1;
            stubs_comparison[0].nStub_emul = i+1;
            stubs_comparison[0].has_data = false;
            stubs_comparison[0].has_emul = true;
            stubs_comparison[0].trknmb_emul = alctV_emul[i].getTrknmb();
            stubs_comparison[0].quality_emul = alctV_emul[i].getQuality();
            stubs_comparison[0].key_WG_emul = alctV_emul[i].getKeyWG();
            stubs_comparison[0].bx_emul = alctV_emul[i].getBX();
            stubs_comparison[0].fullbx_emul = alctV_emul[i].getFullBX();
            GlobalPoint gp_alct_emul(getGlobalPosition(detid.rawId(), alctV_emul[i].getKeyWG(), 60));
            stubs_comparison[0].eta_emul = gp_alct_emul.eta();
            // Emulator BX re-calculated for comparison with BX in the data.
            emul_corr_bx = emul_bx - rawhit_tbin_offset + register_delay;
            if (dataIsAnotherMC_)
              emul_corr_bx = emul_bx;
            stubs_comparison[0].bx_corr_emul = emul_corr_bx;
            stub_tree[0]->Fill();
            if (debug)
              LogTrace("CSCTriggerPrimitivesReader")
                <<"not matched ALCT from emulation "<< alctV_emul[i];

          }//loop emul
        }
      }
    }
  }
}// compare ALCT

void CSCTriggerPrimitivesReader::compareCLCTs(const CSCCLCTDigiCollection* clcts_data,
                                              const CSCCLCTDigiCollection* clcts_emul,
                                              const CSCCLCTPreTriggerDigiCollection* pretrigs_emul) {
  // Number of Tbins before pre-trigger for raw cathode hits.
  const int tbin_cathode_offset = 7;
  //const int tbin_cathode_offset = 8;//in MC, it became 8, Tao
  const int pretrig_trig_zone = 5;// max distance between CLCT key hs and pretrigger hs

  // Loop over all chambers in search for CLCTs.
  std::vector<CSCCLCTDigi>::const_iterator pd, pe;
  std::vector<CSCCLCTPreTriggerDigi>::const_iterator pretrig;
  perStub[2].init(RUN_, Event_);
  perStub[3].init(RUN_, Event_);
  for (int endc = 1; endc <= 2; endc++) {
    for (int stat = 1; stat <= 4; stat++) {
      for (int ring = 1; ring <= maxRing(stat); ring++) {
        for (int cham = 1; cham <= 36; cham++) {
          // Calculate DetId.  0th layer means whole chamber.
          CSCDetId detid(endc, stat, ring, cham, 0);

          // Skip chambers marked as bad.
          if (checkBadChambers_ && badChambers_->isInBadChamber(detid)) continue;

          std::vector<CSCCLCTDigi> clctV_data, clctV_emul;
          const auto& drange = clcts_data->get(detid);
          for (auto digiIt = drange.first; digiIt != drange.second; digiIt++) {
            if ((*digiIt).isValid()) {
              clctV_data.push_back(*digiIt);
            }
          }

          const auto& erange = clcts_emul->get(detid);
          std::vector<bool> bookedclctV_emul;
          for (auto digiIt = erange.first; digiIt != erange.second; digiIt++) {
            if ((*digiIt).isValid()) {
              for (auto pclct : clctV_emul){
                if (digiIt->getBX() != pclct.getBX() and abs(digiIt->getBX() - pclct.getBX())< 5)
                  LogTrace("CSCTriggerPrimitivesReader") <<"Two CLCTs very close in timing!!! Special event: first clct "<< pclct <<" second clct "<<*digiIt << std::endl;
              }
              clctV_emul.push_back(*digiIt);
              bookedclctV_emul.push_back(false);
            }
          }

          std::vector<CSCCLCTPreTriggerDigi> pretrigV_emul;
          const auto& pretrigrange = pretrigs_emul->get(detid);
          for (auto pretrigIt = pretrigrange.first; pretrigIt != pretrigrange.second; pretrigIt++){
            if ((*pretrigIt).isValid()){
              pretrigV_emul.push_back(*pretrigIt);
            }
          }

          int ndata = clctV_data.size();
          int nemul = clctV_emul.size();
          if (ndata == 0 && nemul == 0) continue;
          bool firstfill = true;
          int nemul_readout = 0;
          for (pe = clctV_emul.begin(); pe != clctV_emul.end(); pe++) {
            for (pd = clctV_data.begin(); pd != clctV_data.end(); pd++) {
              int emul_bx = (*pe).getBX();
              int corr_bx =
                ((*pd).getFullBX() + emul_bx - tbin_cathode_offset) & 0x03;
              int bx_data = pd->getBX();
              //if (corr_bx == bx_data or abs(pe->getKeyStrip() - pd->getKeyStrip())<=1){//if emulated BX after correction is same as data bx, it will be readout
              if (corr_bx == bx_data){//if emulated BX after correction is same as data bx, it will be readout
                nemul_readout++;
                break;
              }
            }
          }

          if (debug or nemul>ndata or (ndata != nemul) or ndata != nemul_readout) {
            LogTrace("CSCTriggerPrimitivesReader") << " CLCTs from data "<< ndata <<" CLCTs from emul "<< nemul <<" readout "<< nemul_readout;
            ostringstream strstrm;
            strstrm << "\n--- ME" << ((detid.endcap() == 1) ? "+" : "-")
                    << detid.station() << "/" << detid.ring() << "/"
                    << detid.chamber()
                    << " (sector "  << detid.triggerSector()
                    << " trig id. " << detid.triggerCscId() << "):\n";
            strstrm << "  **** " << ndata << " valid data CLCTs found:\n";
            for (pd = clctV_data.begin(); pd != clctV_data.end(); pd++) {
              strstrm << "     " << (*pd)
                      << " Full BX = " << (*pd).getFullBX() << "\n";
            }
            strstrm << "  **** " << nemul << " valid emul CLCTs found:\n";
            for (pe = clctV_emul.begin(); pe != clctV_emul.end(); pe++) {
              strstrm << "     " << (*pe);
              for (pd = clctV_data.begin(); pd != clctV_data.end(); pd++) {
                if ((*pd).getTrknmb() == (*pe).getTrknmb() or abs((*pe).getKeyStrip() - (*pd).getKeyStrip())<5) {
                  int emul_bx = (*pe).getBX();
                  int corr_bx =
                    ((*pd).getFullBX() + emul_bx - tbin_cathode_offset) & 0x03;
                  strstrm << " Corr BX = " << corr_bx;
                  break;
                }
              }
              strstrm << "\n";
            }
            LogTrace("CSCTriggerPrimitivesReader") << strstrm.str();
          }
          if (nemul_readout > 2){ //reduce nemul_readout to 2 by hand
            //cout <<"CLCT matching nemul readout is "<< nemul_readout <<", larger than 2. reduce it to 2 by hand"<<endl;
            nemul_readout = 2;
          }


          //Per event Fill, From Luca
          //Data, add HS quality later
          int perEv_nStub_data = 0;
          int perEv_nStub_emul = 0;
          perStub[2].init(RUN_, Event_);
          perStub[3].init(RUN_, Event_);
          for (pd = clctV_data.begin(); pd != clctV_data.end(); pd++){
            perEv_nStub_data++;
          }
          perStub[2].t_endcap = endc;
          perStub[2].t_station = stat;
          perStub[2].t_chamber = cham;
          perStub[2].t_ring = ring;
          perStub[2].t_EventNumberAnalyzed = eventsAnalyzed;
          perStub[2].t_nStubs              = ndata;
          perStub[2].t_nStubs_readout              = ndata;
          event_tree[2]->Fill();
          //Emul
          for (pe = clctV_emul.begin(); pe != clctV_emul.end(); pe++){
            perEv_nStub_emul++;
          }
          perStub[3].t_endcap = endc;
          perStub[3].t_station = stat;
          perStub[3].t_chamber = cham;
          perStub[3].t_ring = ring;
          perStub[3].t_EventNumberAnalyzed = eventsAnalyzed;
          perStub[3].t_nStubs              = nemul;
          perStub[3].t_nStubs_readout      = nemul_readout;
          event_tree[3]->Fill();


          int csctype = getCSCType(detid);
          hClctCompFoundCsc[endc-1][csctype]->Fill(cham);
          int ix = chamberIX(detid);
          int ix2 = chamberIXi(detid);
          if(detid.station()>1 && detid.ring()==1) {
            hClctCompFound2x->Fill(ix,detid.chamber()*2);
          }
          else {
            hClctCompFound2->Fill(ix,detid.chamber());
          }
          hClctCompFound2i->Fill(ix2,detid.chamber());
          if (ndata != nemul) {
            LogTrace("CSCTriggerPrimitivesReader")
              //cerr
              << "   +++ Different numbers of CLCTs found in ME"
              << ((endc == 1) ? "+" : "-") << stat << "/"
              << ring << "/" << cham
              << ": data = " << ndata << " emulator = " << nemul << " +++\n";
          }
          else {
            hClctCompSameNCsc[endc-1][csctype]->Fill(cham);
            if(detid.station()>1 && detid.ring()==1) {
              hClctCompSameN2x->Fill(ix,detid.chamber()*2);
            }
            else {
              hClctCompSameN2->Fill(ix,detid.chamber());
            }
            hClctCompSameN2i->Fill(ix2,detid.chamber());
          }
          int i=-1;
          int testwg = 20;
          for (pd = clctV_data.begin(); pd != clctV_data.end(); pd++) {
            i++;
            if ((*pd).isValid() == 0) continue;
            int data_trknmb    = (*pd).getTrknmb();
            int data_quality   = (*pd).getQuality();
            int data_pattern   = (*pd).getPattern();
            int data_striptype = (*pd).getStripType();
            int data_bend      = (*pd).getBend();
            int data_keystrip  = (*pd).getKeyStrip();
            int data_cfeb      = (*pd).getCFEB();
            int data_bx        = (*pd).getBX();
            int fullBX = (*pd).getFullBX(); // 12-bit full BX

            if (data_keystrip >= 128 and stat == 1 and ring == 1) testwg=5;
            else testwg = 20;

            stubs_comparison[1].init(RUN_, Event_);
            stubs_comparison[1].firstfill = firstfill;
            if (firstfill) firstfill = false;
            stubs_comparison[1].nEvents = eventsAnalyzed;
            stubs_comparison[1].endcap = endc;
            stubs_comparison[1].station = stat;
            stubs_comparison[1].ring = ring;
            stubs_comparison[1].chamber = cham;
            stubs_comparison[1].chambertype = detid.iChamberType();
            stubs_comparison[1].totStubs_data = ndata;
            stubs_comparison[1].totStubs_emul = nemul;
            stubs_comparison[1].totStubs_emul_readout = nemul_readout;
            stubs_comparison[1].nStub_data = i+1;
            stubs_comparison[1].has_data = true;
            stubs_comparison[1].quality_data = (*pd).getQuality();
            stubs_comparison[1].key_hs_data = (*pd).getKeyStrip();
            stubs_comparison[1].bend_data = (*pd).getBend();
            stubs_comparison[1].pattern_data = (*pd).getPattern();
            stubs_comparison[1].bx_data = (*pd).getBX();
            stubs_comparison[1].fullbx_data = (*pd).getFullBX();
            stubs_comparison[1].trknmb_data = data_trknmb;
            GlobalPoint gp_clct_data(getGlobalPosition(detid.rawId(), testwg, (*pd).getKeyStrip()));
            stubs_comparison[1].phi_data = gp_clct_data.phi();

            int j=-1;
            for (pe = clctV_emul.begin(); pe != clctV_emul.end(); pe++) {
              j++;
              if ((*pe).isValid() == 0) continue;
              if (bookedclctV_emul[j]) continue; //used alct
              int emul_trknmb    = (*pe).getTrknmb();
              int emul_quality   = (*pe).getQuality();
              int emul_pattern   = (*pe).getPattern();
              int emul_striptype = (*pe).getStripType();
              int emul_bend      = (*pe).getBend();
              int emul_keystrip  = (*pe).getKeyStrip();
              int emul_cfeb      = (*pe).getCFEB();
              int emul_bx        = (*pe).getBX();

              if (abs(data_keystrip - emul_keystrip) <= 2){
                // Emulator BX re-calculated using 12-bit full BX number.
                // Used for comparison with BX in the data.
                int emul_corr_bx =
                  (fullBX + emul_bx - tbin_cathode_offset) & 0x03;
                //std::cout <<"CLCT data_bx "<< data_bx <<" emul_corr_bx "<< emul_corr_bx << std::endl;
                if (dataIsAnotherMC_)
                  emul_corr_bx = (emul_bx & 0x03);
                if (ndata == nemul) {
                  hClctCompTotalCsc[endc-1][csctype]->Fill(cham);
                  if(detid.station()>1 && detid.ring()==1) {
                    hClctCompTotal2x->Fill(ix,detid.chamber()*2);
                  }
                  else {
                    hClctCompTotal2->Fill(ix,detid.chamber());
                  }
                  hClctCompTotal2i->Fill(ix2,detid.chamber());
                }
                if (data_quality   == emul_quality   &&
                    data_pattern   == emul_pattern   &&
                    data_striptype == emul_striptype &&
                    data_bend      == emul_bend      &&
                    data_keystrip  == emul_keystrip  &&
                    data_cfeb      == emul_cfeb      ){
                  if (ndata == nemul) {
                    hClctCompMatchCsc[endc-1][csctype]->Fill(cham);
                    if(detid.station()>1 && detid.ring()==1) {
                      hClctCompMatch2x->Fill(ix,detid.chamber()*2);
                    }
                    else {
                      hClctCompMatch2->Fill(ix,detid.chamber());
                    }
                    hClctCompMatch2i->Fill(ix2,detid.chamber());
                  }
                  if (debug) LogTrace("CSCTriggerPrimitivesReader")
                               << "       Identical CLCTs #" << data_trknmb;
                }
                else {
                  LogTrace("CSCTriggerPrimitivesReader")
                    << "       Different CLCTs #" << data_trknmb << " in ME"
                    << ((endc == 1) ? "+" : "-") << stat << "/"
                    << ring << "/" << cham <<" data_bx "<< data_bx <<" emul_corr_bx "<< emul_corr_bx;
                }
                stubs_comparison[1].nStub_emul = j+1;
                stubs_comparison[1].has_emul = true;
                stubs_comparison[1].quality_emul = (*pe).getQuality();
                stubs_comparison[1].key_hs_emul = (*pe).getKeyStrip();
                stubs_comparison[1].bend_emul = (*pe).getBend();
                stubs_comparison[1].pattern_emul = (*pe).getPattern();
                stubs_comparison[1].bx_emul = (*pe).getBX();
                stubs_comparison[1].bx_corr_emul = emul_corr_bx;
                stubs_comparison[1].trknmb_emul = emul_trknmb;
                GlobalPoint gp_clct_emul(getGlobalPosition(detid.rawId(), testwg, (*pe).getKeyStrip()));
                stubs_comparison[1].phi_emul = gp_clct_emul.phi();
                bookedclctV_emul[j]  = true;

                //int mingap_trig_pretrig = pretrig_trig_zone*2;
                int num_pretrig = 0;
                for (pretrig = pretrigV_emul.begin(); pretrig != pretrigV_emul.end(); pretrig++){
                  if ((*pretrig).getBX() != (*pe).getBX()) continue;
                  int hsgap = std::abs((*pretrig).getKeyStrip() - (*pe).getKeyStrip());
                  bool samechamber = true;
                  if (stat == 1 and ring == 1)
                    samechamber = (((*pretrig).getKeyStrip() <128 and (*pe).getKeyStrip() <128 ) || ((*pretrig).getKeyStrip() >= 128 and (*pe).getKeyStrip() >=128 ));
                  if (not samechamber)
                    continue;

                  if (hsgap <= pretrig_trig_zone)
                    num_pretrig ++;
                  else
                    continue;
                  if ((*pretrig).getPattern() > stubs_comparison[1].maxpattern_pretrig)
                    stubs_comparison[1].maxpattern_pretrig = (*pretrig).getPattern();
                  if ((*pretrig).getQuality() > stubs_comparison[1].maxquality_pretrig){
                    stubs_comparison[1].quality_pretrig = (*pretrig).getQuality();
                    stubs_comparison[1].key_hs_pretrig = (*pretrig).getKeyStrip();
                    stubs_comparison[1].bend_pretrig = (*pretrig).getBend();
                    stubs_comparison[1].pattern_pretrig = (*pretrig).getPattern();
                    stubs_comparison[1].bx_pretrig = (*pretrig).getBX();
                  }
                }
                stubs_comparison[1].npretrig = num_pretrig;
                break;
              }//if (data_trknmb == emul_trknmb)
            }//loop emul
            if (debug and stubs_comparison[1].key_hs_data != stubs_comparison[1].key_hs_emul)
              LogTrace("CSCTriggerPrimitivesReader")
                <<"CSCDetid "<< detid <<" not matched CLCT from data "<<  (*pd)
                <<" stubs_comparison 1 key_hs_data "<<stubs_comparison[1].key_hs_data <<" key_hs_emul "<< stubs_comparison[1].key_hs_emul << endl;
            //cout <<"CLCT data BX "<< stubs_comparison[1].bx_data <<" emul BX "<< stubs_comparison[1].bx_emul<<" emul BX corrected "<< stubs_comparison[1].bx_corr_emul << endl;
            stub_tree[1]->Fill();
          }//loop data
          for (int k=0; k<nemul; k++){
            if (clctV_emul[i].isValid() == 0) continue;
            if (bookedclctV_emul[k]) continue;
            stubs_comparison[1].init(RUN_, Event_);
            stubs_comparison[1].firstfill = firstfill;
            if (firstfill) firstfill = false;
            stubs_comparison[1].nEvents = eventsAnalyzed;
            stubs_comparison[1].endcap = endc;
            stubs_comparison[1].station = stat;
            stubs_comparison[1].ring = ring;
            stubs_comparison[1].chamber = cham;
            stubs_comparison[1].chambertype = detid.iChamberType();
            stubs_comparison[1].totStubs_data = ndata;
            stubs_comparison[1].totStubs_emul = nemul;
            stubs_comparison[1].totStubs_emul_readout = nemul_readout;
            stubs_comparison[1].trknmb_emul = clctV_emul[i].getTrknmb();
            stubs_comparison[1].nStub_data =-1;
            stubs_comparison[1].has_data = false;
            stubs_comparison[1].nStub_emul = k+1;
            stubs_comparison[1].has_emul = true;
            stubs_comparison[1].quality_emul = clctV_emul[k].getQuality();
            stubs_comparison[1].key_hs_emul = clctV_emul[k].getKeyStrip();
            stubs_comparison[1].bend_emul = clctV_emul[k].getBend();
            stubs_comparison[1].pattern_emul = clctV_emul[k].getPattern();
            stubs_comparison[1].bx_emul = clctV_emul[k].getBX();
            stubs_comparison[1].fullbx_emul = clctV_emul[k].getFullBX();
            if (clctV_emul[k].getKeyStrip()>= 128 and stat == 1 and ring == 1) testwg=5;
            else testwg = 20;
            // Emulator BX NOT Known from  the data.
            GlobalPoint gp_clct_emul(getGlobalPosition(detid.rawId(), testwg, clctV_emul[k].getKeyStrip()));
            stubs_comparison[1].phi_emul = gp_clct_emul.phi();
            bookedclctV_emul[k] = true;
            stub_tree[1]->Fill();
          }
        }
      }
    }
  } // for all chambers
} // compareCLCTs

void CSCTriggerPrimitivesReader::compareLCTs(const CSCCorrelatedLCTDigiCollection* lcts_data,
                                             const CSCCorrelatedLCTDigiCollection* lcts_emul,
                                             const CSCALCTDigiCollection* alcts_data,
                                             const CSCCLCTDigiCollection* clcts_data) {
  // Need ALCT and CLCT digi collections to convert emulator bx into
  // hardware bx.
  // Loop over all chambers in search for correlated LCTs.
  std::vector<CSCCorrelatedLCTDigi>::const_iterator pd, pe;
  for (int endc = 1; endc <= 2; endc++) {
    for (int stat = 1; stat <= 4; stat++) {
      for (int ring = 1; ring <= maxRing(stat); ring++) {
        for (int cham = 1; cham <= 36; cham++) {
          // Calculate DetId.  0th layer means whole chamber.
          CSCDetId detid(endc, stat, ring, cham, 0);

          // Skip chambers marked as bad.
          if (checkBadChambers_ && badChambers_->isInBadChamber(detid)) continue;

          std::vector<CSCCorrelatedLCTDigi> lctV_data, lctV_emul;
          std::vector<bool> bookedlctV_emul;
          const auto& drange = lcts_data->get(detid);
          for (auto digiIt = drange.first; digiIt != drange.second; digiIt++) {
            if ((*digiIt).isValid()) {
              lctV_data.push_back(*digiIt);
            }
          }

          const auto& erange = lcts_emul->get(detid);
          for (auto digiIt = erange.first; digiIt != erange.second; digiIt++) {
            if ((*digiIt).isValid()) {
              lctV_emul.push_back(*digiIt);
              bookedlctV_emul.push_back(false);
            }
          }

          int ndata = lctV_data.size();
          int nemul = lctV_emul.size();
          if (ndata == 0 && nemul == 0) continue;
          bool firstfill = true;
          int nemul_readout = 0;
          for (pe = lctV_emul.begin(); pe != lctV_emul.end(); pe++) {
            for (pd = lctV_data.begin(); pd != lctV_data.end(); pd++) {
              int bx_data = pd->getBX();
              int bx_corr = convertBXofLCT((*pe).getBX(), detid,
                                           alcts_data, clcts_data);
              if ((bx_data == bx_corr) or (abs(pe->getKeyWG() - pd->getKeyWG())<=2 ) ){
                nemul_readout++;
                break;
              }

            }
          }

          if (debug or nemul>ndata or (nemul_readout != ndata) or nemul != ndata) {

            LogTrace("CSCTriggerPrimitivesReader") << " LCTs from data "<< ndata <<" LCTs from emul "<< nemul <<" readout "<< nemul_readout;
            ostringstream strstrm;
            strstrm << "\n--- ME" << ((detid.endcap() == 1) ? "+" : "-")
                    << detid.station() << "/" << detid.ring() << "/"
                    << detid.chamber()
                    << " (sector "  << detid.triggerSector()
                    << " trig id. " << detid.triggerCscId() << "):\n";
            strstrm << "  **** " << ndata << " valid data LCTs found:\n";
            for (pd = lctV_data.begin(); pd != lctV_data.end(); pd++) {
              strstrm << "     " << (*pd);
            }
            strstrm << "\n  **** " << nemul << " valid emul LCTs found:\n";
            for (pe = lctV_emul.begin(); pe != lctV_emul.end(); pe++) {
              strstrm << "     " << (*pe);
              strstrm << "    corr BX = "
                      << convertBXofLCT((*pe).getBX(), detid,
                                        alcts_data, clcts_data);
              strstrm << " LCT pattern = " << (*pe).getPattern();
              strstrm << "\n";

            }
            LogTrace("CSCTriggerPrimitivesReader") << strstrm.str();
            //std::cout <<"CompareLCTs "<< strstrm.str()<< std::endl;
          }
          if (nemul_readout > 2){//reduce nemul_readout to 2 by hand
            //cout <<"LCT matching nemul readout is "<< nemul_readout <<", larger than 2. reduce it to 2 by hand"<<endl;
            nemul_readout = 2;
          }

          //
          //Data, add HS quality later
          int perEv_nStub_data = 0;
          int perEv_nStub_emul = 0;
          perStub[4].init(RUN_, Event_);
          perStub[5].init(RUN_, Event_);
          for (pd = lctV_data.begin(); pd != lctV_data.end(); pd++){
            perEv_nStub_data++;
          }
          perStub[4].t_endcap = endc;
          perStub[4].t_station = stat;
          perStub[4].t_chamber = cham;
          perStub[4].t_ring = ring;
          perStub[4].t_EventNumberAnalyzed = eventsAnalyzed;
          perStub[4].t_nStubs              = ndata;
          perStub[4].t_nStubs_readout              = ndata;
          event_tree[4]->Fill();
          //Emul
          for (pe = lctV_emul.begin(); pe != lctV_emul.end(); pe++){
            perEv_nStub_emul++;
          }
          perStub[5].t_endcap = endc;
          perStub[5].t_station = stat;
          perStub[5].t_chamber = cham;
          perStub[5].t_ring = ring;
          perStub[5].t_EventNumberAnalyzed = eventsAnalyzed;
          perStub[5].t_nStubs              = nemul;
          perStub[5].t_nStubs_readout      = nemul_readout;
          event_tree[5]->Fill();

          int csctype = getCSCType(detid);
          hLctCompFoundCsc[endc-1][csctype]->Fill(cham);
          int ix = chamberIX(detid);
          int ix2 = chamberIXi(detid);
          if(detid.station()>1 && detid.ring()==1) {
            hLCTCompFound2x->Fill(ix,detid.chamber()*2);
          }
          else {
            hLCTCompFound2->Fill(ix,detid.chamber());
          }
          hLCTCompFound2i->Fill(ix2,detid.chamber());
          if (ndata != nemul) {
            LogTrace("CSCTriggerPrimitivesReader")
              //cerr
              << "   +++ Different numbers of LCTs found in ME"
              << ((endc == 1) ? "+" : "-") << stat << "/"
              << ring << "/" << cham
              << ": data = " << ndata << " emulator = " << nemul << " +++\n";
          }
          else {
            hLctCompSameNCsc[endc-1][csctype]->Fill(cham);
            if(detid.station()>1 && detid.ring()==1) {
              hLCTCompSameN2x->Fill(ix,detid.chamber()*2);
            }
            else {
              hLCTCompSameN2->Fill(ix,detid.chamber());
            }
            hLCTCompSameN2i->Fill(ix2,detid.chamber());
          }
          int i =-1;
          for (pd = lctV_data.begin(); pd != lctV_data.end(); pd++) {
            i++;
            if ((*pd).isValid() == 0) continue;
            int data_trknmb    = (*pd).getTrknmb();
            int data_quality   = (*pd).getQuality();
            int data_wiregroup = (*pd).getKeyWG();
            int data_keystrip  = (*pd).getStrip();
            int data_pattern   = (*pd).getCLCTPattern();
            int data_striptype = (*pd).getStripType();
            int data_bend      = (*pd).getBend();
            int data_bx        = (*pd).getBX();

            stubs_comparison[2].init(RUN_, Event_);
            stubs_comparison[2].firstfill = firstfill;
            if (firstfill) firstfill = false;
            stubs_comparison[2].nEvents = eventsAnalyzed;
            stubs_comparison[2].endcap = endc;
            stubs_comparison[2].station = stat;
            stubs_comparison[2].ring = ring;
            stubs_comparison[2].chamber = cham;
            stubs_comparison[2].chambertype = detid.iChamberType();
            stubs_comparison[2].totStubs_data = ndata;
            stubs_comparison[2].totStubs_emul = nemul;
            stubs_comparison[2].totStubs_emul_readout = nemul_readout;
            stubs_comparison[2].nStub_data = i+1;
            stubs_comparison[2].has_data = true;
            stubs_comparison[2].quality_data = (*pd).getQuality();
            stubs_comparison[2].key_WG_data = (*pd).getKeyWG();
            stubs_comparison[2].key_hs_data = (*pd).getStrip();
            //cout <<" CSC detid "<< detid <<" LCT in data "<< (*pd) <<" do doesALCTCrossCLCT "<< endl;
            stubs_comparison[2].WGcrossHS_data = doesALCTCrossCLCT(detid,  (*pd).getKeyWG(),  (*pd).getStrip());
            stubs_comparison[2].bend_data = (*pd).getBend();
            stubs_comparison[2].pattern_data = (*pd).getCLCTPattern();
            stubs_comparison[2].bx_data = (*pd).getBX();
            stubs_comparison[2].trknmb_data = data_trknmb;
            GlobalPoint gp_lct_data(getGlobalPosition(detid.rawId(), (*pd).getKeyWG(), (*pd).getStrip()));
            stubs_comparison[2].eta_data = gp_lct_data.eta();
            stubs_comparison[2].phi_data = gp_lct_data.phi();
            int j=-1;
            for (pe = lctV_emul.begin(); pe != lctV_emul.end(); pe++) {
              j++;
              if ((*pe).isValid() == 0) continue;
              if (bookedlctV_emul[j]) continue; //used alct
              int emul_trknmb    = (*pe).getTrknmb();
              int emul_quality   = (*pe).getQuality();
              int emul_wiregroup = (*pe).getKeyWG();
              int emul_keystrip  = (*pe).getStrip();
              int emul_pattern   = (*pe).getCLCTPattern();
              int emul_striptype = (*pe).getStripType();
              int emul_bend      = (*pe).getBend();
              int emul_bx        = (*pe).getBX();
              if (abs(data_keystrip - emul_keystrip) <= 2 || abs(data_wiregroup - emul_wiregroup) <= 2){
                // Convert emulator BX into hardware BX using full 12-bit
                // BX words in ALCT and CLCT digi collections.
                int emul_corr_bx = convertBXofLCT(emul_bx, detid,
                                                  alcts_data, clcts_data);
                //std::cout <<"LCT data_bx "<< data_bx <<" emul_corr_bx "<< emul_corr_bx << std::endl;
                if (dataIsAnotherMC_)
                  emul_corr_bx = (emul_bx & 0x01);

                if (ndata == nemul) {
                  hLctCompTotalCsc[endc-1][csctype]->Fill(cham);
                  if(detid.station()>1 && detid.ring()==1) {
                    hLCTCompTotal2x->Fill(ix,detid.chamber()*2);
                  }
                  else {
                    hLCTCompTotal2->Fill(ix,detid.chamber());
                  }
                  hLCTCompTotal2i->Fill(ix2,detid.chamber());
                }
                if (data_quality   == emul_quality   &&
                    data_wiregroup == emul_wiregroup &&
                    data_keystrip  == emul_keystrip  &&
                    data_pattern   == emul_pattern   &&
                    data_striptype == emul_striptype &&
                    data_bend      == emul_bend      ){
                  if (ndata == nemul) {
                    hLctCompMatchCsc[endc-1][csctype]->Fill(cham);
                    if(detid.station()>1 && detid.ring()==1) {
                      hLCTCompMatch2x->Fill(ix,detid.chamber()*2);
                    }
                    else {
                      hLCTCompMatch2->Fill(ix,detid.chamber());
                    }
                    hLCTCompMatch2i->Fill(ix2,detid.chamber());
                  }
                  if (debug) LogTrace("CSCTriggerPrimitivesReader")
                               << "       Identical LCTs #" << data_trknmb;
                }
                else {
                  LogTrace("CSCTriggerPrimitivesReader")
                    << "       Different LCTs #" << data_trknmb << " in ME"
                    << ((endc == 1) ? "+" : "-") << stat << "/"
                    << ring << "/" << cham <<" data_lct bx "<< data_bx<< " emul bx "<< emul_corr_bx;
                }
                stubs_comparison[2].nStub_emul = j+1;
                stubs_comparison[2].has_emul = true;
                stubs_comparison[2].quality_emul = (*pe).getQuality();
                stubs_comparison[2].key_WG_emul = (*pe).getKeyWG();
                stubs_comparison[2].key_hs_emul = (*pe).getStrip();
                stubs_comparison[2].bend_emul = (*pe).getBend();
                stubs_comparison[2].pattern_emul = (*pe).getCLCTPattern();
                stubs_comparison[2].bx_emul = (*pe).getBX();
                //cout <<" CSC detid "<< detid <<" LCT in emul "<< (*pe) <<" do doesALCTCrossCLCT "<< endl;
                stubs_comparison[2].WGcrossHS_emul = doesALCTCrossCLCT(detid,  (*pe).getKeyWG(),  (*pe).getStrip());
                stubs_comparison[2].bx_corr_emul = emul_corr_bx;
                stubs_comparison[2].trknmb_emul = emul_trknmb;
                GlobalPoint gp_lct_emul(getGlobalPosition(detid.rawId(), (*pe).getKeyWG(), (*pe).getStrip()));
                stubs_comparison[2].eta_emul = gp_lct_emul.eta();
                stubs_comparison[2].phi_emul = gp_lct_emul.phi();
                bookedlctV_emul[j]  = true;
                break;
              }//if (data_trknmb == emul_trknmb)
            }//loop emul
            if (debug and stubs_comparison[2].key_hs_data != stubs_comparison[2].key_hs_emul)
              LogTrace("CSCTriggerPrimitivesReader")
                <<"CSCDetid "<< detid <<" not matched LCT from Data "<< (*pd)
                <<" stubs_comparison 2 key_hs_data "<<stubs_comparison[2].key_hs_data <<" key_hs_emul "<< stubs_comparison[2].key_hs_emul << endl;
            stub_tree[2]->Fill();
          }//loop data
          for (int k=0; k<nemul; k++){
            if (bookedlctV_emul[k]) continue;
            if (lctV_emul[i].isValid() == 0) continue;
            stubs_comparison[2].init(RUN_, Event_);
            stubs_comparison[2].nEvents = eventsAnalyzed;
            stubs_comparison[2].endcap = endc;
            stubs_comparison[2].station = stat;
            stubs_comparison[2].ring = ring;
            stubs_comparison[2].chamber = cham;
            stubs_comparison[2].chambertype = detid.iChamberType();
            stubs_comparison[2].totStubs_data = ndata;
            stubs_comparison[2].totStubs_emul = nemul;
            stubs_comparison[2].totStubs_emul_readout = nemul_readout;
            stubs_comparison[2].trknmb_emul = lctV_emul[i].getTrknmb();
            stubs_comparison[2].nStub_data =-1;
            stubs_comparison[2].has_data = false;
            stubs_comparison[2].nStub_emul = k+1;
            stubs_comparison[2].has_emul = true;
            stubs_comparison[2].WGcrossHS_emul = doesALCTCrossCLCT(detid,  lctV_emul[k].getKeyWG(),  lctV_emul[k].getStrip());
            stubs_comparison[2].quality_emul = lctV_emul[k].getQuality();
            stubs_comparison[2].key_WG_emul = lctV_emul[k].getKeyWG();
            stubs_comparison[2].key_hs_emul = lctV_emul[k].getStrip();
            stubs_comparison[2].bend_emul = lctV_emul[k].getBend();
            stubs_comparison[2].pattern_emul = lctV_emul[k].getCLCTPattern();
            stubs_comparison[2].bx_emul = lctV_emul[k].getBX();
            //stubs_comparison[2].fullbx_emul = lctV_emul[k].getFullBX();
            // Emulator BX NOT Known from  the data.
            GlobalPoint gp_lct_emul(getGlobalPosition(detid.rawId(), lctV_emul[k].getKeyWG(), lctV_emul[k].getStrip()));
            stubs_comparison[2].eta_emul = gp_lct_emul.eta();
            stubs_comparison[2].phi_emul = gp_lct_emul.phi();
            bookedlctV_emul[k] = true;
            stub_tree[2]->Fill();

          }
        }
      }
    }
  }
}


void CSCTriggerPrimitivesReader::compareMPCLCTs(
                                                const CSCCorrelatedLCTDigiCollection* lcts_data,
                                                const CSCCorrelatedLCTDigiCollection* lcts_emul,
                                                const CSCALCTDigiCollection* alcts_data,
                                                const CSCCLCTDigiCollection* clcts_data) {

  // Need ALCT and CLCT digi collections to convert emulator bx into
  // hardware bx.
  // Loop over all chambers in search for correlated LCTs.
  std::vector<CSCCorrelatedLCTDigi>::const_iterator pd, pe;
  for (int endc = 1; endc <= 2; endc++) {
    for (int stat = 1; stat <= 4; stat++) {
      for (int ring = 1; ring <= maxRing(stat); ring++) {
        for (int cham = 1; cham <= 36; cham++) {
          // Calculate DetId.  0th layer means whole chamber.
          CSCDetId detid(endc, stat, ring, cham, 0);

          // Skip chambers marked as bad.
          if (checkBadChambers_ && badChambers_->isInBadChamber(detid)) continue;

          std::vector<CSCCorrelatedLCTDigi> lctV_data, lctV_emul;
          std::vector<bool> bookedlctV_emul;
          const auto& drange = lcts_data->get(detid);
          for (auto digiIt = drange.first; digiIt != drange.second; digiIt++) {
            if ((*digiIt).isValid()) {
              lctV_data.push_back(*digiIt);
            }
          }

          const auto& erange = lcts_emul->get(detid);
          for (auto digiIt = erange.first; digiIt != erange.second; digiIt++) {
            if ((*digiIt).isValid()) {
              lctV_emul.push_back(*digiIt);
              bookedlctV_emul.push_back(false);
            }
          }

          int ndata = lctV_data.size();
          int nemul = lctV_emul.size();
          if (ndata == 0 && nemul == 0) continue;
          bool firstfill = true;

          if (debug) {
            ostringstream strstrm;
            strstrm << "\n--- ME" << ((detid.endcap() == 1) ? "+" : "-")
                    << detid.station() << "/" << detid.ring() << "/"
                    << detid.chamber()
                    << " (sector "  << detid.triggerSector()
                    << " trig id. " << detid.triggerCscId() << "):\n";
            strstrm << "  **** " << ndata << " valid data LCTs found:\n";
            for (pd = lctV_data.begin(); pd != lctV_data.end(); pd++) {
              strstrm << "     " << (*pd);
            }
            strstrm << "\n  **** " << nemul << " valid emul LCTs found:\n";
            for (pe = lctV_emul.begin(); pe != lctV_emul.end(); pe++) {
              strstrm << "     " << (*pe);
              strstrm << "    corr BX = "
                      << convertBXofLCT((*pe).getBX(), detid,
                                        alcts_data, clcts_data);
              strstrm << " LCT pattern = " << (*pe).getPattern();
              strstrm << "\n";
            }
            LogTrace("CSCTriggerPrimitivesReader") << strstrm.str();
            //std::cout <<"CompareLCTs "<< strstrm.str()<< std::endl;
          }

          //int ix = chamberIX(detid);
          //int ix2 = chamberIXi(detid);
          if (ndata != nemul) {
            LogTrace("CSCTriggerPrimitivesReader")
              << "   +++ Different numbers of MPC LCTs found in ME"
              << ((endc == 1) ? "+" : "-") << stat << "/"
              << ring << "/" << cham
              << ": data = " << ndata << " emulator = " << nemul << " +++\n";
          }
          /*else {
            hLctCompSameNCsc[endc-1][csctype]->Fill(cham);
            if(detid.station()>1 && detid.ring()==1) {
            hLCTCompSameN2x->Fill(ix,detid.chamber()*2);
            }
            else {
            hLCTCompSameN2->Fill(ix,detid.chamber());
            }
            hLCTCompSameN2i->Fill(ix2,detid.chamber());
            }*/
          int i =-1;
          for (pd = lctV_data.begin(); pd != lctV_data.end(); pd++) {
            i++;
            if ((*pd).isValid() == 0) continue;
            int data_trknmb    = (*pd).getTrknmb();
            int data_quality   = (*pd).getQuality();
            int data_wiregroup = (*pd).getKeyWG();
            int data_keystrip  = (*pd).getStrip();
            int data_pattern   = (*pd).getCLCTPattern();
            int data_striptype = (*pd).getStripType();
            int data_bend      = (*pd).getBend();
            int data_bx        = (*pd).getBX();

            stubs_comparison[3].init(RUN_, Event_);
            stubs_comparison[3].firstfill = firstfill;
            if (firstfill) firstfill = false;
            stubs_comparison[3].nEvents = eventsAnalyzed;
            stubs_comparison[3].endcap = endc;
            stubs_comparison[3].station = stat;
            stubs_comparison[3].ring = ring;
            stubs_comparison[3].chamber = cham;
            stubs_comparison[3].chambertype = detid.iChamberType();
            stubs_comparison[3].totStubs_data = ndata;
            stubs_comparison[3].totStubs_emul = nemul;
            stubs_comparison[3].nStub_data = i+1;
            stubs_comparison[3].has_data = true;
            stubs_comparison[3].WGcrossHS_data = doesALCTCrossCLCT(detid,  (*pd).getKeyWG(),  (*pd).getStrip());
            stubs_comparison[3].quality_data = (*pd).getQuality();
            stubs_comparison[3].key_WG_data = (*pd).getKeyWG();
            stubs_comparison[3].key_hs_data = (*pd).getStrip();
            stubs_comparison[3].bend_data = (*pd).getBend();
            stubs_comparison[3].pattern_data = (*pd).getCLCTPattern();
            stubs_comparison[3].bx_data = (*pd).getBX();
            stubs_comparison[3].trknmb_data = data_trknmb;
            GlobalPoint gp_lct_data(getGlobalPosition(detid.rawId(), (*pd).getKeyWG(), (*pd).getStrip()));
            stubs_comparison[3].eta_data = gp_lct_data.eta();
            stubs_comparison[3].phi_data = gp_lct_data.phi();
            int j=-1;
            for (pe = lctV_emul.begin(); pe != lctV_emul.end(); pe++) {
              j++;
              if ((*pe).isValid() == 0) continue;
              int emul_trknmb    = (*pe).getTrknmb();
              int emul_quality   = (*pe).getQuality();
              int emul_wiregroup = (*pe).getKeyWG();
              int emul_keystrip  = (*pe).getStrip();
              int emul_pattern   = (*pe).getCLCTPattern();
              int emul_striptype = (*pe).getStripType();
              int emul_bend      = (*pe).getBend();
              int emul_bx        = (*pe).getBX();
              if (data_trknmb == emul_trknmb) {
                // Convert emulator BX into hardware BX using full 12-bit
                // BX words in ALCT and CLCT digi collections.
                int emul_corr_bx = convertBXofLCT(emul_bx, detid,
                                                  alcts_data, clcts_data);
                //std::cout <<"LCT data_bx "<< data_bx <<" emul_corr_bx "<< emul_corr_bx << std::endl;
                if (dataIsAnotherMC_)
                  emul_corr_bx = (emul_bx & 0x01);

                /*if (ndata == nemul) {
                  hLctCompTotalCsc[endc-1][csctype]->Fill(cham);
                  if(detid.station()>1 && detid.ring()==1) {
                  hLCTCompTotal2x->Fill(ix,detid.chamber()*2);
                  }
                  else {
                  hLCTCompTotal2->Fill(ix,detid.chamber());
                  }
                  hLCTCompTotal2i->Fill(ix2,detid.chamber());
                  }*/
                if (data_quality   == emul_quality   &&
                    data_wiregroup == emul_wiregroup &&
                    data_keystrip  == emul_keystrip  &&
                    data_pattern   == emul_pattern   &&
                    data_striptype == emul_striptype &&
                    data_bend      == emul_bend      ){
                  //data_bx        == emul_corr_bx) {
                  /*if (ndata == nemul) {
                    hLctCompMatchCsc[endc-1][csctype]->Fill(cham);
                    if(detid.station()>1 && detid.ring()==1) {
                    hLCTCompMatch2x->Fill(ix,detid.chamber()*2);
                    }
                    else {
                    hLCTCompMatch2->Fill(ix,detid.chamber());
                    }
                    hLCTCompMatch2i->Fill(ix2,detid.chamber());
                    }*/
                  if (debug) LogTrace("CSCTriggerPrimitivesReader")
                               << "       Identical LCTs #" << data_trknmb;
                  stubs_comparison[3].nStub_emul = j+1;
                  stubs_comparison[3].has_emul = true;
                  stubs_comparison[3].WGcrossHS_emul = doesALCTCrossCLCT(detid,  (*pe).getKeyWG(),  (*pe).getStrip());
                  stubs_comparison[3].quality_emul = (*pe).getQuality();
                  stubs_comparison[3].key_WG_emul = (*pe).getKeyWG();
                  stubs_comparison[3].key_hs_emul = (*pe).getStrip();
                  stubs_comparison[3].bend_emul = (*pe).getBend();
                  stubs_comparison[3].pattern_emul = (*pe).getCLCTPattern();
                  stubs_comparison[3].bx_emul = (*pe).getBX();
                  stubs_comparison[3].bx_corr_emul = emul_corr_bx;
                  stubs_comparison[3].trknmb_emul = emul_trknmb;
                  GlobalPoint gp_lct_emul(getGlobalPosition(detid.rawId(), (*pe).getKeyWG(), (*pe).getStrip()));
                  stubs_comparison[3].eta_emul = gp_lct_emul.eta();
                  stubs_comparison[3].phi_emul = gp_lct_emul.phi();
                  bookedlctV_emul[j]  = true;
                  break;
                }
                else {
                  LogTrace("CSCTriggerPrimitivesReader")
                    //cerr
                    << "       Different LCTs #" << data_trknmb << " in ME"
                    << ((endc == 1) ? "+" : "-") << stat << "/"
                    << ring << "/" << cham <<" data_lct bx "<< data_bx<< " emul bx "<< emul_corr_bx;
                }
              }
            }//loop emul
            if (debug and stubs_comparison[3].key_hs_data != stubs_comparison[3].key_hs_emul)
              LogTrace("CSCTriggerPrimitivesReader")
                <<"stubs_comparison 2 key_hs_data "<<stubs_comparison[3].key_hs_data <<" key_hs_emul "<< stubs_comparison[3].key_hs_emul;
            stub_tree[3]->Fill();
          }//loop data

          for (int k=0; k<nemul; k++){
            if (bookedlctV_emul[k]) continue;
            stubs_comparison[3].init(RUN_, Event_);
            stubs_comparison[3].firstfill = firstfill;
            if (firstfill) firstfill = false;
            stubs_comparison[3].nEvents = eventsAnalyzed;
            stubs_comparison[3].endcap = endc;
            stubs_comparison[3].station = stat;
            stubs_comparison[3].ring = ring;
            stubs_comparison[3].chamber = cham;
            stubs_comparison[3].chambertype = detid.iChamberType();
            stubs_comparison[3].totStubs_data = ndata;
            stubs_comparison[3].totStubs_emul = nemul;
            stubs_comparison[3].trknmb_emul = lctV_emul[i].getTrknmb();
            stubs_comparison[3].nStub_data =-1;
            stubs_comparison[3].has_data = false;
            stubs_comparison[3].nStub_emul = k+1;
            stubs_comparison[3].has_emul = true;
            stubs_comparison[3].WGcrossHS_emul = doesALCTCrossCLCT(detid, lctV_emul[k].getKeyWG(),  lctV_emul[k].getStrip());
            stubs_comparison[3].quality_emul = lctV_emul[k].getQuality();
            stubs_comparison[3].key_WG_emul = lctV_emul[k].getKeyWG();
            stubs_comparison[3].key_hs_emul = lctV_emul[k].getStrip();
            stubs_comparison[3].bend_emul = lctV_emul[k].getBend();
            stubs_comparison[3].pattern_emul = lctV_emul[k].getCLCTPattern();
            stubs_comparison[3].bx_emul = lctV_emul[k].getBX();
            //stubs_comparison[2].fullbx_emul = lctV_emul[k].getFullBX();
            // Emulator BX NOT Known from  the data.
            GlobalPoint gp_lct_emul(getGlobalPosition(detid.rawId(), lctV_emul[k].getKeyWG(), lctV_emul[k].getStrip()));
            stubs_comparison[3].eta_emul = gp_lct_emul.eta();
            stubs_comparison[3].phi_emul = gp_lct_emul.phi();
            bookedlctV_emul[k] = true;
            stub_tree[3]->Fill();

          }
        }//end loop of chambers
      }
    }
  }
}




int CSCTriggerPrimitivesReader::convertBXofLCT(
                                               const int emul_bx, const CSCDetId& detid,
                                               const CSCALCTDigiCollection* alcts_data,
                                               const CSCCLCTDigiCollection* clcts_data) {
  int full_anode_bx = -999;
  //int full_cathode_bx = -999;
  int lct_bx = -999;
  int tbin_anode_offset = 5; // 2007, run 14419.

  // Extract full 12-bit anode BX word from ALCT collections.
  const auto& arange = alcts_data->get(detid);
  for (auto digiIt = arange.first; digiIt != arange.second; digiIt++) {
    if ((*digiIt).isValid()) {
      full_anode_bx = (*digiIt).getFullBX();
      break;
    }
  }

  // Extract full 12-bit cathode BX word from CLCT collections.
  const auto& crange = clcts_data->get(detid);
  for (auto digiIt = crange.first; digiIt != crange.second; digiIt++) {
    if ((*digiIt).isValid()) {
      //full_cathode_bx = (*digiIt).getFullBX();
      break;
    }
  }

  // Use these 12-bit BX's to convert emulator BX into hardware BX.
  if (full_anode_bx == -999) {
    // What to do???
    edm::LogWarning("L1CSCTPEmulatorWrongInput")
      << "+++ Warning in convertBXofLCT(): full anode BX is not available!"
      << " +++\n";
  }
  else {
    // LCT BX has two bits: the least-significant bit is the LSB of ALCT BX;
    // the MSB is 1/0 depending on whether the 12-bit full cathode BX is 0
    // or not.
    lct_bx = (full_anode_bx + emul_bx - tbin_anode_offset) & 0x01;
    // SV, 12/Jun/08: it looks like this bit is never set - docu must be
    // wrong.
    //lct_bx = lct_bx | ((full_cathode_bx == 0) << 1);
  }
  //std::cout <<"convertBXofLCT old emul_bx "<< emul_bx <<" new lct bx "<< lct_bx <<" ful_anode_bx "<< full_anode_bx << std::endl;
  return lct_bx;
}


void CSCTriggerPrimitivesReader::HotWires(const edm::Event& iEvent) {
  if (!bookedHotWireHistos) bookHotWireHistos();
  edm::Handle<CSCWireDigiCollection> wires;
  //  iEvent.getByLabel(wireDigiProducer_.label(), wireDigiProducer_.instance(), wires);
  iEvent.getByToken(wireDigi_token_, wires);

  int serial_old=-1;
  for (auto dWDiter=wires->begin(); dWDiter!=wires->end(); dWDiter++) {
    CSCDetId id = (CSCDetId)(*dWDiter).first;
    int serial = chamberSerial(id)-1;
    //     printf("serial %i\n",serial);
    std::vector<CSCWireDigi>::const_iterator wireIter = (*dWDiter).second.first;
    std::vector<CSCWireDigi>::const_iterator lWire = (*dWDiter).second.second;
    bool has_layer=false;
    for( ; wireIter != lWire; ++wireIter) {
      has_layer=true;
      int i_layer= id.layer()-1;
      int i_wire = wireIter->getWireGroup()-1;
      int nbins = wireIter->getTimeBinsOn().size();
      int serial2=serial*(6*112)+i_layer*112+i_wire;
      /*
        printf("endcap %i, station %i, ring %i, chamber %i, serial %i, serial2 %i, layer %i, wg %i\n",
        id.endcap(), id.station(), id.ring(), id.chamber(), serial, serial2, i_layer, i_wire);
      */
      hHotWire1->Fill(serial2,nbins);
      /*
        if(id.station()==1) {
        printf("ring %i, chamber type %i, wg %i\n",id.ring(),id.iChamberType(),i_wire);
        }
      */
    }
    if(serial_old!=serial && has_layer) {
      //       nHotCham[serial]++;
      hHotCham1->Fill(serial);
      //    printf("serial %i filled\n",serial);
      serial_old=serial;
    }
  }
}

void CSCTriggerPrimitivesReader::MCStudies(const edm::Event& ev,
                                           const CSCALCTDigiCollection* alcts,
                                           const CSCCLCTDigiCollection* clcts) {
  // MC particles, if any.
  //edm::Handle<edm::HepMCProduct> mcp;
  //ev.getByLabel("source", mcp);
  //ev.getByType(mcp);
  vector<edm::Handle<edm::HepMCProduct> > allhepmcp;
  // Use "getManyByType" to be able to check the existence of MC info.
  ev.getManyByType(allhepmcp);

  //cout << "HepMC info: " << allhepmcp.size() << endl;
  if (allhepmcp.size() > 0) {
    const HepMC::GenEvent& mc = allhepmcp[0]->getHepMCData();
    int i = 0;
    for (HepMC::GenEvent::particle_const_iterator p = mc.particles_begin();
         p != mc.particles_end(); ++p) {
      int id = (*p)->pdg_id();
      double phitmp = (*p)->momentum().phi();
      if (phitmp < 0) phitmp += 2.*M_PI;
      if (debug) LogDebug("CSCTriggerPrimitivesReader")
                   << "MC part #" << ++i << ": id = "  << id
                   << ", status = " << (*p)->status()
                   << "\n   pX = " << (*p)->momentum().x()
                   << ", pY = " << (*p)->momentum().y()
                   << ", pT = " << (*p)->momentum().perp() << " GeV"
                   << ", p =  " << (*p)->momentum().rho()  << " GeV"
                   << "\n   eta = " << (*p)->momentum().pseudoRapidity()
                   << ", phi = " << phitmp << " (" << phitmp*180./M_PI << " deg)";
    }

    // If hepMC info is there, try to get wire and comparator digis,
    // and SimHits.
    edm::Handle<CSCWireDigiCollection>       wireDigis;
    edm::Handle<CSCComparatorDigiCollection> compDigis;
    edm::Handle<edm::PSimHitContainer>       simHitsH;
    //    ev.getByLabel(wireDigiProducer_.label(), wireDigiProducer_.instance(),
    //        wireDigis);
    //    ev.getByLabel(compDigiProducer_.label(), compDigiProducer_.instance(),
    //        compDigis);
    //    ev.getByLabel(simHitProducer_.label(), simHitProducer_.instance(),
    //        simHits);
    ev.getByToken(wireDigi_token_, wireDigis);
    ev.getByToken(compDigi_token_, compDigis);
    ev.getByToken(simHit_token_, simHitsH);
    const edm::PSimHitContainer& simHits = *simHitsH.product();

    if (!wireDigis.isValid()) {
      edm::LogWarning("L1CSCTPEmulatorWrongInput")
        << "+++ Warning: Collection of wire digis with label"
        << wireDigiProducer_.label()
        << " requested, but not found in the event... Skipping the rest +++\n";
      return;
    }
    if (!compDigis.isValid()) {
      edm::LogWarning("L1CSCTPEmulatorWrongInput")
        << "+++ Warning: Collection of comparator digis with label"
        << compDigiProducer_.label()
        << " requested, but not found in the event... Skipping the rest +++\n";
      return;
    }
    if (!simHitsH.isValid()) {
      edm::LogWarning("L1CSCTPEmulatorWrongInput")
        << "+++ Warning: Collection of SimHits with label"
        << simHitProducer_.label()
        << " requested, but not found in the event... Skipping the rest +++\n";
      return;
    }


    if (debug) LogTrace("CSCTriggerPrimitivesReader")
                 << "   #CSC SimHits: " << simHits.size();

    // MC-based resolution studies.
    calcResolution(alcts, clcts, wireDigis.product(), compDigis.product(),
                   simHits);

    // MC-based efficiency studies.
    calcEfficiency(alcts, clcts, simHits);
  }
}

void CSCTriggerPrimitivesReader::calcResolution(const CSCALCTDigiCollection* alcts, const CSCCLCTDigiCollection* clcts,
                                                const CSCWireDigiCollection* wiredc,
                                                const CSCComparatorDigiCollection* compdc,
                                                const edm::PSimHitContainer& allSimHits) {

  // Book histos when called for the first time.
  if (!bookedResolHistos) bookResolHistos();

  // ALCT resolution
  CSCAnodeLCTAnalyzer alct_analyzer;
  alct_analyzer.setGeometry(geom_);

  for (auto adetUnitIt = alcts->begin(); adetUnitIt != alcts->end(); adetUnitIt++) {
    const CSCDetId& id = (*adetUnitIt).first;
    if (checkBadChambers_ && badChambers_->isInBadChamber(id)) continue;
    const auto& range = (*adetUnitIt).second;
    for (auto digiIt = range.first;
         digiIt != range.second; digiIt++) {

      bool alct_valid = (*digiIt).isValid();
      if (alct_valid) {
        vector<CSCAnodeLayerInfo> alctInfo =
          alct_analyzer.getSimInfo(*digiIt, id, wiredc, &allSimHits);

        double hitPhi = -999.0, hitEta = -999.0;
        int hitWG = alct_analyzer.nearestWG(alctInfo, hitPhi, hitEta);
        if (hitWG >= 0.) {
          // Key wire group and key layer id.
          int wiregroup = (*digiIt).getKeyWG();

          CSCDetId layerId(id.endcap(), id.station(), id.ring(),
                           id.chamber(), CSCConstants::KEY_ALCT_LAYER);
          int endc    = id.endcap();
          int stat    = id.station();
          int csctype = getCSCType(id);

          double alctEta  = alct_analyzer.getWGEta(layerId, wiregroup);
          double deltaEta = alctEta - hitEta;
          hResolDeltaEta->Fill(deltaEta);

          double deltaWG = wiregroup - hitWG;
          if (debug) LogTrace("CSCTriggerPrimitivesReader")
                       << "WG: MC = " << hitWG << " rec = " << wiregroup
                       << " delta = " << deltaWG;
          hResolDeltaWG->Fill(deltaWG);

          hEtaRecVsSim->Fill(fabs(hitEta), fabs(alctEta));
          hEtaDiffCsc[csctype][0]->Fill(deltaEta);
          hEtaDiffCsc[csctype][endc]->Fill(deltaEta);
          hAlctVsEta[stat-1]->Fill(fabs(alctEta));
          hEtaDiffVsEta[stat-1]->Fill(fabs(alctEta), fabs(deltaEta));
          hEtaDiffVsWireCsc[csctype]->Fill(wiregroup, deltaEta);
        }
        // should I comment out this "else"?
        //else {
        //  edm::LogWarning("L1CSCTPEmulatorWrongInput")
        //    << "+++ Warning in calcResolution(): no matched SimHit"
        //    << " found! +++\n";
        //}
      }
    }
  }

  // CLCT resolution
  static const int key_layer = CSCConstants::KEY_CLCT_LAYER;
  CSCCathodeLCTAnalyzer clct_analyzer;
  clct_analyzer.setGeometry(geom_);

  for (auto cdetUnitIt = clcts->begin(); cdetUnitIt != clcts->end(); cdetUnitIt++) {
    const CSCDetId& id = (*cdetUnitIt).first;
    if (checkBadChambers_ && badChambers_->isInBadChamber(id)) continue;
    const auto& range = (*cdetUnitIt).second;
    for (auto digiIt = range.first;
         digiIt != range.second; digiIt++) {

      bool clct_valid = (*digiIt).isValid();
      if (clct_valid) {
        vector<CSCCathodeLayerInfo> clctInfo =
          clct_analyzer.getSimInfo(*digiIt, id, compdc, &allSimHits);

        double hitPhi = -999.0, hitEta = -999.0, deltaStrip = -999.0;
        int hitHS = clct_analyzer.nearestHS(clctInfo, hitPhi, hitEta);
        if (hitHS >= 0.) {
          // Key strip and key layer id.
          int halfstrip = (*digiIt).getKeyStrip();
          int strip     = halfstrip/2;
          int distrip   = halfstrip/4;
          int stripType = (*digiIt).getStripType();

          CSCDetId layerId(id.endcap(), id.station(), id.ring(),
                           id.chamber(), key_layer);
          int endc    = id.endcap();
          int stat    = id.station();
          int csctype = getCSCType(id);

          // 'float strip' is in the units of 'strip', i.e., angular
          // widths of each strip. The counting is from 0.0 at the extreme
          // edge of the 'first' strip at one edge of the detector.
          float fstrip = -999.;
          if (stripType == 0) { // di-strip CLCT
            fstrip = strip + 1.;
          }
          else {                // half-strip CLCT
            fstrip = strip + 0.5*(halfstrip%2) + 0.25;
          }
          double clctPhi = clct_analyzer.getStripPhi(layerId, fstrip);
          double deltaPhi = clctPhi - hitPhi;
          if      (deltaPhi < -M_PI) deltaPhi += 2.*M_PI;
          else if (deltaPhi >  M_PI) deltaPhi -= 2.*M_PI;
          deltaPhi *= 1000; // in mrad
          if      (hitPhi  < 0) hitPhi  += 2.*M_PI;
          if      (clctPhi < 0) clctPhi += 2.*M_PI;
          if (debug) LogTrace("CSCTriggerPrimitivesReader")
                       << " clctPhi = " << clctPhi << " hitPhi = " << hitPhi
                       << " deltaPhi = " << deltaPhi;

          hResolDeltaPhi->Fill(deltaPhi);
          if (stripType == 0) { // di-strip CLCT
            deltaStrip = distrip - hitHS/4;
            hResolDeltaDS->Fill(deltaStrip);
            hResolDeltaPhiDS->Fill(deltaPhi);
            hPhiDiffVsStripCsc[csctype][0]->Fill(distrip,   deltaPhi);
          }
          else {                // half-strip CLCT
            deltaStrip = halfstrip - hitHS;
            hResolDeltaHS->Fill(deltaStrip);
            hResolDeltaPhiHS->Fill(deltaPhi);
            hPhiDiffVsStripCsc[csctype][1]->Fill(halfstrip, deltaPhi);
          }
          if (debug) LogTrace("CSCTriggerPrimitivesReader")
                       << "Half-strip: MC = " << hitHS << " rec = " << halfstrip
                       << " pattern type = " << stripType << " delta = " << deltaStrip;

          hPhiRecVsSim->Fill(hitPhi, clctPhi);
          hPhiDiffCsc[csctype][0]->Fill(deltaPhi);
          hPhiDiffCsc[csctype][endc]->Fill(deltaPhi);
          hPhiDiffCsc[csctype][stripType+3]->Fill(deltaPhi);
          hClctVsPhi[stat-1]->Fill(clctPhi);
          hPhiDiffVsPhi[stat-1]->Fill(clctPhi, fabs(deltaPhi));

          // Histograms to check phi offsets for various pattern types
          if (stripType == 1) { // half-strips
            double hsperrad = getHsPerRad(csctype); // halfstrips-per-radian
            if((endc == 1 && (stat == 1 || stat == 2)) ||
               (endc == 2 && (stat == 3 || stat == 4))) {
              hPhiDiffPattern[(*digiIt).getPattern()]->Fill(deltaPhi/1000*hsperrad);
            }
          }
        }
        // should I comment out this "else"?
        //else {
        //  edm::LogWarning("L1CSCTPEmulatorWrongInput")
        //    << "+++ Warning in calcResolution(): no matched SimHit"
        //    << " found! +++\n";
        //}

        // "True bend", defined as difference in phi between muon hit
        // positions in the first and in the sixth layers.
        double phi1 = -999.0, phi6 = -999.0;
        vector<CSCCathodeLayerInfo>::const_iterator pli;
        for (pli = clctInfo.begin(); pli != clctInfo.end(); pli++) {
          CSCDetId layerId = pli->getId();
          int layer = layerId.layer();
          if (layer == 1 || layer == 6) {
            // Get simHits in this layer.
            for (const auto& psh : allSimHits) {
              // Find detId where simHit is located.
              CSCDetId hitId = (CSCDetId)psh.detUnitId();
              if (hitId == layerId &&
                  abs(psh.particleType()) == 13) { // muon hits only
                const CSCLayer* csclayer = geom_->layer(layerId);
                const auto& thisPoint = csclayer->toGlobal(psh.localPosition());
                double phi = thisPoint.phi();
                if (layer == 1)      phi1 = phi;
                else if (layer == 6) phi6 = phi;
                break; // simply take the first suitable hit.
              }
            }
          }
        }
        if (phi1 > -99. && phi6 > -99.) {
          double deltaPhi = phi1 - phi6;
          if (deltaPhi > M_PI)       deltaPhi -= 2.*M_PI;
          else if (deltaPhi < -M_PI) deltaPhi += 2.*M_PI;
          int csctype = getCSCType(id);
          hTrueBendCsc[csctype]->Fill(deltaPhi*1000.); // in mrad
        }
      }
    }
  }
}

void CSCTriggerPrimitivesReader::calcEfficiency(const CSCALCTDigiCollection* alcts,
                                                const CSCCLCTDigiCollection* clcts,
                                                const edm::PSimHitContainer& allSimHits) {

  // Book histos when called for the first time.
  if (!bookedEfficHistos) bookEfficHistos();

  // Create list of chambers having SimHits.
  vector<CSCDetId> chamberIds;
  vector<CSCDetId>::const_iterator chamberIdIt;
  for (const auto& simHitIt : allSimHits) {
    // Find detId where simHit is located.
    bool sameId = false;
    CSCDetId hitId = (CSCDetId) simHitIt.detUnitId();
    // Skip chambers marked as bad (includes most of ME4/2 chambers).
    if (checkBadChambers_ && badChambers_->isInBadChamber(hitId)) continue;
    //if (hitId.ring() == 4) continue; // skip ME1/A for now.
    if (!plotME1A && hitId.ring() == 4) continue;
    for (const auto& chamberId : chamberIds) {
      if (chamberId.endcap()  == hitId.endcap() &&
          chamberId.station() == hitId.station() &&
          chamberId.ring()    == hitId.ring() &&
          chamberId.chamber() == hitId.chamber()) {
        sameId = true;
        break;
      }
    }
    if (!sameId) {
      CSCDetId newChamberId(hitId.endcap(), hitId.station(), hitId.ring(),
                            hitId.chamber(), 0);
      chamberIds.push_back(newChamberId);
    }
  }
  LogTrace("CSCTriggerPrimitivesReader")
    << "Found SimHits in " << chamberIds.size() << " CSCs";

  bool used[CSCConstants::NUM_LAYERS];
  vector<PSimHit> simHitsV[CSCConstants::NUM_LAYERS];
  for (const auto& chamberId : chamberIds) {
    // Find out how many layers of this chamber have SimHits.
    int nLayers = 0;
    for (int ilayer = 0; ilayer < CSCConstants::NUM_LAYERS; ilayer++) {
      used[ilayer] = false;
      simHitsV[ilayer].clear();
    }

    int endcap  = chamberId.endcap();
    int station = chamberId.station();
    int ring    = chamberId.ring();
    int chamber = chamberId.chamber();
    for (const auto& simHitIt : allSimHits) {
      CSCDetId hitId = (CSCDetId) simHitIt.detUnitId();
      if (hitId.endcap() == endcap && hitId.station() == station &&
          hitId.ring()   == ring   && hitId.chamber() == chamber) {
        int layer = hitId.layer() - 1;
        if (!used[layer] && abs(simHitIt.particleType()) == 13) {
          nLayers++;
          used[layer] = true;
          simHitsV[layer].push_back(simHitIt);
        }
      }
    }
    LogTrace("CSCTriggerPrimitivesReader")
      << "CSC in ME" << ((endcap == 1) ? "+" : "-")
      << station << "/" << ring << "/" << chamber
      << " has muon hits in " << nLayers << " layers";

    // If the number of layers with hits is above threshold, look for
    // a presence of LCTs.
    if (nLayers > 3) { // Should be a parameter.
      // Start with the key layer and take the eta of the first hit.
      // Really crude; should be improved.
      double hitEta = -999.;
      for (int ilayer = 2; ilayer < CSCConstants::NUM_LAYERS; ilayer++) {
        vector<PSimHit> layerSimHitsV = simHitsV[ilayer];
        if (layerSimHitsV.size() > 0) {
          LocalPoint hitLP = layerSimHitsV[0].localPosition();
          CSCDetId layerId = (CSCDetId)(layerSimHitsV[0]).detUnitId();
          const CSCLayer* csclayer = geom_->layer(layerId);
          GlobalPoint hitGP = csclayer->toGlobal(hitLP);
          hitEta = hitGP.eta();
          break;
        }
      }
      if (hitEta < -3.) {
        edm::LogWarning("L1CSCTPEmulatorWrongValues")
          << "+++ Warning in calcEfficiency(): no SimHit found"
          << " where there must be at least " << nLayers << "! +++\n";
        continue;
      }
      int csctype = getCSCType(*chamberIdIt);
      hEfficHitsEta[station-1]->Fill(fabs(hitEta));
      hEfficHitsEtaCsc[csctype]->Fill(fabs(hitEta));

      bool isALCT = false;
      for (auto adetUnitIt = alcts->begin(); adetUnitIt != alcts->end();
           adetUnitIt++) {
        const CSCDetId& id = (*adetUnitIt).first;
        if (id == (*chamberIdIt)) {
          const auto& range = (*adetUnitIt).second;
          for (auto digiIt = range.first;
               digiIt != range.second; digiIt++) {
            if (digiIt->isValid()) {
              // Check the distance??
              LogTrace("CSCTriggerPrimitivesReader") << "ALCT was found";
              isALCT = true;
              break;
            }
          }
        }
        if (isALCT) break;
      }
      if (isALCT) {
        hEfficALCTEta[station-1]->Fill(fabs(hitEta));
        hEfficALCTEtaCsc[csctype]->Fill(fabs(hitEta));
      }
      else {
        LogTrace("CSCTriggerPrimitivesReader") << "ALCT was not found";
      }

      bool isCLCT = false;
      for (auto cdetUnitIt = clcts->begin(); cdetUnitIt != clcts->end();
           cdetUnitIt++) {
        const CSCDetId& id = (*cdetUnitIt).first;
        if (id == (*chamberIdIt)) {
          const auto& range = (*cdetUnitIt).second;
          for (auto digiIt = range.first;
               digiIt != range.second; digiIt++) {
            if (digiIt->isValid()) {
              // Check the distance??
              LogTrace("CSCTriggerPrimitivesReader") << "CLCT was found";
              isCLCT = true;
              break;
            }
          }
        }
        if (isCLCT) break;
      }
      if (isCLCT) {
        hEfficCLCTEta[station-1]->Fill(fabs(hitEta));
        hEfficCLCTEtaCsc[csctype]->Fill(fabs(hitEta));
      }
      else {
        LogTrace("CSCTriggerPrimitivesReader") << "CLCT was not found";
      }

    }
  }
}

void CSCTriggerPrimitivesReader::drawALCTHistos() {
  TCanvas *c1 = new TCanvas("c1", "", 0, 0, 500, 640);
  string fname = resultsFileNamesPrefix_+"alcts.ps";
  TPostScript *ps = new TPostScript(fname.c_str(), 111);

  TPad *pad[MAXPAGES];
  for (int i_page = 0; i_page < MAXPAGES; i_page++) {
    pad[i_page] = new TPad("", "", .05, .05, .93, .93);
  }

  int page = 1;
  TText t;
  t.SetTextFont(32);
  t.SetTextSize(0.025);
  char pagenum[7], titl[50];
  TPaveLabel *title;

  int max_idh = plotME42 ? CSC_TYPES : CSC_TYPES-1;

  ps->NewPage();
  c1->Clear();  c1->cd(0);
  title = new TPaveLabel(0.1, 0.94, 0.9, 0.98, "Number of ALCTs");
  title->SetFillColor(10);  title->Draw();
  sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
  gStyle->SetOptStat(111110);
  pad[page]->Draw();
  pad[page]->Divide(1,3);
  pad[page]->cd(1);  hAlctPerEvent->Draw();
  pad[page]->cd(2);  hAlctPerChamber->Draw();
  for (int i = 0; i < CSC_TYPES; i++) {
    hAlctPerCSC->GetXaxis()->SetBinLabel(i+1, csc_type[i].c_str());
  }
  // Should be multiplied by 40/nevents to convert to MHz
  pad[page]->cd(3);  hAlctPerCSC->Draw();
  page++;  c1->Update();

  for (int endc = 0; endc < MAX_ENDCAPS; endc++) {
    ps->NewPage();
    c1->Clear();  c1->cd(0);
    sprintf(titl, "ALCTs per chamber, endcap %d", endc+1);
    title = new TPaveLabel(0.1, 0.94, 0.9, 0.98, titl);
    title->SetFillColor(10);  title->Draw();
    sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
    gStyle->SetOptStat(10);
    pad[page]->Draw();
    pad[page]->Divide(2,5);
    for (int idh = 0; idh < max_idh; idh++) {
      if (!plotME1A && idh == 3) continue;
      hAlctCsc[endc][idh]->SetMinimum(0.0);
      pad[page]->cd(idh+1);  hAlctCsc[endc][idh]->Draw();
    }
    page++;  c1->Update();
  }

  ps->NewPage();
  c1->Clear();  c1->cd(0);
  title = new TPaveLabel(0.1, 0.94, 0.9, 0.98, "ALCT quantities");
  title->SetFillColor(10);  title->Draw();
  sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
  gStyle->SetOptStat(111110);
  pad[page]->Draw();
  pad[page]->Divide(2,3);
  pad[page]->cd(1);  hAlctValid->Draw();
  pad[page]->cd(2);  hAlctQuality->Draw();
  pad[page]->cd(3);  hAlctAccel->Draw();
  pad[page]->cd(4);  hAlctCollis->Draw();
  pad[page]->cd(5);  hAlctKeyGroup->Draw();
  pad[page]->cd(6);  hAlctBXN->Draw();
  page++;  c1->Update();

  ps->Close();
  delete c1;
}

void CSCTriggerPrimitivesReader::drawCLCTHistos() {
  TCanvas *c1 = new TCanvas("c1", "", 0, 0, 500, 640);
  string fname = resultsFileNamesPrefix_+"clcts.ps";
  TPostScript *ps = new TPostScript(fname.c_str(), 111);

  TPad *pad[MAXPAGES];
  for (int i_page = 0; i_page < MAXPAGES; i_page++) {
    pad[i_page] = new TPad("", "", .05, .05, .93, .93);
  }

  int page = 1;
  TText t;
  t.SetTextFont(32);
  t.SetTextSize(0.025);
  char pagenum[7], titl[50];
  TPaveLabel *title;

  ps->NewPage();
  c1->Clear();  c1->cd(0);
  title = new TPaveLabel(0.1, 0.94, 0.9, 0.98, "Number of CLCTs");
  title->SetFillColor(10);  title->Draw();
  sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
  gStyle->SetOptStat(111110);
  pad[page]->Draw();
  pad[page]->Divide(1,3);
  pad[page]->cd(1);  hClctPerEvent->Draw();

  int max_idh = plotME42 ? CSC_TYPES : CSC_TYPES-1;

  edm::LogInfo("CSCTriggerPrimitivesReader") << "\n";
  int nbins = hClctPerChamber->GetNbinsX();
  for (int ibin = 1; ibin <= nbins; ibin++) {
    double f_bin = hClctPerChamber->GetBinContent(ibin);
    edm::LogInfo("CSCTriggerPrimitivesReader")
      << "  # CLCTs/chamber: " << ibin-1 << "; events: " << f_bin << endl;
  }

  pad[page]->cd(2);  hClctPerChamber->Draw();
  for (int i = 0; i < CSC_TYPES; i++) {
    hClctPerCSC->GetXaxis()->SetBinLabel(i+1, csc_type[i].c_str());
  }
  // Should be multiplied by 40/nevents to convert to MHz
  pad[page]->cd(3);  hClctPerCSC->Draw();
  page++;  c1->Update();

  for (int endc = 0; endc < MAX_ENDCAPS; endc++) {
    ps->NewPage();
    c1->Clear();  c1->cd(0);
    sprintf(titl, "CLCTs per chamber, endcap %d", endc+1);
    title = new TPaveLabel(0.1, 0.94, 0.9, 0.98, titl);
    title->SetFillColor(10);  title->Draw();
    sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
    gStyle->SetOptStat(10);
    pad[page]->Draw();
    pad[page]->Divide(2,5);
    for (int idh = 0; idh < max_idh; idh++) {
      if (!plotME1A && idh == 3) continue;
      hClctCsc[endc][idh]->SetMinimum(0.0);
      pad[page]->cd(idh+1);  hClctCsc[endc][idh]->Draw();
    }
    page++;  c1->Update();
  }

  ps->NewPage();
  c1->Clear();  c1->cd(0);
  title = new TPaveLabel(0.1, 0.94, 0.9, 0.98, "CLCT quantities");
  title->SetFillColor(10);  title->Draw();
  sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
  gStyle->SetOptStat(111110);
  pad[page]->Draw();
  pad[page]->Divide(2,4);
  pad[page]->cd(1);  hClctValid->Draw();
  pad[page]->cd(2);  hClctQuality->Draw();
  pad[page]->cd(3);  hClctSign->Draw();
  hClctPattern[1]->SetTitle("CLCT pattern #");
  pad[page]->cd(4);  hClctPattern[1]->Draw();
  pad[page]->cd(5);  hClctCFEB->Draw();
  pad[page]->cd(7);  hClctKeyStrip[1]->Draw();
  pad[page]->cd(8);  hClctBXN->Draw();
  page++;  c1->Update();

  ps->NewPage();
  c1->Clear();  c1->cd(0);
  title = new TPaveLabel(0.1, 0.94, 0.9, 0.98,
                         "CLCT bend for various chamber types, halfstrips");
  title->SetFillColor(10);  title->Draw();
  sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
  gStyle->SetOptStat(110);
  pad[page]->Draw();
  pad[page]->Divide(2,5);
  for (int idh = 0; idh < max_idh; idh++) {
    if (!plotME1A && idh == 3) continue;
    pad[page]->cd(idh+1);
    hClctBendCsc[idh][1]->GetXaxis()->SetTitle("Pattern bend");
    hClctBendCsc[idh][1]->GetYaxis()->SetTitle("Number of LCTs");
    hClctBendCsc[idh][1]->Draw();
  }
  page++;  c1->Update();

  ps->Close();
  delete c1;
}

void CSCTriggerPrimitivesReader::drawLCTTMBHistos() {
  TCanvas *c1 = new TCanvas("c1", "", 0, 0, 500, 640);
  string fname = resultsFileNamesPrefix_+"lcts_tmb.ps";
  TPostScript *ps = new TPostScript(fname.c_str(), 111);

  TPad *pad[MAXPAGES];
  for (int i_page = 0; i_page < MAXPAGES; i_page++) {
    pad[i_page] = new TPad("", "", .05, .05, .93, .93);
  }

  int page = 1;
  TText t;
  t.SetTextFont(32);
  t.SetTextSize(0.025);
  char pagenum[7], titl[50];
  TPaveLabel *title;

  int max_idh = plotME42 ? CSC_TYPES : CSC_TYPES-1;

  ps->NewPage();
  c1->Clear();  c1->cd(0);
  title = new TPaveLabel(0.1, 0.94, 0.9, 0.98, "Number of LCTs");
  title->SetFillColor(10);  title->Draw();
  sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
  gStyle->SetOptStat(111110);
  pad[page]->Draw();
  pad[page]->Divide(2,2);
  pad[page]->cd(1);  hLctTMBPerEvent->Draw();
  pad[page]->cd(2);  hLctTMBPerChamber->Draw();
  c1->Update();
  for (int i = 0; i < CSC_TYPES; i++) {
    hLctTMBPerCSC->GetXaxis()->SetBinLabel(i+1, csc_type[i].c_str());
    hCorrLctTMBPerCSC->GetXaxis()->SetBinLabel(i+1, csc_type[i].c_str());
  }
  // Should be multiplied by 40/nevents to convert to MHz
  pad[page]->cd(3);  hLctTMBPerCSC->Draw();
  pad[page]->cd(4);  hCorrLctTMBPerCSC->Draw();
  gStyle->SetOptStat(1110);
  page++;  c1->Update();

  ps->NewPage();
  c1->Clear();  c1->cd(0);
  title = new TPaveLabel(0.1, 0.94, 0.9, 0.98, "LCT geometry");
  title->SetFillColor(10);  title->Draw();
  sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
  gStyle->SetOptStat(110110);
  pad[page]->Draw();
  pad[page]->Divide(2,4);
  pad[page]->cd(1);  hLctTMBEndcap->Draw();
  pad[page]->cd(2);  hLctTMBStation->Draw();
  pad[page]->cd(3);  hLctTMBSector->Draw();
  pad[page]->cd(4);  hLctTMBRing->Draw();
  for (int istat = 0; istat < MAX_STATIONS; istat++) {
    pad[page]->cd(istat+5);  hLctTMBChamber[istat]->Draw();
  }
  page++;  c1->Update();

  for (int endc = 0; endc < MAX_ENDCAPS; endc++) {
    ps->NewPage();
    c1->Clear();  c1->cd(0);
    sprintf(titl, "LCTs per chamber, endcap %d", endc+1);
    title = new TPaveLabel(0.1, 0.94, 0.9, 0.98, titl);
    title->SetFillColor(10);  title->Draw();
    sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
    gStyle->SetOptStat(10);
    pad[page]->Draw();
    pad[page]->Divide(2,5);
    for (int idh = 0; idh < max_idh; idh++) {
      if (!plotME1A && idh == 3) continue;
      hLctTMBCsc[endc][idh]->SetMinimum(0.0);
      pad[page]->cd(idh+1);  hLctTMBCsc[endc][idh]->Draw();
    }
    page++;  c1->Update();
  }

  ps->NewPage();
  c1->Clear();  c1->cd(0);
  title = new TPaveLabel(0.1, 0.94, 0.9, 0.98, "LCT quantities");
  title->SetFillColor(10);  title->Draw();
  sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
  gStyle->SetOptStat(110110);
  pad[page]->Draw();
  pad[page]->Divide(2,4);
  pad[page]->cd(1);  hLctTMBValid->Draw();
  pad[page]->cd(2);  hLctTMBQuality->Draw();
  pad[page]->cd(3);  hLctTMBKeyGroup->Draw();
  pad[page]->cd(4);  hLctTMBKeyStrip->Draw();
  pad[page]->cd(5);  hLctTMBStripType->Draw();
  pad[page]->cd(6);  hLctTMBPattern->Draw();
  pad[page]->cd(7);  hLctTMBBend->Draw();
  pad[page]->cd(8);  hLctTMBBXN->Draw();
  page++;  c1->Update();

  ps->Close();
  delete c1;
}

void CSCTriggerPrimitivesReader::drawLCTMPCHistos() {
  TCanvas *c1 = new TCanvas("c1", "", 0, 0, 500, 640);
  string fname = resultsFileNamesPrefix_+"lcts_mpc.ps";
  TPostScript *ps = new TPostScript(fname.c_str(), 111);

  TPad *pad[MAXPAGES];
  for (int i_page = 0; i_page < MAXPAGES; i_page++) {
    pad[i_page] = new TPad("", "", .05, .05, .93, .93);
  }

  int page = 1;
  TText t;
  t.SetTextFont(32);
  t.SetTextSize(0.025);
  char pagenum[7];
  TPaveLabel *title;

  ps->NewPage();
  c1->Clear();  c1->cd(0);
  title = new TPaveLabel(0.1, 0.94, 0.9, 0.98, "Number of LCTs");
  title->SetFillColor(10);  title->Draw();
  sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
  gStyle->SetOptStat(111110);
  pad[page]->Draw();
  pad[page]->Divide(1,3);
  pad[page]->cd(1);  hLctMPCPerEvent->Draw();
  for (int i = 0; i < CSC_TYPES; i++) {
    hLctMPCPerCSC->GetXaxis()->SetBinLabel(i+1, csc_type[i].c_str());
    hCorrLctMPCPerCSC->GetXaxis()->SetBinLabel(i+1, csc_type[i].c_str());
  }
  // Should be multiplied by 40/nevents to convert to MHz
  pad[page]->cd(2);  hLctMPCPerCSC->Draw();
  pad[page]->cd(3);  hCorrLctMPCPerCSC->Draw();
  page++;  c1->Update();

  ps->NewPage();
  c1->Clear();  c1->cd(0);
  title = new TPaveLabel(0.1, 0.94, 0.9, 0.98, "LCT geometry");
  title->SetFillColor(10);  title->Draw();
  sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
  gStyle->SetOptStat(110110);
  pad[page]->Draw();
  pad[page]->Divide(2,4);
  pad[page]->cd(1);  hLctMPCEndcap->Draw();
  pad[page]->cd(2);  hLctMPCStation->Draw();
  pad[page]->cd(3);  hLctMPCSector->Draw();
  pad[page]->cd(4);  hLctMPCRing->Draw();
  for (int istat = 0; istat < MAX_STATIONS; istat++) {
    pad[page]->cd(istat+5);  hLctMPCChamber[istat]->Draw();
  }
  page++;  c1->Update();

  ps->NewPage();
  c1->Clear();  c1->cd(0);
  title = new TPaveLabel(0.1, 0.94, 0.9, 0.98, "LCT quantities");
  title->SetFillColor(10);  title->Draw();
  sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
  pad[page]->Draw();
  pad[page]->Divide(2,4);
  pad[page]->cd(1);  hLctMPCValid->Draw();
  pad[page]->cd(2);  hLctMPCQuality->Draw();
  pad[page]->cd(3);  hLctMPCKeyGroup->Draw();
  pad[page]->cd(4);  hLctMPCKeyStrip->Draw();
  pad[page]->cd(5);  hLctMPCStripType->Draw();
  pad[page]->cd(6);  hLctMPCPattern->Draw();
  pad[page]->cd(7);  hLctMPCBend->Draw();
  pad[page]->cd(8);  hLctMPCBXN->Draw();
  page++;  c1->Update();

  ps->Close();
  delete c1;
}

void CSCTriggerPrimitivesReader::drawCompHistos() {
  TCanvas *c1 = new TCanvas("c1", "", 0, 0, 500, 640);
  string fname = resultsFileNamesPrefix_+"lcts_comp.ps";
  TPostScript *ps = new TPostScript(fname.c_str(), 111);

  TPad *pad[MAXPAGES];
  for (int i_page = 0; i_page < MAXPAGES; i_page++) {
    pad[i_page] = new TPad("", "", .05, .05, .93, .93);
  }

  int page = 1;
  TText t;
  t.SetTextFont(32);
  t.SetTextSize(0.025);
  char pagenum[15];
  TPaveLabel *title;
  Int_t nbins;

  TText teff;
  teff.SetTextFont(32);
  teff.SetTextSize(0.08);
  char eff[25], titl[60];

  int max_idh = plotME42 ? CSC_TYPES : CSC_TYPES-1;

  for (int endc = 0; endc < MAX_ENDCAPS; endc++) { // endcaps
    ps->NewPage();
    c1->Clear();  c1->cd(0);
    sprintf(titl, "ALCT firmware-emulator: match in number found, endcap %d",
            endc+1);
    title = new TPaveLabel(0.1, 0.94, 0.9, 0.98, titl);
    title->SetFillColor(10);  title->Draw();
    sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
    //gStyle->SetOptStat(110010);
    gStyle->SetOptStat(0);
    pad[page]->Draw();
    pad[page]->Divide(2,5);
    TH1F *hAlctFoundEffVsCsc[MAX_ENDCAPS][CSC_TYPES];
    for (int idh = 0; idh < max_idh; idh++) {
      if (!plotME1A && idh == 3) continue;
      hAlctFoundEffVsCsc[endc][idh] =
        (TH1F*)hAlctCompFoundCsc[endc][idh]->Clone();
      hAlctFoundEffVsCsc[endc][idh]->Divide(hAlctCompSameNCsc[endc][idh],
                                            hAlctCompFoundCsc[endc][idh],
                                            1., 1., "B");
      nbins = hAlctCompFoundCsc[endc][idh]->GetNbinsX();
      for (Int_t ibin = 1; ibin <= nbins; ibin++) {
        if (hAlctCompFoundCsc[endc][idh]->GetBinContent(ibin) == 0) {
          hAlctFoundEffVsCsc[endc][idh]->SetBinContent(ibin, -1.);
          hAlctFoundEffVsCsc[endc][idh]->SetBinError(ibin, 0.);
        }
      }
      gPad->Update();  gStyle->SetStatX(0.65);
      hAlctFoundEffVsCsc[endc][idh]->SetMinimum(-0.05);
      hAlctFoundEffVsCsc[endc][idh]->SetMaximum(1.05);
      hAlctFoundEffVsCsc[endc][idh]->GetXaxis()->SetTitleOffset(0.7);
      hAlctFoundEffVsCsc[endc][idh]->GetYaxis()->SetTitleOffset(0.8);
      hAlctFoundEffVsCsc[endc][idh]->GetXaxis()->SetLabelSize(0.06); // default=0.04
      hAlctFoundEffVsCsc[endc][idh]->GetYaxis()->SetLabelSize(0.06); // default=0.04
      hAlctFoundEffVsCsc[endc][idh]->GetXaxis()->SetTitleSize(0.07); // default=0.05
      hAlctFoundEffVsCsc[endc][idh]->GetYaxis()->SetTitleSize(0.07); // default=0.05
      hAlctFoundEffVsCsc[endc][idh]->GetXaxis()->SetTitle("CSC id");
      hAlctFoundEffVsCsc[endc][idh]->GetYaxis()->SetTitle("% of same number found");
      pad[page]->cd(idh+1);  hAlctFoundEffVsCsc[endc][idh]->Draw("e");
      double numer = hAlctCompSameNCsc[endc][idh]->Integral();
      double denom = hAlctCompFoundCsc[endc][idh]->Integral();
      double ratio = 0.0, error = 0.0;
      if (denom > 0.) {
        ratio = numer/denom;
        error = sqrt(ratio*(1.-ratio)/denom);
      }
      sprintf(eff, "eff = (%4.1f +/- %4.1f)%%", ratio*100., error*100.);
      teff.DrawTextNDC(0.3, 0.5, eff);
    }
    page++;  c1->Update();
  }

  for (int endc = 0; endc < MAX_ENDCAPS; endc++) { // endcaps
    ps->NewPage();
    c1->Clear();  c1->cd(0);
    sprintf(titl, "ALCT firmware-emulator: exact match, endcap %d", endc+1);
    title = new TPaveLabel(0.1, 0.94, 0.9, 0.98, titl);
    title->SetFillColor(10);  title->Draw();
    sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
    //gStyle->SetOptStat(110010);
    gStyle->SetOptStat(0);
    pad[page]->Draw();
    pad[page]->Divide(2,5);
    TH1F *hAlctMatchEffVsCsc[MAX_ENDCAPS][CSC_TYPES];
    for (int idh = 0; idh < max_idh; idh++) {
      if (!plotME1A && idh == 3) continue;
      hAlctMatchEffVsCsc[endc][idh] =
        (TH1F*)hAlctCompTotalCsc[endc][idh]->Clone();
      hAlctMatchEffVsCsc[endc][idh]->Divide(hAlctCompMatchCsc[endc][idh],
                                            hAlctCompTotalCsc[endc][idh],
                                            1., 1., "B");
      nbins = hAlctCompTotalCsc[endc][idh]->GetNbinsX();
      for (Int_t ibin = 1; ibin <= nbins; ibin++) {
        if (hAlctCompTotalCsc[endc][idh]->GetBinContent(ibin) == 0) {
          hAlctMatchEffVsCsc[endc][idh]->SetBinContent(ibin, -1.);
          hAlctMatchEffVsCsc[endc][idh]->SetBinError(ibin, 0.);
        }
      }
      gPad->Update();  gStyle->SetStatX(0.65);
      hAlctMatchEffVsCsc[endc][idh]->SetMinimum(-0.05);
      hAlctMatchEffVsCsc[endc][idh]->SetMaximum(1.05);
      hAlctMatchEffVsCsc[endc][idh]->GetXaxis()->SetTitleOffset(0.7);
      hAlctMatchEffVsCsc[endc][idh]->GetYaxis()->SetTitleOffset(0.8);
      hAlctMatchEffVsCsc[endc][idh]->GetXaxis()->SetLabelSize(0.06);
      hAlctMatchEffVsCsc[endc][idh]->GetYaxis()->SetLabelSize(0.06);
      hAlctMatchEffVsCsc[endc][idh]->GetXaxis()->SetTitleSize(0.07);
      hAlctMatchEffVsCsc[endc][idh]->GetYaxis()->SetTitleSize(0.07);
      hAlctMatchEffVsCsc[endc][idh]->GetXaxis()->SetTitle("CSC id");
      hAlctMatchEffVsCsc[endc][idh]->GetYaxis()->SetTitle("% of exact match");
      pad[page]->cd(idh+1);  hAlctMatchEffVsCsc[endc][idh]->Draw("e");
      double numer = hAlctCompMatchCsc[endc][idh]->Integral();
      double denom = hAlctCompTotalCsc[endc][idh]->Integral();
      double ratio = 0.0, error = 0.0;
      if (denom > 0.) {
        ratio = numer/denom;
        error = sqrt(ratio*(1.-ratio)/denom);
      }
      sprintf(eff, "eff = (%4.1f +/- %4.1f)%%", ratio*100., error*100.);
      teff.DrawTextNDC(0.3, 0.5, eff);
    }
    page++;  c1->Update();
  }

  for (int endc = 0; endc < MAX_ENDCAPS; endc++) {
    ps->NewPage();
    c1->Clear();  c1->cd(0);
    sprintf(titl, "CLCT firmware-emulator: match in number found, endcap %d",
            endc+1);
    title = new TPaveLabel(0.1, 0.94, 0.9, 0.98, titl);
    title->SetFillColor(10);  title->Draw();
    sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
    //gStyle->SetOptStat(110010);
    gStyle->SetOptStat(0);
    pad[page]->Draw();
    pad[page]->Divide(2,5);
    TH1F *hClctFoundEffVsCsc[MAX_ENDCAPS][CSC_TYPES];
    for (int idh = 0; idh < max_idh; idh++) {
      if (!plotME1A && idh == 3) continue;
      hClctFoundEffVsCsc[endc][idh] =
        (TH1F*)hClctCompFoundCsc[endc][idh]->Clone();
      hClctFoundEffVsCsc[endc][idh]->Divide(hClctCompSameNCsc[endc][idh],
                                            hClctCompFoundCsc[endc][idh],
                                            1., 1., "B");
      nbins = hClctCompFoundCsc[endc][idh]->GetNbinsX();
      for (Int_t ibin = 1; ibin <= nbins; ibin++) {
        if (hClctCompFoundCsc[endc][idh]->GetBinContent(ibin) == 0) {
          hClctFoundEffVsCsc[endc][idh]->SetBinContent(ibin, -1.);
          hClctFoundEffVsCsc[endc][idh]->SetBinError(ibin, 0.);
        }
      }
      gPad->Update();  gStyle->SetStatX(0.65);
      hClctFoundEffVsCsc[endc][idh]->SetMinimum(-0.05);
      hClctFoundEffVsCsc[endc][idh]->SetMaximum(1.05);
      hClctFoundEffVsCsc[endc][idh]->GetXaxis()->SetTitleOffset(0.7);
      hClctFoundEffVsCsc[endc][idh]->GetYaxis()->SetTitleOffset(0.8);
      hClctFoundEffVsCsc[endc][idh]->GetXaxis()->SetLabelSize(0.06);
      hClctFoundEffVsCsc[endc][idh]->GetYaxis()->SetLabelSize(0.06);
      hClctFoundEffVsCsc[endc][idh]->GetXaxis()->SetTitleSize(0.07);
      hClctFoundEffVsCsc[endc][idh]->GetYaxis()->SetTitleSize(0.07);
      hClctFoundEffVsCsc[endc][idh]->GetXaxis()->SetTitle("CSC id");
      hClctFoundEffVsCsc[endc][idh]->GetYaxis()->SetTitle("% of same number found");
      pad[page]->cd(idh+1);  hClctFoundEffVsCsc[endc][idh]->Draw("e");
      double numer = hClctCompSameNCsc[endc][idh]->Integral();
      double denom = hClctCompFoundCsc[endc][idh]->Integral();
      double ratio = 0.0, error = 0.0;
      if (denom > 0.) {
        ratio = numer/denom;
        error = sqrt(ratio*(1.-ratio)/denom);
      }
      sprintf(eff, "eff = (%4.1f +/- %4.1f)%%", ratio*100., error*100.);
      teff.DrawTextNDC(0.3, 0.5, eff);
    }
    page++;  c1->Update();
  }

  for (int endc = 0; endc < MAX_ENDCAPS; endc++) {
    ps->NewPage();
    c1->Clear();  c1->cd(0);
    sprintf(titl, "CLCT firmware-emulator: exact match, endcap %d", endc+1);
    title = new TPaveLabel(0.1, 0.94, 0.9, 0.98, titl);
    title->SetFillColor(10);  title->Draw();
    sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
    //gStyle->SetOptStat(110010);
    gStyle->SetOptStat(0);
    pad[page]->Draw();
    pad[page]->Divide(2,5);
    TH1F *hClctMatchEffVsCsc[MAX_ENDCAPS][CSC_TYPES];
    for (int idh = 0; idh < max_idh; idh++) {
      if (!plotME1A && idh == 3) continue;
      hClctMatchEffVsCsc[endc][idh] =
        (TH1F*)hClctCompTotalCsc[endc][idh]->Clone();
      hClctMatchEffVsCsc[endc][idh]->Divide(hClctCompMatchCsc[endc][idh],
                                            hClctCompTotalCsc[endc][idh],
                                            1., 1., "B");
      nbins = hClctCompTotalCsc[endc][idh]->GetNbinsX();
      for (Int_t ibin = 1; ibin <= nbins; ibin++) {
        if (hClctCompTotalCsc[endc][idh]->GetBinContent(ibin) == 0) {
          hClctMatchEffVsCsc[endc][idh]->SetBinContent(ibin, -1.);
          hClctMatchEffVsCsc[endc][idh]->SetBinError(ibin, 0.);
        }
      }
      gPad->Update();  gStyle->SetStatX(0.65);
      hClctMatchEffVsCsc[endc][idh]->SetMinimum(-0.05);
      hClctMatchEffVsCsc[endc][idh]->SetMaximum(1.05);
      hClctMatchEffVsCsc[endc][idh]->GetXaxis()->SetTitleOffset(0.7);
      hClctMatchEffVsCsc[endc][idh]->GetYaxis()->SetTitleOffset(0.8);
      hClctMatchEffVsCsc[endc][idh]->GetXaxis()->SetLabelSize(0.06);
      hClctMatchEffVsCsc[endc][idh]->GetYaxis()->SetLabelSize(0.06);
      hClctMatchEffVsCsc[endc][idh]->GetXaxis()->SetTitleSize(0.07);
      hClctMatchEffVsCsc[endc][idh]->GetYaxis()->SetTitleSize(0.07);
      hClctMatchEffVsCsc[endc][idh]->GetXaxis()->SetTitle("CSC id");
      hClctMatchEffVsCsc[endc][idh]->GetYaxis()->SetTitle("% of exact match");
      pad[page]->cd(idh+1);  hClctMatchEffVsCsc[endc][idh]->Draw("e");
      double numer = hClctCompMatchCsc[endc][idh]->Integral();
      double denom = hClctCompTotalCsc[endc][idh]->Integral();
      double ratio = 0.0, error = 0.0;
      if (denom > 0.) {
        ratio = numer/denom;
        error = sqrt(ratio*(1.-ratio)/denom);
      }
      sprintf(eff, "eff = (%4.1f +/- %4.1f)%%", ratio*100., error*100.);
      teff.DrawTextNDC(0.3, 0.5, eff);
    }
    page++;  c1->Update();
  }

  for (int endc = 0; endc < MAX_ENDCAPS; endc++) {
    ps->NewPage();
    c1->Clear();  c1->cd(0);
    sprintf(titl, "LCT firmware-emulator: match in number found, endcap %d",
            endc+1);
    title = new TPaveLabel(0.1, 0.94, 0.9, 0.98, titl);
    title->SetFillColor(10);  title->Draw();
    sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
    //gStyle->SetOptStat(110010);
    gStyle->SetOptStat(0);
    pad[page]->Draw();
    pad[page]->Divide(2,5);
    TH1F *hLctFoundEffVsCsc[MAX_ENDCAPS][CSC_TYPES];
    for (int idh = 0; idh < max_idh; idh++) {
      if (!plotME1A && idh == 3) continue;
      hLctFoundEffVsCsc[endc][idh] =
        (TH1F*)hLctCompFoundCsc[endc][idh]->Clone();
      hLctFoundEffVsCsc[endc][idh]->Divide(hLctCompSameNCsc[endc][idh],
                                           hLctCompFoundCsc[endc][idh],
                                           1., 1., "B");
      nbins = hLctCompFoundCsc[endc][idh]->GetNbinsX();
      for (Int_t ibin = 1; ibin <= nbins; ibin++) {
        if (hLctCompFoundCsc[endc][idh]->GetBinContent(ibin) == 0) {
          hLctFoundEffVsCsc[endc][idh]->SetBinContent(ibin, -1.);
          hLctFoundEffVsCsc[endc][idh]->SetBinError(ibin, 0.);
        }
      }
      gPad->Update();  gStyle->SetStatX(0.65);
      hLctFoundEffVsCsc[endc][idh]->SetMinimum(-0.05);
      hLctFoundEffVsCsc[endc][idh]->SetMaximum(1.05);
      hLctFoundEffVsCsc[endc][idh]->GetXaxis()->SetTitleOffset(0.7);
      hLctFoundEffVsCsc[endc][idh]->GetYaxis()->SetTitleOffset(0.8);
      hLctFoundEffVsCsc[endc][idh]->GetXaxis()->SetLabelSize(0.06);
      hLctFoundEffVsCsc[endc][idh]->GetYaxis()->SetLabelSize(0.06);
      hLctFoundEffVsCsc[endc][idh]->GetXaxis()->SetTitleSize(0.07);
      hLctFoundEffVsCsc[endc][idh]->GetYaxis()->SetTitleSize(0.07);
      hLctFoundEffVsCsc[endc][idh]->GetXaxis()->SetTitle("CSC id");
      hLctFoundEffVsCsc[endc][idh]->GetYaxis()->SetTitle("% of same number found");
      pad[page]->cd(idh+1);  hLctFoundEffVsCsc[endc][idh]->Draw("e");
      double numer = hLctCompSameNCsc[endc][idh]->Integral();
      double denom = hLctCompFoundCsc[endc][idh]->Integral();
      double ratio = 0.0, error = 0.0;
      if (denom > 0.) {
        ratio = numer/denom;
        error = sqrt(ratio*(1.-ratio)/denom);
      }
      sprintf(eff, "eff = (%4.1f +/- %4.1f)%%", ratio*100., error*100.);
      teff.DrawTextNDC(0.3, 0.5, eff);
    }
    page++;  c1->Update();
  }

  for (int endc = 0; endc < MAX_ENDCAPS; endc++) {
    ps->NewPage();
    c1->Clear();  c1->cd(0);
    sprintf(titl, "LCT firmware-emulator: exact match, endcap %d", endc+1);
    title = new TPaveLabel(0.1, 0.94, 0.9, 0.98, titl);
    title->SetFillColor(10);  title->Draw();
    sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
    //gStyle->SetOptStat(110010);
    gStyle->SetOptStat(0);
    pad[page]->Draw();
    pad[page]->Divide(2,5);
    TH1F *hLctMatchEffVsCsc[MAX_ENDCAPS][CSC_TYPES];
    for (int idh = 0; idh < max_idh; idh++) {
      if (!plotME1A && idh == 3) continue;
      hLctMatchEffVsCsc[endc][idh] =
        (TH1F*)hLctCompTotalCsc[endc][idh]->Clone();
      hLctMatchEffVsCsc[endc][idh]->Divide(hLctCompMatchCsc[endc][idh],
                                           hLctCompTotalCsc[endc][idh],
                                           1., 1., "B");
      nbins = hLctCompTotalCsc[endc][idh]->GetNbinsX();
      for (Int_t ibin = 1; ibin <= nbins; ibin++) {
        if (hLctCompTotalCsc[endc][idh]->GetBinContent(ibin) == 0) {
          hLctMatchEffVsCsc[endc][idh]->SetBinContent(ibin, -1.);
          hLctMatchEffVsCsc[endc][idh]->SetBinError(ibin, 0.);
        }
      }
      gPad->Update();  gStyle->SetStatX(0.65);
      hLctMatchEffVsCsc[endc][idh]->SetMinimum(-0.05);
      hLctMatchEffVsCsc[endc][idh]->SetMaximum(1.05);
      hLctMatchEffVsCsc[endc][idh]->GetXaxis()->SetTitleOffset(0.7);
      hLctMatchEffVsCsc[endc][idh]->GetYaxis()->SetTitleOffset(0.8);
      hLctMatchEffVsCsc[endc][idh]->GetXaxis()->SetLabelSize(0.06);
      hLctMatchEffVsCsc[endc][idh]->GetYaxis()->SetLabelSize(0.06);
      hLctMatchEffVsCsc[endc][idh]->GetXaxis()->SetTitleSize(0.07);
      hLctMatchEffVsCsc[endc][idh]->GetYaxis()->SetTitleSize(0.07);
      hLctMatchEffVsCsc[endc][idh]->GetXaxis()->SetTitle("CSC id");
      hLctMatchEffVsCsc[endc][idh]->GetYaxis()->SetTitle("% of exact match");
      pad[page]->cd(idh+1);  hLctMatchEffVsCsc[endc][idh]->Draw("e");
      double numer = hLctCompMatchCsc[endc][idh]->Integral();
      double denom = hLctCompTotalCsc[endc][idh]->Integral();
      double ratio = 0.0, error = 0.0;
      if (denom > 0.) {
        ratio = numer/denom;
        error = sqrt(ratio*(1.-ratio)/denom);
      }
      sprintf(eff, "eff = (%4.1f +/- %4.1f)%%", ratio*100., error*100.);
      teff.DrawTextNDC(0.3, 0.5, eff);
    }
    page++;  c1->Update();
  }

  ps->Close();
  delete c1;
}

void CSCTriggerPrimitivesReader::drawResolHistos() {
  TCanvas *c1 = new TCanvas("c1", "", 0, 0, 500, 640);
  string fname = resultsFileNamesPrefix_+"lcts_resol.ps";
  TPostScript *ps = new TPostScript(fname.c_str(), 111);

  TPad *pad[MAXPAGES];
  for (int i_page = 0; i_page < MAXPAGES; i_page++) {
    pad[i_page] = new TPad("", "", .05, .05, .93, .93);
  }

  int page = 1;
  TText t;
  t.SetTextFont(32);
  t.SetTextSize(0.025);
  char pagenum[7];
  TPaveLabel *title;

  int max_idh = plotME42 ? CSC_TYPES : CSC_TYPES-1;

  ps->NewPage();
  c1->Clear();  c1->cd(0);
  title = new TPaveLabel(0.1, 0.94, 0.9, 0.98, "ALCT resolution");
  title->SetFillColor(10);  title->Draw();
  sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
  gStyle->SetOptStat(111110);
  pad[page]->Draw();
  pad[page]->Divide(2,2);
  gStyle->SetStatX(1.00);  gStyle->SetStatY(0.65);
  pad[page]->cd(1);  hEtaRecVsSim->SetMarkerSize(0.2);  hEtaRecVsSim->Draw();
  gPad->Update();  gStyle->SetStatX(1.00);  gStyle->SetStatY(0.995);
  hResolDeltaWG->GetXaxis()->SetTitle("WG_{rec} - WG_{sim}");
  hResolDeltaWG->GetXaxis()->SetTitleOffset(1.2);
  hResolDeltaWG->GetYaxis()->SetTitle("Entries");
  hResolDeltaWG->GetYaxis()->SetTitleOffset(1.9);
  hResolDeltaWG->GetXaxis()->SetLabelSize(0.03);
  hResolDeltaWG->GetYaxis()->SetLabelSize(0.03);
  pad[page]->cd(3);  hResolDeltaWG->Draw();
  hResolDeltaEta->GetXaxis()->SetNdivisions(505); // twice fewer divisions
  hResolDeltaEta->GetXaxis()->SetLabelSize(0.04);
  hResolDeltaEta->GetYaxis()->SetLabelSize(0.04);
  pad[page]->cd(4);  hResolDeltaEta->Draw();  hResolDeltaEta->Fit("gaus","Q");
  page++;  c1->Update();

  ps->NewPage();
  c1->Clear();  c1->cd(0);
  title = new TPaveLabel(0.1, 0.94, 0.9, 0.98, "#eta_rec-#eta_sim");
  title->SetFillColor(10);  title->Draw();
  sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
  pad[page]->Draw();
  pad[page]->Divide(2,5);
  for (int idh = 0; idh < max_idh; idh++) {
    if (!plotME1A && idh == 3) continue;
    hEtaDiffCsc[idh][0]->GetXaxis()->SetLabelSize(0.07); // default=0.04
    hEtaDiffCsc[idh][0]->GetYaxis()->SetLabelSize(0.07);
    pad[page]->cd(idh+1);  hEtaDiffCsc[idh][0]->Draw();
    if (hEtaDiffCsc[idh][0]->GetEntries() > 1)
      hEtaDiffCsc[idh][0]->Fit("gaus","Q");
  }
  page++;  c1->Update();

  ps->NewPage();
  c1->Clear();  c1->cd(0);
  title = new TPaveLabel(0.1, 0.94, 0.9, 0.98, "#eta_rec-#eta_sim, endcap1");
  title->SetFillColor(10);  title->Draw();
  sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
  pad[page]->Draw();
  pad[page]->Divide(2,5);
  for (int idh = 0; idh < max_idh; idh++) {
    if (!plotME1A && idh == 3) continue;
    hEtaDiffCsc[idh][1]->GetXaxis()->SetLabelSize(0.07); // default=0.04
    hEtaDiffCsc[idh][1]->GetYaxis()->SetLabelSize(0.07);
    pad[page]->cd(idh+1);  hEtaDiffCsc[idh][1]->Draw();
    if (hEtaDiffCsc[idh][1]->GetEntries() > 1)
      hEtaDiffCsc[idh][1]->Fit("gaus","Q");
  }
  page++;  c1->Update();

  ps->NewPage();
  c1->Clear();  c1->cd(0);
  title = new TPaveLabel(0.1, 0.94, 0.9, 0.98, "#eta_rec-#eta_sim, endcap2");
  title->SetFillColor(10);  title->Draw();
  sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
  pad[page]->Draw();
  pad[page]->Divide(2,5);
  for (int idh = 0; idh < max_idh; idh++) {
    if (!plotME1A && idh == 3) continue;
    hEtaDiffCsc[idh][2]->GetXaxis()->SetLabelSize(0.07); // default=0.04
    hEtaDiffCsc[idh][2]->GetYaxis()->SetLabelSize(0.07);
    pad[page]->cd(idh+1);  hEtaDiffCsc[idh][2]->Draw();
    if (hEtaDiffCsc[idh][2]->GetEntries() > 1)
      hEtaDiffCsc[idh][2]->Fit("gaus","Q");
  }
  page++;  c1->Update();

  ps->NewPage();
  c1->Clear();  c1->cd(0);
  title = new TPaveLabel(0.1, 0.94, 0.9, 0.98,
                         "#LT#eta_rec-#eta_sim#GT vs #eta_rec");
  title->SetFillColor(10);  title->Draw();
  sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
  //gStyle->SetOptStat(0);
  pad[page]->Draw();
  pad[page]->Divide(2,2);
  TH1F *hMeanEtaDiffVsEta[MAX_STATIONS];
  for (int istation = 0; istation < MAX_STATIONS; istation++) {
    hMeanEtaDiffVsEta[istation] = (TH1F*)hEtaDiffVsEta[istation]->Clone();
    hMeanEtaDiffVsEta[istation]->Divide(hEtaDiffVsEta[istation],
                                        hAlctVsEta[istation], 1., 1.);
    hMeanEtaDiffVsEta[istation]->GetXaxis()->SetTitleOffset(1.2);
    hMeanEtaDiffVsEta[istation]->GetXaxis()->SetTitle("#eta");
    hMeanEtaDiffVsEta[istation]->SetMaximum(0.05);
    pad[page]->cd(istation+1);  hMeanEtaDiffVsEta[istation]->Draw();
  }
  page++;  c1->Update();

  c1->Clear();  c1->cd(0);
  title = new TPaveLabel(0.1, 0.94, 0.9, 0.98,
                         "#eta_rec-#eta_sim vs wiregroup");
  title->SetFillColor(10);  title->Draw();
  sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
  pad[page]->Draw();
  pad[page]->Divide(2,5);
  for (int idh = 0; idh < max_idh; idh++) {
    if (!plotME1A && idh == 3) continue;
    pad[page]->cd(idh+1);  hEtaDiffVsWireCsc[idh]->SetMarkerSize(0.2);
    hEtaDiffVsWireCsc[idh]->GetXaxis()->SetTitle("Wiregroup");
    hEtaDiffVsWireCsc[idh]->GetXaxis()->SetTitleSize(0.07);
    hEtaDiffVsWireCsc[idh]->GetXaxis()->SetTitleOffset(1.2);
    hEtaDiffVsWireCsc[idh]->GetYaxis()->SetTitle("#eta_rec-#eta_sim");
    hEtaDiffVsWireCsc[idh]->GetYaxis()->SetTitleSize(0.07);
    hEtaDiffVsWireCsc[idh]->GetXaxis()->SetLabelSize(0.07); // default=0.04
    hEtaDiffVsWireCsc[idh]->GetYaxis()->SetLabelSize(0.07);
    hEtaDiffVsWireCsc[idh]->Draw();
  }
  page++;  c1->Update();

  ps->NewPage();
  c1->Clear();  c1->cd(0);
  title = new TPaveLabel(0.1, 0.94, 0.9, 0.98, "#phi resolution");
  title->SetFillColor(10);  title->Draw();
  sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
  gStyle->SetOptStat(111110);
  pad[page]->Draw();
  pad[page]->Divide(2,2);
  gStyle->SetStatX(1.00);  gStyle->SetStatY(0.65);
  pad[page]->cd(1);  hPhiRecVsSim->SetMarkerSize(0.2);  hPhiRecVsSim->Draw();
  gPad->Update();  gStyle->SetStatX(1.00);  gStyle->SetStatY(0.995);
  hResolDeltaPhi->GetXaxis()->SetLabelSize(0.04);
  hResolDeltaPhi->GetYaxis()->SetLabelSize(0.04);
  pad[page]->cd(2);  hResolDeltaPhi->Draw();  hResolDeltaPhi->Fit("gaus","Q");
  hResolDeltaHS->GetXaxis()->SetTitle("HS_{rec} - HS_{sim}");
  hResolDeltaHS->GetXaxis()->SetTitleOffset(1.2);
  hResolDeltaHS->GetYaxis()->SetTitle("Entries");
  hResolDeltaHS->GetYaxis()->SetTitleOffset(1.7);
  hResolDeltaHS->GetXaxis()->SetLabelSize(0.03);
  hResolDeltaHS->GetYaxis()->SetLabelSize(0.03);
  pad[page]->cd(3);  hResolDeltaHS->Draw();
  hResolDeltaDS->GetXaxis()->SetTitle("DS_{rec} - DS_{sim}");
  hResolDeltaDS->GetXaxis()->SetTitleOffset(1.2);
  hResolDeltaDS->GetYaxis()->SetTitle("Entries");
  hResolDeltaDS->GetYaxis()->SetTitleOffset(1.6);
  hResolDeltaDS->GetXaxis()->SetLabelSize(0.04);
  hResolDeltaDS->GetYaxis()->SetLabelSize(0.04);
  pad[page]->cd(4);  hResolDeltaDS->Draw();
  page++;  c1->Update();

  ps->NewPage();
  c1->Clear();  c1->cd(0);
  title = new TPaveLabel(0.1, 0.94, 0.9, 0.98, "#phi_rec-#phi_sim (mrad)");
  title->SetFillColor(10);  title->Draw();
  sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
  pad[page]->Draw();
  pad[page]->Divide(2,5);
  for (int idh = 0; idh < max_idh; idh++) {
    if (!plotME1A && idh == 3) continue;
    hPhiDiffCsc[idh][0]->GetXaxis()->SetLabelSize(0.07); // default=0.04
    hPhiDiffCsc[idh][0]->GetYaxis()->SetLabelSize(0.07);
    pad[page]->cd(idh+1);  hPhiDiffCsc[idh][0]->Draw();
    if (hPhiDiffCsc[idh][0]->GetEntries() > 1)
      hPhiDiffCsc[idh][0]->Fit("gaus","Q");
  }
  page++;  c1->Update();

  ps->NewPage();
  c1->Clear();  c1->cd(0);
  title = new TPaveLabel(0.1, 0.94, 0.9, 0.98,
                         "#phi_rec-#phi_sim (mrad), endcap1");
  title->SetFillColor(10);  title->Draw();
  sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
  pad[page]->Draw();
  pad[page]->Divide(2,5);
  for (int idh = 0; idh < max_idh; idh++) {
    if (!plotME1A && idh == 3) continue;
    hPhiDiffCsc[idh][1]->GetXaxis()->SetLabelSize(0.07); // default=0.04
    hPhiDiffCsc[idh][1]->GetYaxis()->SetLabelSize(0.07);
    pad[page]->cd(idh+1);  hPhiDiffCsc[idh][1]->Draw();
    if (hPhiDiffCsc[idh][1]->GetEntries() > 1)
      hPhiDiffCsc[idh][1]->Fit("gaus","Q");
  }
  page++;  c1->Update();

  ps->NewPage();
  c1->Clear();  c1->cd(0);
  title = new TPaveLabel(0.1, 0.94, 0.9, 0.98,
                         "#phi_rec-#phi_sim (mrad), endcap2");
  title->SetFillColor(10);  title->Draw();
  sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
  pad[page]->Draw();
  pad[page]->Divide(2,5);
  for (int idh = 0; idh < max_idh; idh++) {
    if (!plotME1A && idh == 3) continue;
    hPhiDiffCsc[idh][2]->GetXaxis()->SetLabelSize(0.07); // default=0.04
    hPhiDiffCsc[idh][2]->GetYaxis()->SetLabelSize(0.07);
    pad[page]->cd(idh+1);  hPhiDiffCsc[idh][2]->Draw();
    if (hPhiDiffCsc[idh][2]->GetEntries() > 1)
      hPhiDiffCsc[idh][2]->Fit("gaus","Q");
  }
  page++;  c1->Update();

  ps->NewPage();
  c1->Clear();  c1->cd(0);
  title = new TPaveLabel(0.1, 0.94, 0.9, 0.98,
                         "#LT#phi_rec-#phi_sim#GT (mrad) vs #phi_rec");
  title->SetFillColor(10);  title->Draw();
  sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
  gStyle->SetOptStat(0);
  pad[page]->Draw();
  pad[page]->Divide(2,2);
  TH1F *hMeanPhiDiffVsPhi[MAX_STATIONS];
  for (int istation = 0; istation < MAX_STATIONS; istation++) {
    hMeanPhiDiffVsPhi[istation] = (TH1F*)hPhiDiffVsPhi[istation]->Clone();
    hMeanPhiDiffVsPhi[istation]->Divide(hPhiDiffVsPhi[istation],
                                        hClctVsPhi[istation], 1., 1.);
    hMeanPhiDiffVsPhi[istation]->GetXaxis()->SetTitleOffset(1.2);
    hMeanPhiDiffVsPhi[istation]->GetYaxis()->SetTitleOffset(1.7);
    hMeanPhiDiffVsPhi[istation]->GetXaxis()->SetTitle("#phi");
    hMeanPhiDiffVsPhi[istation]->GetYaxis()->SetTitle("#LT#phi_rec-#phi_sim#GT (mrad)");
    hMeanPhiDiffVsPhi[istation]->SetMaximum(5.);
    pad[page]->cd(istation+1);  hMeanPhiDiffVsPhi[istation]->Draw();
  }
  page++;  c1->Update();

  c1->Clear();  c1->cd(0);
  title = new TPaveLabel(0.1, 0.94, 0.9, 0.98,
                         "#phi_rec-#phi_sim (mrad) vs halfstrip #");
  title->SetFillColor(10);  title->Draw();
  sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
  gStyle->SetOptStat(0);
  pad[page]->Draw();
  pad[page]->Divide(2,5);
  for (int idh = 0; idh < max_idh; idh++) {
    if (!plotME1A && idh == 3) continue;
    pad[page]->cd(idh+1);  hPhiDiffVsStripCsc[idh][1]->SetMarkerSize(0.2);
    hPhiDiffVsStripCsc[idh][1]->GetXaxis()->SetTitle("Halfstrip");
    hPhiDiffVsStripCsc[idh][1]->GetXaxis()->SetTitleOffset(1.4);
    hPhiDiffVsStripCsc[idh][1]->GetYaxis()->SetTitle("#phi_rec-#phi_sim (mrad)");
    hPhiDiffVsStripCsc[idh][1]->GetXaxis()->SetLabelSize(0.07); // default=0.04
    hPhiDiffVsStripCsc[idh][1]->GetYaxis()->SetLabelSize(0.07);
    hPhiDiffVsStripCsc[idh][1]->Draw();
  }
  page++;  c1->Update();

  c1->Clear();  c1->cd(0);
  title = new TPaveLabel(0.1, 0.94, 0.9, 0.98,
                         "#phi_rec-#phi_sim, halfstrips only, different patterns");
  title->SetFillColor(10);  title->Draw();
  sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
  gStyle->SetOptStat(111110);
  pad[page]->Draw();
  pad[page]->Divide(3,3);
  int min_pattern, max_pattern;
  min_pattern = 2;
  max_pattern = CSCConstants::NUM_CLCT_PATTERNS;
  for (int idh = min_pattern; idh < max_pattern; idh++) {
    hPhiDiffPattern[idh]->GetXaxis()->SetTitle("Halfstrip");
    hPhiDiffPattern[idh]->GetXaxis()->SetTitleOffset(1.2);
    pad[page]->cd(idh-min_pattern+1);  hPhiDiffPattern[idh]->Draw();
    // if (hPhiDiffPattern[idh]->GetEntries() > 1)
    //   hPhiDiffPattern[idh]->Fit("gaus","Q");
  }
  page++;  c1->Update();

  ps->NewPage();
  c1->Clear();  c1->cd(0);
  title = new TPaveLabel(0.1, 0.94, 0.9, 0.98,
                         "#phi_1-#phi_6 (mrad), muon SimHits");
  title->SetFillColor(10);  title->Draw();
  sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
  pad[page]->Draw();
  pad[page]->Divide(2,5);
  for (int idh = 0; idh < max_idh; idh++) {
    if (!plotME1A && idh == 3) continue;
    hTrueBendCsc[idh]->GetYaxis()->SetTitle("Entries");
    hTrueBendCsc[idh]->GetYaxis()->SetTitleSize(0.07);
    hTrueBendCsc[idh]->GetXaxis()->SetLabelSize(0.07); // default=0.04
    hTrueBendCsc[idh]->GetYaxis()->SetLabelSize(0.07);
    pad[page]->cd(idh+1);  hTrueBendCsc[idh]->Draw();
  }
  page++;  c1->Update();

  ps->Close();
  delete c1;
}

void CSCTriggerPrimitivesReader::drawEfficHistos() {
  TCanvas *c1 = new TCanvas("c1", "", 0, 0, 500, 700);
  string fname = resultsFileNamesPrefix_+"lcts_effic.ps";
  TPostScript *ps = new TPostScript(fname.c_str(), 111);

  TPad *pad[MAXPAGES];
  for (int i_page = 0; i_page < MAXPAGES; i_page++) {
    pad[i_page] = new TPad("", "", .05, .05, .93, .93);
  }

  int page = 1;
  TText t;
  t.SetTextFont(32);
  t.SetTextSize(0.025);
  char pagenum[7];
  TPaveLabel *title;
  char histtitle[60];

  gStyle->SetOptDate(0);
  gStyle->SetTitleSize(0.1, "");   // size for pad title; default is 0.02

  int max_idh = plotME42 ? CSC_TYPES : CSC_TYPES-1;

  ps->NewPage();
  c1->Clear();  c1->cd(0);
  title = new TPaveLabel(0.1, 0.94, 0.9, 0.98, "ALCT efficiency vs #eta");
  title->SetFillColor(10);  title->Draw();
  sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
  gStyle->SetOptStat(0);
  pad[page]->Draw();
  pad[page]->Divide(2,2);
  TH1F *hALCTEffVsEta[MAX_STATIONS];
  for (int istation = 0; istation < MAX_STATIONS; istation++) {
    hALCTEffVsEta[istation] = (TH1F*)hEfficHitsEta[istation]->Clone();
    hALCTEffVsEta[istation]->Divide(hEfficALCTEta[istation],
                                    hEfficHitsEta[istation], 1., 1., "B");
    hALCTEffVsEta[istation]->GetXaxis()->SetTitleOffset(1.2);
    hALCTEffVsEta[istation]->GetXaxis()->SetTitle("#eta");
    hALCTEffVsEta[istation]->SetMaximum(1.05);
    sprintf(histtitle, "ALCT efficiency vs #eta, station %d", istation+1);
    hALCTEffVsEta[istation]->SetTitle(histtitle);
    pad[page]->cd(istation+1);  hALCTEffVsEta[istation]->Draw();
  }
  page++;  c1->Update();

  c1->Clear();  c1->cd(0);
  title = new TPaveLabel(0.1, 0.94, 0.9, 0.98, "ALCT efficiency vs #eta");
  title->SetFillColor(10);  title->Draw();
  sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
  gStyle->SetOptStat(11111);
  pad[page]->Draw();
  pad[page]->Divide(2,5);
  TH1F *hALCTEffVsEtaCsc[CSC_TYPES];
  for (int idh = 0; idh < max_idh; idh++) {
    if (!plotME1A && idh == 3) continue;
    hALCTEffVsEtaCsc[idh] = (TH1F*)hEfficHitsEtaCsc[idh]->Clone();
    hALCTEffVsEtaCsc[idh]->Divide(hEfficALCTEtaCsc[idh],
                                  hEfficHitsEtaCsc[idh], 1., 1., "B");
    if (idh == 3 || idh == 4 || idh == 6 || idh == 8) {
      gPad->Update();  gStyle->SetStatX(0.43);
    }
    else {
      gPad->Update();  gStyle->SetStatX(1.00);
    }
    hALCTEffVsEtaCsc[idh]->GetXaxis()->SetTitle("#eta");
    hALCTEffVsEtaCsc[idh]->GetXaxis()->SetTitleOffset(0.8);
    hALCTEffVsEtaCsc[idh]->GetXaxis()->SetTitleSize(0.07); // default=0.05
    hALCTEffVsEtaCsc[idh]->GetXaxis()->SetLabelSize(0.10); // default=0.04
    hALCTEffVsEtaCsc[idh]->GetYaxis()->SetLabelSize(0.10);
    hALCTEffVsEtaCsc[idh]->SetLabelOffset(0.012, "XY");
    hALCTEffVsEtaCsc[idh]->SetMinimum(0.50);
    hALCTEffVsEtaCsc[idh]->SetMaximum(1.05);
    hALCTEffVsEtaCsc[idh]->SetTitle(csc_type[idh].c_str());
    hALCTEffVsEtaCsc[idh]->SetTitleSize(0.1, "");
    hALCTEffVsEtaCsc[idh]->SetLineWidth(2);
    hALCTEffVsEtaCsc[idh]->SetLineColor(4);
    pad[page]->cd(idh+1);  gPad->SetGrid(1);  hALCTEffVsEtaCsc[idh]->Draw();
  }
  page++;  c1->Update();

  ps->NewPage();
  c1->Clear();  c1->cd(0);
  title = new TPaveLabel(0.1, 0.94, 0.9, 0.98, "CLCT efficiency vs #eta");
  title->SetFillColor(10);  title->Draw();
  sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
  gStyle->SetOptStat(0);
  pad[page]->Draw();
  pad[page]->Divide(2,2);
  TH1F *hCLCTEffVsEta[MAX_STATIONS];
  for (int istation = 0; istation < MAX_STATIONS; istation++) {
    hCLCTEffVsEta[istation] = (TH1F*)hEfficHitsEta[istation]->Clone();
    hCLCTEffVsEta[istation]->Divide(hEfficCLCTEta[istation],
                                    hEfficHitsEta[istation], 1., 1., "B");
    hCLCTEffVsEta[istation]->GetXaxis()->SetTitleOffset(1.2);
    hCLCTEffVsEta[istation]->GetXaxis()->SetTitle("#eta");
    hCLCTEffVsEta[istation]->SetMaximum(1.05);
    sprintf(histtitle, "CLCT efficiency vs #eta, station %d", istation+1);
    hCLCTEffVsEta[istation]->SetTitle(histtitle);
    pad[page]->cd(istation+1);  hCLCTEffVsEta[istation]->Draw();
  }
  page++;  c1->Update();

  c1->Clear();  c1->cd(0);
  title = new TPaveLabel(0.1, 0.94, 0.9, 0.98, "CLCT efficiency vs #eta");
  title->SetFillColor(10);  title->Draw();
  sprintf(pagenum, "- %d -", page);  t.DrawText(0.9, 0.02, pagenum);
  gStyle->SetOptStat(111110);
  pad[page]->Draw();
  pad[page]->Divide(2,5);
  TH1F *hCLCTEffVsEtaCsc[CSC_TYPES];
  for (int idh = 0; idh < max_idh; idh++) {
    if (!plotME1A && idh == 3) continue;
    hCLCTEffVsEtaCsc[idh] = (TH1F*)hEfficHitsEtaCsc[idh]->Clone();
    hCLCTEffVsEtaCsc[idh]->Divide(hEfficCLCTEtaCsc[idh],
                                  hEfficHitsEtaCsc[idh], 1., 1., "B");
    if (idh == 3 || idh == 4 || idh == 6 || idh == 8) {
      gPad->Update();  gStyle->SetStatX(0.43);
    }
    else {
      gPad->Update();  gStyle->SetStatX(1.00);
    }
    hCLCTEffVsEtaCsc[idh]->GetXaxis()->SetTitle("#eta");
    hCLCTEffVsEtaCsc[idh]->GetXaxis()->SetTitleOffset(0.8);
    hCLCTEffVsEtaCsc[idh]->GetXaxis()->SetTitleSize(0.07); // default=0.05
    hCLCTEffVsEtaCsc[idh]->GetXaxis()->SetLabelSize(0.10); // default=0.04
    hCLCTEffVsEtaCsc[idh]->GetYaxis()->SetLabelSize(0.10);
    hCLCTEffVsEtaCsc[idh]->SetLabelOffset(0.012, "XY");
    hCLCTEffVsEtaCsc[idh]->SetMinimum(0.50);
    hCLCTEffVsEtaCsc[idh]->SetMaximum(1.05);
    hCLCTEffVsEtaCsc[idh]->SetTitle(csc_type[idh].c_str());
    hCLCTEffVsEtaCsc[idh]->SetLineWidth(2);
    hCLCTEffVsEtaCsc[idh]->SetLineColor(4);
    pad[page]->cd(idh+1);  gPad->SetGrid(1);  hCLCTEffVsEtaCsc[idh]->Draw();
  }
  page++;  c1->Update();

  ps->Close();
  delete c1;
}

void CSCTriggerPrimitivesReader::drawHistosForTalks() {
  TCanvas *c1 = new TCanvas("c1", "", 0, 0, 500, 640);
  TCanvas *c2 = new TCanvas("c2", "", 0, 0, 540, 540);

  TPad *pad[MAXPAGES];
  for (int i_page = 0; i_page < MAXPAGES; i_page++) {
    pad[i_page] = new TPad("", "", .07, .07, .93, .93);
  }

  int page = 1;
  TPaveLabel *title;
  gStyle->SetOptDate(0);

  int max_idh = plotME42 ? CSC_TYPES : CSC_TYPES-1;

  TPostScript *eps1 = new TPostScript("clcts.eps", 113);
  eps1->NewPage();
  c1->Clear();  c1->cd(0);
  title = new TPaveLabel(0.1, 0.94, 0.9, 0.98, "CLCT quantities");
  title->SetFillColor(10);  title->Draw();
  pad[page]->Draw();
  pad[page]->Divide(2,3);
  pad[page]->cd(1);  hClctQuality->Draw();
  pad[page]->cd(2);  hClctSign->Draw();
  TH1F* hClctPatternTot = (TH1F*)hClctPattern[0]->Clone();
  hClctPatternTot->SetTitle("CLCT pattern #");
  hClctPatternTot->Add(hClctPattern[0], hClctPattern[1], 1., 1.);
  pad[page]->cd(3);  hClctPatternTot->Draw();
  hClctPattern[0]->SetLineStyle(2);  hClctPattern[0]->Draw("same");
  hClctPattern[1]->SetLineStyle(3);  hClctPattern[1]->Draw("same");
  pad[page]->cd(4);  hClctCFEB->Draw();
  pad[page]->cd(5);  hClctStripType->Draw();
  TH1F* hClctKeyStripTot = (TH1F*)hClctKeyStrip[0]->Clone();
  hClctKeyStripTot->SetTitle("CLCT key strip #");
  hClctKeyStripTot->Add(hClctKeyStrip[0], hClctKeyStrip[1], 1., 1.);
  pad[page]->cd(6);  hClctKeyStripTot->Draw();
  page++;  c1->Update();
  c1->Print("asdf.png");
  eps1->Close();

  // Resolution histograms.
  if (bookedResolHistos) {
    gStyle->SetTitleSize(0.055, "");   // size for pad title; default is 0.02

    TPostScript *eps2 = new TPostScript("alct_deltaWG.eps", 113);
    eps2->NewPage();
    c2->Clear();  c2->cd(0);
    gStyle->SetOptStat(0);
    pad[page]->Draw();
    pad[page]->Divide(1,1);
    hResolDeltaWG->GetXaxis()->SetTitle("WG_{rec} - WG_{sim}");
    hResolDeltaWG->GetXaxis()->SetTitleOffset(1.2);
    hResolDeltaWG->GetYaxis()->SetTitle("Entries");
    hResolDeltaWG->GetYaxis()->SetTitleOffset(1.9);
    hResolDeltaWG->GetXaxis()->SetLabelSize(0.04);
    hResolDeltaWG->GetYaxis()->SetLabelSize(0.04);
    pad[page]->cd(1);  hResolDeltaWG->Draw();
    page++;  c2->Update();
    eps2->Close();

    TPostScript *eps3 = new TPostScript("clct_deltaHS.eps", 113);
    eps3->NewPage();
    c2->Clear();  c2->cd(0);
    pad[page]->Draw();
    pad[page]->Divide(1,1);
    hResolDeltaHS->GetXaxis()->SetTitle("HS_{rec} - HS_{sim}");
    hResolDeltaHS->GetXaxis()->SetTitleOffset(1.2);
    hResolDeltaHS->GetYaxis()->SetTitle("Entries");
    hResolDeltaHS->GetYaxis()->SetTitleOffset(1.7);
    hResolDeltaHS->GetXaxis()->SetLabelSize(0.04); // default=0.04
    hResolDeltaHS->GetYaxis()->SetLabelSize(0.04);
    pad[page]->cd(1);  hResolDeltaHS->Draw();
    page++;  c2->Update();
    eps3->Close();

    TPostScript *eps4 = new TPostScript("clct_deltaDS.eps", 113);
    eps4->NewPage();
    c2->Clear();  c2->cd(0);
    gStyle->SetOptStat(0);
    pad[page]->Draw();
    pad[page]->Divide(1,1);
    hResolDeltaDS->GetXaxis()->SetTitle("DS_{rec} - DS_{sim}");
    hResolDeltaDS->GetXaxis()->SetTitleOffset(1.2);
    hResolDeltaDS->GetYaxis()->SetTitle("Entries");
    hResolDeltaDS->GetYaxis()->SetTitleOffset(1.6);
    hResolDeltaDS->GetXaxis()->SetLabelSize(0.04); // default=0.04
    hResolDeltaDS->GetYaxis()->SetLabelSize(0.04);
    pad[page]->cd(1);  hResolDeltaDS->Draw();
    page++;  c2->Update();
    eps4->Close();

    TPostScript *eps5 = new TPostScript("clct_deltaPhi_hs.eps", 113);
    eps5->NewPage();
    c1->Clear();  c1->cd(0);
    title = new TPaveLabel(0.1, 0.94, 0.9, 0.98,
                           "#phi_rec-#phi_sim (mrad), halfstrips only");
    title->SetFillColor(10);  title->Draw();
    gStyle->SetTitleSize(0.1, "");   // size for pad title; default is 0.02
    gStyle->SetOptStat(111110);
    pad[page]->Draw();
    pad[page]->Divide(2,5);
    for (int idh = 0; idh < max_idh; idh++) {
      if (!plotME1A && idh == 3) continue;
      pad[page]->cd(idh+1);  hPhiDiffCsc[idh][4]->Draw();
      //if (hPhiDiffCsc[idh][4]->GetEntries() > 1)
      //hPhiDiffCsc[idh][4]->Fit("gaus","Q");
      //hPhiDiffCsc[idh][4]->GetXaxis()->SetTitle("#phi_{rec} - #phi_{sim} (mrad)");
      //hPhiDiffCsc[idh][4]->GetXaxis()->SetTitleSize(0.06);
      //hPhiDiffCsc[idh][4]->GetXaxis()->SetTitleOffset(0.9);
      hPhiDiffCsc[idh][4]->GetYaxis()->SetTitle("Entries        ");
      hPhiDiffCsc[idh][4]->GetYaxis()->SetTitleSize(0.07);
      hPhiDiffCsc[idh][4]->GetYaxis()->SetTitleOffset(1.0);
      hPhiDiffCsc[idh][4]->GetXaxis()->SetLabelSize(0.10); // default=0.04
      hPhiDiffCsc[idh][4]->GetYaxis()->SetLabelSize(0.10);
      hPhiDiffCsc[idh][4]->SetLabelOffset(0.012, "XY");
    }
    page++;  c1->Update();
    eps5->Close();
  }

  // Efficiency histograms.
  if (bookedEfficHistos) {
    TPostScript *eps6 = new TPostScript("alct_effic.eps", 113);
    eps6->NewPage();
    c1->Clear();  c1->cd(0);
    title = new TPaveLabel(0.1, 0.94, 0.9, 0.98, "ALCT efficiency vs #eta");
    title->SetFillColor(10);  title->Draw();
    gStyle->SetOptStat(0);
    pad[page]->Draw();
    pad[page]->Divide(2,5);
    TH1F *hALCTEffVsEtaCsc[CSC_TYPES];
    for (int idh = 0; idh < max_idh; idh++) {
      if (!plotME1A && idh == 3) continue;
      hALCTEffVsEtaCsc[idh] = (TH1F*)hEfficHitsEtaCsc[idh]->Clone();
      hALCTEffVsEtaCsc[idh]->Divide(hEfficALCTEtaCsc[idh],
                                    hEfficHitsEtaCsc[idh], 1., 1., "B");
      hALCTEffVsEtaCsc[idh]->GetXaxis()->SetTitle("#eta");
      hALCTEffVsEtaCsc[idh]->GetXaxis()->SetTitleOffset(0.8);
      hALCTEffVsEtaCsc[idh]->GetXaxis()->SetTitleSize(0.07); // default=0.05
      hALCTEffVsEtaCsc[idh]->GetXaxis()->SetLabelSize(0.10); // default=0.04
      hALCTEffVsEtaCsc[idh]->GetYaxis()->SetLabelSize(0.10);
      hALCTEffVsEtaCsc[idh]->SetLabelOffset(0.012, "XY");
      hALCTEffVsEtaCsc[idh]->SetMinimum(0.50);
      hALCTEffVsEtaCsc[idh]->SetMaximum(1.05);
      hALCTEffVsEtaCsc[idh]->SetTitle(csc_type[idh].c_str());
      hALCTEffVsEtaCsc[idh]->SetTitleSize(0.1, "");
      hALCTEffVsEtaCsc[idh]->SetLineWidth(2);
      hALCTEffVsEtaCsc[idh]->SetLineColor(4);
      hALCTEffVsEtaCsc[idh]->SetLineColor(4);
      pad[page]->cd(idh+1);  gPad->SetGrid(1);  hALCTEffVsEtaCsc[idh]->Draw();
    }
    page++;  c1->Update();
    eps6->Close();

    TPostScript *eps7 = new TPostScript("clct_effic.eps", 113);
    eps7->NewPage();
    c1->Clear();  c1->cd(0);
    title = new TPaveLabel(0.1, 0.94, 0.9, 0.98, "CLCT efficiency vs #eta");
    title->SetFillColor(10);  title->Draw();
    gStyle->SetOptStat(0);
    pad[page]->Draw();
    pad[page]->Divide(2,5);
    TH1F *hCLCTEffVsEtaCsc[CSC_TYPES];
    for (int idh = 0; idh < max_idh; idh++) {
      if (!plotME1A && idh == 3) continue;
      hCLCTEffVsEtaCsc[idh] = (TH1F*)hEfficHitsEtaCsc[idh]->Clone();
      hCLCTEffVsEtaCsc[idh]->Divide(hEfficCLCTEtaCsc[idh],
                                    hEfficHitsEtaCsc[idh], 1., 1., "B");
      hCLCTEffVsEtaCsc[idh]->GetXaxis()->SetTitle("#eta");
      hCLCTEffVsEtaCsc[idh]->GetXaxis()->SetTitleOffset(0.8);
      hCLCTEffVsEtaCsc[idh]->GetXaxis()->SetTitleSize(0.07); // default=0.05
      hCLCTEffVsEtaCsc[idh]->GetXaxis()->SetLabelSize(0.10); // default=0.04
      hCLCTEffVsEtaCsc[idh]->GetYaxis()->SetLabelSize(0.10);
      hCLCTEffVsEtaCsc[idh]->SetLabelOffset(0.012, "XY");
      hCLCTEffVsEtaCsc[idh]->SetMinimum(0.50);
      hCLCTEffVsEtaCsc[idh]->SetMaximum(1.05);
      hCLCTEffVsEtaCsc[idh]->SetTitle(csc_type[idh].c_str());
      hCLCTEffVsEtaCsc[idh]->SetLineWidth(2);
      hCLCTEffVsEtaCsc[idh]->SetLineColor(4);
      pad[page]->cd(idh+1);  gPad->SetGrid(1);  hCLCTEffVsEtaCsc[idh]->Draw();
    }
    page++;  c1->Update();
    eps7->Close();
  }
  delete c1;
  delete c2;
}

GlobalPoint
CSCTriggerPrimitivesReader::getGlobalPosition(unsigned int rawId, int keyWg, int keyHS) const
{
    // taken from https://github.com/cms-sw/cmssw/blob/dc9f78b6af4ad56c9342cf14041b6485a60b0691/L1Trigger/CSCTriggerPrimitives/src/CSCTriggerPrimitivesReaderGEM.cc
    CSCDetId cscId = CSCDetId(rawId);
    int ring = cscId.ring();
    if (cscId.station() == 1 and cscId.ring() == 1 and (lut_wg_vs_hs_me1b[keyWg][0] <0 || keyHS>=128)){
      ring =4;
      if (keyHS >= 128)
        keyHS = keyHS-128;
    }

    CSCDetId key_id(cscId.endcap(), cscId.station(), ring,
                  cscId.chamber(), CSCConstants::KEY_CLCT_LAYER);
    auto cscChamber = geom_->chamber(cscId);
    float fractional_strip = 0.5 * (keyHS + 1) - 0.25;
    auto layer_geo = cscChamber->layer(CSCConstants::KEY_CLCT_LAYER)->geometry();
    // LCT::getKeyWG() also starts from 0
    float wire = layer_geo->middleWireOfGroup(keyWg + 1);
    LocalPoint csc_intersect = layer_geo->intersectionOfStripAndWire(fractional_strip, wire);
    GlobalPoint csc_gp = geom_->idToDet(key_id)->surface().toGlobal(csc_intersect);
    return csc_gp;

}


bool CSCTriggerPrimitivesReader::doesALCTCrossCLCT(CSCDetId id, int key_wg, int key_hs) const
{

  bool isME11 = (id.station() ==1 and (id.ring()==1 or id.ring()==4));
  if (not isME11) return true;

  int theEndcap = (id.endcap() == 1)?1 : 2;

  if (key_hs > CSCConstants::MAX_HALF_STRIP_ME1B)
    {
      key_hs = key_hs - CSCConstants::MAX_HALF_STRIP_ME1B -1;//convert it from 128-223 -> 0-95
      if ( !gangedME1a )
        {
          // wrap around ME11 HS number for -z endcap
          if (theEndcap==2) key_hs = CSCConstants::MAX_HALF_STRIP_ME1A_UNGANGED - key_hs;
          if ( key_hs >= lut_wg_vs_hs_me1a[key_wg][0] &&
               key_hs <= lut_wg_vs_hs_me1a[key_wg][1]    ) return true;
          return false;
        }
      else
        {
          if (theEndcap==2) key_hs = CSCConstants::MAX_HALF_STRIP_ME1A_GANGED - key_hs;
          if ( key_hs >= lut_wg_vs_hs_me1ag[key_wg][0] &&
               key_hs <= lut_wg_vs_hs_me1ag[key_wg][1]    ) return true;
          return false;
        }
    }
  if ( key_hs <= CSCConstants::MAX_HALF_STRIP_ME1B)
    {
      if (theEndcap==2) key_hs = CSCConstants::MAX_HALF_STRIP_ME1B - key_hs;
      if ( key_hs >= lut_wg_vs_hs_me1b[key_wg][0] &&
           key_hs <= lut_wg_vs_hs_me1b[key_wg][1]      ) return true;
    }
  return false;
}




// Returns chamber type (0-9) according to the station and ring number
int CSCTriggerPrimitivesReader::getCSCType(const CSCDetId& id) {
  int type = -999;

  if (id.station() == 1) {
    type = (id.triggerCscId()-1)/3;
    if (id.ring() == 4) {
      type = 3;
    }
  }
  else { // stations 2-4
    type = 3 + id.ring() + 2*(id.station()-2);
  }

  assert(type >= 0 && type < CSC_TYPES);
  return type;
}

// Returns halfstrips-per-radian for different CSC types
double CSCTriggerPrimitivesReader::getHsPerRad(const int idh) {
  return (NCHAMBERS[idh]*MAX_HS[idh]/TWOPI);
}

DEFINE_FWK_MODULE(CSCTriggerPrimitivesReader);
//-------------------------------------------------