/data/refman/pasoursint/CMSSW_4_2_9_HLT1_bphpatch4/src/Validation/GlobalHits/src/GlobalHitsProdHist.cc

Go to the documentation of this file.

#include "Validation/GlobalHits/interface/GlobalHitsProdHist.h"
#include "Geometry/Records/interface/CaloGeometryRecord.h"

GlobalHitsProdHist::GlobalHitsProdHist(const edm::ParameterSet& iPSet) :
  fName(""), verbosity(0), frequency(0), vtxunit(0), 
  getAllProvenances(false), printProvenanceInfo(false), count(0)
{
  std::string MsgLoggerCat = "GlobalHitsProdHist_GlobalHitsProdHist";

  // get information from parameter set
  fName = iPSet.getUntrackedParameter<std::string>("Name");
  verbosity = iPSet.getUntrackedParameter<int>("Verbosity");
  frequency = iPSet.getUntrackedParameter<int>("Frequency");
  vtxunit = iPSet.getUntrackedParameter<int>("VtxUnit");
  edm::ParameterSet m_Prov =
    iPSet.getParameter<edm::ParameterSet>("ProvenanceLookup");
  getAllProvenances = 
    m_Prov.getUntrackedParameter<bool>("GetAllProvenances");
  printProvenanceInfo = 
    m_Prov.getUntrackedParameter<bool>("PrintProvenanceInfo");

  //get Labels to use to extract information
  PxlBrlLowSrc_ = iPSet.getParameter<edm::InputTag>("PxlBrlLowSrc");
  PxlBrlHighSrc_ = iPSet.getParameter<edm::InputTag>("PxlBrlHighSrc");
  PxlFwdLowSrc_ = iPSet.getParameter<edm::InputTag>("PxlFwdLowSrc");
  PxlFwdHighSrc_ = iPSet.getParameter<edm::InputTag>("PxlFwdHighSrc");

  SiTIBLowSrc_ = iPSet.getParameter<edm::InputTag>("SiTIBLowSrc");
  SiTIBHighSrc_ = iPSet.getParameter<edm::InputTag>("SiTIBHighSrc");
  SiTOBLowSrc_ = iPSet.getParameter<edm::InputTag>("SiTOBLowSrc");
  SiTOBHighSrc_ = iPSet.getParameter<edm::InputTag>("SiTOBHighSrc");
  SiTIDLowSrc_ = iPSet.getParameter<edm::InputTag>("SiTIDLowSrc");
  SiTIDHighSrc_ = iPSet.getParameter<edm::InputTag>("SiTIDHighSrc");
  SiTECLowSrc_ = iPSet.getParameter<edm::InputTag>("SiTECLowSrc");
  SiTECHighSrc_ = iPSet.getParameter<edm::InputTag>("SiTECHighSrc");

  MuonCscSrc_ = iPSet.getParameter<edm::InputTag>("MuonCscSrc");
  MuonDtSrc_ = iPSet.getParameter<edm::InputTag>("MuonDtSrc");
  MuonRpcSrc_ = iPSet.getParameter<edm::InputTag>("MuonRpcSrc");

  ECalEBSrc_ = iPSet.getParameter<edm::InputTag>("ECalEBSrc");
  ECalEESrc_ = iPSet.getParameter<edm::InputTag>("ECalEESrc");
  ECalESSrc_ = iPSet.getParameter<edm::InputTag>("ECalESSrc");

  HCalSrc_ = iPSet.getParameter<edm::InputTag>("HCalSrc");

  // use value of first digit to determine default output level (inclusive)
  // 0 is none, 1 is basic, 2 is fill output, 3 is gather output
  verbosity %= 10;

  // print out Parameter Set information being used
  if (verbosity >= 0) {
    edm::LogInfo(MsgLoggerCat) 
      << "\n===============================\n"
      << "Initialized as EDProducer with parameter values:\n"
      << "    Name          = " << fName << "\n"
      << "    Verbosity     = " << verbosity << "\n"
      << "    Frequency     = " << frequency << "\n"
      << "    VtxUnit       = " << vtxunit << "\n"
      << "    GetProv       = " << getAllProvenances << "\n"
      << "    PrintProv     = " << printProvenanceInfo << "\n"
      << "    PxlBrlLowSrc  = " << PxlBrlLowSrc_.label() 
      << ":" << PxlBrlLowSrc_.instance() << "\n"
      << "    PxlBrlHighSrc = " << PxlBrlHighSrc_.label() 
      << ":" << PxlBrlHighSrc_.instance() << "\n"
      << "    PxlFwdLowSrc  = " << PxlFwdLowSrc_.label() 
      << ":" << PxlBrlLowSrc_.instance() << "\n"
      << "    PxlFwdHighSrc = " << PxlFwdHighSrc_.label() 
      << ":" << PxlBrlHighSrc_.instance() << "\n"
      << "    SiTIBLowSrc   = " << SiTIBLowSrc_.label() 
      << ":" << SiTIBLowSrc_.instance() << "\n"
      << "    SiTIBHighSrc  = " << SiTIBHighSrc_.label() 
      << ":" << SiTIBHighSrc_.instance() << "\n"
      << "    SiTOBLowSrc   = " << SiTOBLowSrc_.label() 
      << ":" << SiTOBLowSrc_.instance() << "\n"
      << "    SiTOBHighSrc  = " << SiTOBHighSrc_.label() 
      << ":" << SiTOBHighSrc_.instance() << "\n"
      << "    SiTIDLowSrc   = " << SiTIDLowSrc_.label() 
      << ":" << SiTIDLowSrc_.instance() << "\n"
      << "    SiTIDHighSrc  = " << SiTIDHighSrc_.label() 
      << ":" << SiTIDHighSrc_.instance() << "\n"
      << "    SiTECLowSrc   = " << SiTECLowSrc_.label() 
      << ":" << SiTECLowSrc_.instance() << "\n"
      << "    SiTECHighSrc  = " << SiTECHighSrc_.label() 
      << ":" << SiTECHighSrc_.instance() << "\n"
      << "    MuonCscSrc    = " << MuonCscSrc_.label() 
      << ":" << MuonCscSrc_.instance() << "\n"
      << "    MuonDtSrc     = " << MuonDtSrc_.label() 
      << ":" << MuonDtSrc_.instance() << "\n"
      << "    MuonRpcSrc    = " << MuonRpcSrc_.label() 
      << ":" << MuonRpcSrc_.instance() << "\n"
      << "    ECalEBSrc     = " << ECalEBSrc_.label() 
      << ":" << ECalEBSrc_.instance() << "\n"
      << "    ECalEESrc     = " << ECalEESrc_.label() 
      << ":" << ECalEESrc_.instance() << "\n"
      << "    ECalESSrc     = " << ECalESSrc_.label() 
      << ":" << ECalESSrc_.instance() << "\n"
      << "    HCalSrc       = " << HCalSrc_.label() 
      << ":" << HCalSrc_.instance() << "\n"
      << "===============================\n";
  }

  //create histograms
  Char_t hname[200];
  Char_t htitle[200];
 
  // MCGeant
  sprintf(hname,"hMCRGP1");
  histName_.push_back(hname);
  sprintf(htitle,"RawGenParticles");
  hMCRGP[0] = new TH1F(hname,htitle,100,0.,5000.);
  sprintf(hname,"hMCRGP2");
  histName_.push_back(hname);
  hMCRGP[1] = new TH1F(hname,htitle,100,0.,500.);  
  for (Int_t i = 0; i < 2; ++i) {
    hMCRGP[i]->GetXaxis()->SetTitle("Number of Raw Generated Particles");
    hMCRGP[i]->GetYaxis()->SetTitle("Count");
    histMap_[hMCRGP[i]->GetName()] = hMCRGP[i];
  }
  
  sprintf(hname,"hMCG4Vtx1");
  histName_.push_back(hname);
  sprintf(htitle,"G4 Vertices");
  hMCG4Vtx[0] = new TH1F(hname,htitle,100,0.,50000.);
  sprintf(hname,"hMCG4Vtx2");
  histName_.push_back(hname);
  hMCG4Vtx[1] = new TH1F(hname,htitle,100,-0.5,99.5); 
  for (Int_t i = 0; i < 2; ++i) {
    hMCG4Vtx[i]->GetXaxis()->SetTitle("Number of Vertices");
    hMCG4Vtx[i]->GetYaxis()->SetTitle("Count");
    histMap_[hMCG4Vtx[i]->GetName()] = hMCG4Vtx[i];
  }
  
  sprintf(hname,"hMCG4Trk1");
  histName_.push_back(hname);
  sprintf(htitle,"G4 Tracks");
  hMCG4Trk[0] = new TH1F(hname,htitle,150,0.,15000.);
  sprintf(hname,"hMCG4Trk2");
  histName_.push_back(hname);
  hMCG4Trk[1] = new TH1F(hname,htitle,150,-0.5,99.5);    
  for (Int_t i = 0; i < 2; ++i) {
    hMCG4Trk[i]->GetXaxis()->SetTitle("Number of Tracks");
    hMCG4Trk[i]->GetYaxis()->SetTitle("Count");
    histMap_[hMCG4Trk[i]->GetName()] = hMCG4Trk[i];
  }
  
  sprintf(hname,"hGeantVtxX1");
  histName_.push_back(hname);
  sprintf(htitle,"Geant vertex x/micrometer");
  hGeantVtxX[0] = new TH1F(hname,htitle,100,-8000000.,8000000.);
  sprintf(hname,"hGeantVtxX2");
  histName_.push_back(hname);
  hGeantVtxX[1] = new TH1F(hname,htitle,100,-50.,50.); 
  for (Int_t i = 0; i < 2; ++i) {
    hGeantVtxX[i]->GetXaxis()->SetTitle("x of Vertex (um)");
    hGeantVtxX[i]->GetYaxis()->SetTitle("Count");
    histMap_[hGeantVtxX[i]->GetName()] = hGeantVtxX[i];
  }
  
  sprintf(hname,"hGeantVtxY1");
  histName_.push_back(hname);
  sprintf(htitle,"Geant vertex y/micrometer");
  hGeantVtxY[0] = new TH1F(hname,htitle,100,-8000000,8000000.);
  sprintf(hname,"hGeantVtxY2");
  histName_.push_back(hname);
  hGeantVtxY[1] = new TH1F(hname,htitle,100,-50.,50.); 
  for (Int_t i = 0; i < 2; ++i) {
    hGeantVtxY[i]->GetXaxis()->SetTitle("y of Vertex (um)");
    hGeantVtxY[i]->GetYaxis()->SetTitle("Count");
    histMap_[hGeantVtxY[i]->GetName()] = hGeantVtxY[i];
  }
  
  sprintf(hname,"hGeantVtxZ1");
  histName_.push_back(hname);
  sprintf(htitle,"Geant vertex z/millimeter");
  hGeantVtxZ[0] = new TH1F(hname,htitle,100,-11000.,11000.);
  sprintf(hname,"hGeantVtxZ2");
  histName_.push_back(hname);
  hGeantVtxZ[1] = new TH1F(hname,htitle,100,-250.,250.);
  for (Int_t i = 0; i < 2; ++i) {
    hGeantVtxZ[i]->GetXaxis()->SetTitle("z of Vertex (mm)");
    hGeantVtxZ[i]->GetYaxis()->SetTitle("Count");
    histMap_[hGeantVtxZ[i]->GetName()] = hGeantVtxZ[i];
  }
  
  sprintf(hname,"hGeantTrkPt");
  histName_.push_back(hname);
  sprintf(htitle,"Geant track pt/GeV");
  hGeantTrkPt = new TH1F(hname,htitle,100,0.,200.);
  hGeantTrkPt->GetXaxis()->SetTitle("pT of Track (GeV)");
  hGeantTrkPt->GetYaxis()->SetTitle("Count");
  histMap_[hGeantTrkPt->GetName()] = hGeantTrkPt;
  
  sprintf(hname,"hGeantTrkE");
  histName_.push_back(hname);
  sprintf(htitle,"Geant track E/GeV");
  hGeantTrkE = new TH1F(hname,htitle,100,0.,5000.);
  hGeantTrkE->GetXaxis()->SetTitle("E of Track (GeV)");
  hGeantTrkE->GetYaxis()->SetTitle("Count");
  histMap_[hGeantTrkE->GetName()] = hGeantTrkE;
 
  // ECal
  sprintf(hname,"hCaloEcal1");
  histName_.push_back(hname);
  sprintf(htitle,"Ecal hits");
  hCaloEcal[0] = new TH1F(hname,htitle,100,0.,10000.);
  sprintf(hname,"hCaloEcal2");
  histName_.push_back(hname);
  hCaloEcal[1] = new TH1F(hname,htitle,100,-0.5,99.5);
  
  sprintf(hname,"hCaloEcalE1");
  histName_.push_back(hname);
  sprintf(htitle,"Ecal hits, energy/GeV");
  hCaloEcalE[0] = new TH1F(hname,htitle,100,0.,10.);
  sprintf(hname,"hCaloEcalE2");
  histName_.push_back(hname);
  hCaloEcalE[1] = new TH1F(hname,htitle,100,0.,0.1);
  
  sprintf(hname,"hCaloEcalToF1");
  histName_.push_back(hname);
  sprintf(htitle,"Ecal hits, ToF/ns");
  hCaloEcalToF[0] = new TH1F(hname,htitle,100,0.,1000.);
  sprintf(hname,"hCaloEcalToF2");
  histName_.push_back(hname);
  hCaloEcalToF[1] = new TH1F(hname,htitle,100,0.,100.);
  
  for (Int_t i = 0; i < 2; ++i) {
    hCaloEcal[i]->GetXaxis()->SetTitle("Number of Hits");
    hCaloEcal[i]->GetYaxis()->SetTitle("Count");
    histMap_[hCaloEcal[i]->GetName()] = hCaloEcal[i];
    hCaloEcalE[i]->GetXaxis()->SetTitle("Energy of Hits (GeV)");
    hCaloEcalE[i]->GetYaxis()->SetTitle("Count");
    histMap_[hCaloEcalE[i]->GetName()] = hCaloEcalE[i];
    hCaloEcalToF[i]->GetXaxis()->SetTitle("Time of Flight of Hits (ns)");
    hCaloEcalToF[i]->GetYaxis()->SetTitle("Count");
    histMap_[hCaloEcalToF[i]->GetName()] = hCaloEcalToF[i];
  }
  
  sprintf(hname,"hCaloEcalPhi");
  histName_.push_back(hname);
  sprintf(htitle,"Ecal hits, phi/rad");
  hCaloEcalPhi = new TH1F(hname,htitle,100,-3.2,3.2);
  hCaloEcalPhi->GetXaxis()->SetTitle("Phi of Hits (rad)");
  hCaloEcalPhi->GetYaxis()->SetTitle("Count");
  histMap_[hCaloEcalPhi->GetName()] = hCaloEcalPhi;
  
  sprintf(hname,"hCaloEcalEta");
  histName_.push_back(hname);
  sprintf(htitle,"Ecal hits, eta");
  hCaloEcalEta = new TH1F(hname,htitle,100,-5.5,5.5);
  hCaloEcalEta->GetXaxis()->SetTitle("Eta of Hits");
  hCaloEcalEta->GetYaxis()->SetTitle("Count");
  histMap_[hCaloEcalEta->GetName()] = hCaloEcalEta;

  sprintf(hname,"hCaloPreSh1");
  histName_.push_back(hname);
  sprintf(htitle,"PreSh hits");
  hCaloPreSh[0] = new TH1F(hname,htitle,100,0.,10000.);
  sprintf(hname,"hCaloPreSh2");
  histName_.push_back(hname);
  hCaloPreSh[1] = new TH1F(hname,htitle,100,-0.5,99.5);
  
  sprintf(hname,"hCaloPreShE1");
  histName_.push_back(hname);
  sprintf(htitle,"PreSh hits, energy/GeV");
  hCaloPreShE[0] = new TH1F(hname,htitle,100,0.,10.);
  sprintf(hname,"hCaloPreShE2");
  histName_.push_back(hname);
  hCaloPreShE[1] = new TH1F(hname,htitle,100,0.,0.1);

  sprintf(hname,"hCaloPreShToF1");
  histName_.push_back(hname);
  sprintf(htitle,"PreSh hits, ToF/ns");
  hCaloPreShToF[0] = new TH1F(hname,htitle,100,0.,1000.);
  sprintf(hname,"hCaloPreShToF2");
  histName_.push_back(hname);
  hCaloPreShToF[1] = new TH1F(hname,htitle,100,0.,100.);
  
  for (Int_t i = 0; i < 2; ++i) {
    hCaloPreSh[i]->GetXaxis()->SetTitle("Number of Hits");
    hCaloPreSh[i]->GetYaxis()->SetTitle("Count");
    histMap_[hCaloPreSh[i]->GetName()] = hCaloPreSh[i];
    hCaloPreShE[i]->GetXaxis()->SetTitle("Energy of Hits (GeV)");
    hCaloPreShE[i]->GetYaxis()->SetTitle("Count");
    histMap_[hCaloPreShE[i]->GetName()] = hCaloPreShE[i];
    hCaloPreShToF[i]->GetXaxis()->SetTitle("Time of Flight of Hits (ns)");
    hCaloPreShToF[i]->GetYaxis()->SetTitle("Count");
    histMap_[hCaloPreShToF[i]->GetName()] = hCaloPreShToF[i];
  }
  
  sprintf(hname,"hCaloPreShPhi");
  histName_.push_back(hname);
  sprintf(htitle,"PreSh hits, phi/rad");
  hCaloPreShPhi = new TH1F(hname,htitle,100,-3.2,3.2);
  hCaloPreShPhi->GetXaxis()->SetTitle("Phi of Hits (rad)");
  hCaloPreShPhi->GetYaxis()->SetTitle("Count");
  histMap_[hCaloPreShPhi->GetName()] = hCaloPreShPhi;

  sprintf(hname,"hCaloPreShEta");
  histName_.push_back(hname);
  sprintf(htitle,"PreSh hits, eta");
  hCaloPreShEta = new TH1F(hname,htitle,100,-5.5,5.5);
  hCaloPreShEta->GetXaxis()->SetTitle("Eta of Hits");
  hCaloPreShEta->GetYaxis()->SetTitle("Count");
  histMap_[hCaloPreShEta->GetName()] = hCaloPreShEta;  

  // Hcal
  sprintf(hname,"hCaloHcal1");
  histName_.push_back(hname);
  sprintf(htitle,"Hcal hits");
  hCaloHcal[0] = new TH1F(hname,htitle,100,0.,10000.);
  sprintf(hname,"hCaloHcal2");
  histName_.push_back(hname);
  hCaloHcal[1] = new TH1F(hname,htitle,100,-0.5,99.5);
  
  sprintf(hname,"hCaloHcalE1");
  histName_.push_back(hname);
  sprintf(htitle,"Hcal hits, energy/GeV");
  hCaloHcalE[0] = new TH1F(hname,htitle,100,0.,10.);
  sprintf(hname,"hCaloHcalE2");
  histName_.push_back(hname);
  hCaloHcalE[1] = new TH1F(hname,htitle,100,0.,0.1);
  
  sprintf(hname,"hCaloHcalToF1");
  histName_.push_back(hname);
  sprintf(htitle,"Hcal hits, ToF/ns");
  hCaloHcalToF[0] = new TH1F(hname,htitle,100,0.,1000.);
  sprintf(hname,"hCaloHcalToF2");
  histName_.push_back(hname);
  hCaloHcalToF[1] = new TH1F(hname,htitle,100,0.,100.);
  
  for (Int_t i = 0; i < 2; ++i) {
    hCaloHcal[i]->GetXaxis()->SetTitle("Number of Hits");
    hCaloHcal[i]->GetYaxis()->SetTitle("Count");
    histMap_[hCaloHcal[i]->GetName()] = hCaloHcal[i];
    hCaloHcalE[i]->GetXaxis()->SetTitle("Energy of Hits (GeV)");
    hCaloHcalE[i]->GetYaxis()->SetTitle("Count");
    histMap_[hCaloHcalE[i]->GetName()] = hCaloHcalE[i];
    hCaloHcalToF[i]->GetXaxis()->SetTitle("Time of Flight of Hits (ns)");
    hCaloHcalToF[i]->GetYaxis()->SetTitle("Count");
    histMap_[hCaloHcalToF[i]->GetName()] = hCaloHcalToF[i];
  }
  
  sprintf(hname,"hCaloHcalPhi");
  histName_.push_back(hname);
  sprintf(htitle,"Hcal hits, phi/rad");
  hCaloHcalPhi = new TH1F(hname,htitle,100,-3.2,3.2);
  hCaloHcalPhi->GetXaxis()->SetTitle("Phi of Hits (rad)");
  hCaloHcalPhi->GetYaxis()->SetTitle("Count");
  histMap_[hCaloHcalPhi->GetName()] = hCaloHcalPhi;
  
  sprintf(hname,"hCaloHcalEta");
  histName_.push_back(hname);
  sprintf(htitle,"Hcal hits, eta");
  hCaloHcalEta = new TH1F(hname,htitle,100,-5.5,5.5);
  hCaloHcalEta->GetXaxis()->SetTitle("Eta of Hits");
  hCaloHcalEta->GetYaxis()->SetTitle("Count");
  histMap_[hCaloHcalEta->GetName()] = hCaloHcalEta;

  // tracker
  sprintf(hname,"hTrackerPx1");
  histName_.push_back(hname);
  sprintf(htitle,"Pixel hits");
  hTrackerPx[0] = new TH1F(hname,htitle,100,0.,10000.);
  sprintf(hname,"hTrackerPx2");
  histName_.push_back(hname);
  hTrackerPx[1] = new TH1F(hname,htitle,100,-0.5,99.5);
  for (Int_t i = 0; i < 2; ++i) {
    hTrackerPx[i]->GetXaxis()->SetTitle("Number of Pixel Hits");
    hTrackerPx[i]->GetYaxis()->SetTitle("Count");
    histMap_[hTrackerPx[i]->GetName()] = hTrackerPx[i];
  }
  
  sprintf(hname,"hTrackerPxPhi");
  histName_.push_back(hname);
  sprintf(htitle,"Pixel hits phi/rad");
  hTrackerPxPhi = new TH1F(hname,htitle,100,-3.2,3.2);
  hTrackerPxPhi->GetXaxis()->SetTitle("Phi of Hits (rad)");
  hTrackerPxPhi->GetYaxis()->SetTitle("Count");
  histMap_[hTrackerPxPhi->GetName()] = hTrackerPxPhi;

  sprintf(hname,"hTrackerPxEta");
  histName_.push_back(hname);
  sprintf(htitle,"Pixel hits eta");
  hTrackerPxEta = new TH1F(hname,htitle,100,-3.5,3.5);
  hTrackerPxEta->GetXaxis()->SetTitle("Eta of Hits");
  hTrackerPxEta->GetYaxis()->SetTitle("Count");
  histMap_[hTrackerPxEta->GetName()] = hTrackerPxEta;

  sprintf(hname,"hTrackerPxBToF");
  histName_.push_back(hname);
  sprintf(htitle,"Pixel barrel hits, ToF/ns");
  hTrackerPxBToF = new TH1F(hname,htitle,100,0.,40.);
  hTrackerPxBToF->GetXaxis()->SetTitle("Time of Flight of Hits (ns)");
  hTrackerPxBToF->GetYaxis()->SetTitle("Count");
  histMap_[hTrackerPxBToF->GetName()] = hTrackerPxBToF;

  sprintf(hname,"hTrackerPxBR");
  histName_.push_back(hname);
  sprintf(htitle,"Pixel barrel hits, R/cm");
  hTrackerPxBR = new TH1F(hname,htitle,100,0.,50.);
  hTrackerPxBR->GetXaxis()->SetTitle("R of Hits (cm)");
  hTrackerPxBR->GetYaxis()->SetTitle("Count");
  histMap_[hTrackerPxBR->GetName()] = hTrackerPxBR;

  sprintf(hname,"hTrackerPxFToF");
  histName_.push_back(hname);
  sprintf(htitle,"Pixel forward hits, ToF/ns");
  hTrackerPxFToF = new TH1F(hname,htitle,100,0.,50.);
  hTrackerPxFToF->GetXaxis()->SetTitle("Time of Flight of Hits (ns)");
  hTrackerPxFToF->GetYaxis()->SetTitle("Count");
  histMap_[hTrackerPxFToF->GetName()] = hTrackerPxFToF;

  sprintf(hname,"hTrackerPxFZ");
  histName_.push_back(hname);
  sprintf(htitle,"Pixel forward hits, Z/cm");
  hTrackerPxFZ = new TH1F(hname,htitle,200,-100.,100.);
  hTrackerPxFZ->GetXaxis()->SetTitle("Z of Hits (cm)");
  hTrackerPxFZ->GetYaxis()->SetTitle("Count");
  histMap_[hTrackerPxFZ->GetName()] = hTrackerPxFZ;
  
  sprintf(hname,"hTrackerSi1");
  histName_.push_back(hname);
  sprintf(htitle,"Silicon hits");
  hTrackerSi[0] = new TH1F(hname,htitle,100,0.,10000.);
  sprintf(hname,"hTrackerSi2");
  histName_.push_back(hname);
  hTrackerSi[1] = new TH1F(hname,htitle,100,-0.5,99.5);
  for (Int_t i = 0; i < 2; ++i) { 
    hTrackerSi[i]->GetXaxis()->SetTitle("Number of Silicon Hits");
    hTrackerSi[i]->GetYaxis()->SetTitle("Count");
    histMap_[hTrackerSi[i]->GetName()] = hTrackerSi[i];
  }
  
  sprintf(hname,"hTrackerSiPhi");
  histName_.push_back(hname);
  sprintf(htitle,"Silicon hits phi/rad");
  hTrackerSiPhi = new TH1F(hname,htitle,100,-3.2,3.2);
  hTrackerSiPhi->GetXaxis()->SetTitle("Phi of Hits (rad)");
  hTrackerSiPhi->GetYaxis()->SetTitle("Count");
  histMap_[hTrackerSiPhi->GetName()] = hTrackerSiPhi;

  sprintf(hname,"hTrackerSiEta");
  histName_.push_back(hname);
  sprintf(htitle,"Silicon hits eta");
  hTrackerSiEta = new TH1F(hname,htitle,100,-3.5,3.5);
  hTrackerSiEta->GetXaxis()->SetTitle("Eta of Hits");
  hTrackerSiEta->GetYaxis()->SetTitle("Count");
  histMap_[hTrackerSiEta->GetName()] = hTrackerSiEta;
  
  sprintf(hname,"hTrackerSiBToF");
  histName_.push_back(hname);
  sprintf(htitle,"Silicon barrel hits, ToF/ns");
  hTrackerSiBToF = new TH1F(hname,htitle,100,0.,50.);
  hTrackerSiBToF->GetXaxis()->SetTitle("Time of Flight of Hits (ns)");
  hTrackerSiBToF->GetYaxis()->SetTitle("Count");
  histMap_[hTrackerSiBToF->GetName()] = hTrackerSiBToF;

  sprintf(hname,"hTrackerSiBR");
  histName_.push_back(hname);
  sprintf(htitle,"Silicon barrel hits, R/cm");
  hTrackerSiBR = new TH1F(hname,htitle,100,0.,200.);
  hTrackerSiBR->GetXaxis()->SetTitle("R of Hits (cm)");
  hTrackerSiBR->GetYaxis()->SetTitle("Count");
  histMap_[hTrackerSiBR->GetName()] = hTrackerSiBR;
  
  sprintf(hname,"hTrackerSiFToF");
  histName_.push_back(hname);
  sprintf(htitle,"Silicon forward hits, ToF/ns");
  hTrackerSiFToF = new TH1F(hname,htitle,100,0.,75.);
  hTrackerSiFToF->GetXaxis()->SetTitle("Time of Flight of Hits (ns)");
  hTrackerSiFToF->GetYaxis()->SetTitle("Count");
  histMap_[hTrackerSiFToF->GetName()] = hTrackerSiFToF;

  sprintf(hname,"hTrackerSiFZ");
  histName_.push_back(hname);
  sprintf(htitle,"Silicon forward hits, Z/cm");
  hTrackerSiFZ = new TH1F(hname,htitle,200,-300.,300.);
  hTrackerSiFZ->GetXaxis()->SetTitle("Z of Hits (cm)");
  hTrackerSiFZ->GetYaxis()->SetTitle("Count");
  histMap_[hTrackerSiFZ->GetName()] = hTrackerSiFZ;  

  // muon
  sprintf(hname,"hMuon1");
  histName_.push_back(hname);
  sprintf(htitle,"Muon hits");
  hMuon[0] = new TH1F(hname,htitle,100,0.,10000.);
  sprintf(hname,"hMuon2");
  histName_.push_back(hname);
  hMuon[1] = new TH1F(hname,htitle,100,-0.5,99.5);
  for (Int_t i = 0; i < 2; ++i) { 
    hMuon[i]->GetXaxis()->SetTitle("Number of Muon Hits");
    hMuon[i]->GetYaxis()->SetTitle("Count");
    histMap_[hMuon[i]->GetName()] = hMuon[i];
  }  
  
  sprintf(hname,"hMuonPhi");
  histName_.push_back(hname);
  sprintf(htitle,"Muon hits phi/rad");
  hMuonPhi = new TH1F(hname,htitle,100,-3.2,3.2);
  hMuonPhi->GetXaxis()->SetTitle("Phi of Hits (rad)");
  hMuonPhi->GetYaxis()->SetTitle("Count");
  histMap_[hMuonPhi->GetName()] = hMuonPhi;

  sprintf(hname,"hMuonEta");
  histName_.push_back(hname);
  sprintf(htitle,"Muon hits eta");
  hMuonEta = new TH1F(hname,htitle,100,-3.5,3.5);
  hMuonEta->GetXaxis()->SetTitle("Eta of Hits");
  hMuonEta->GetYaxis()->SetTitle("Count");
  histMap_[hMuonEta->GetName()] = hMuonEta;

  sprintf(hname,"hMuonCscToF1");
  histName_.push_back(hname);
  sprintf(htitle,"Muon CSC hits, ToF/ns");
  hMuonCscToF[0] = new TH1F(hname,htitle,100,0.,250.);
  sprintf(hname,"hMuonCscToF2");
  histName_.push_back(hname);
  hMuonCscToF[1] = new TH1F(hname,htitle,100,0.,50.);
  for (Int_t i = 0; i < 2; ++i) {   
    hMuonCscToF[i]->GetXaxis()->SetTitle("Time of Flight of Hits (ns)");
    hMuonCscToF[i]->GetYaxis()->SetTitle("Count");
    histMap_[hMuonCscToF[i]->GetName()] = hMuonCscToF[i];
  }  
  
  sprintf(hname,"hMuonCscZ");
  histName_.push_back(hname);
  sprintf(htitle,"Muon CSC hits, Z/cm");
  hMuonCscZ = new TH1F(hname,htitle,200,-1500.,1500.);
  hMuonCscZ->GetXaxis()->SetTitle("Z of Hits (cm)");
  hMuonCscZ->GetYaxis()->SetTitle("Count");
  histMap_[hMuonCscZ->GetName()] = hMuonCscZ;
  
  sprintf(hname,"hMuonDtToF1");
  histName_.push_back(hname);
  sprintf(htitle,"Muon DT hits, ToF/ns");
  hMuonDtToF[0] = new TH1F(hname,htitle,100,0.,250.);
  sprintf(hname,"hMuonDtToF2");
  histName_.push_back(hname);
  hMuonDtToF[1] = new TH1F(hname,htitle,100,0.,50.);
  for (Int_t i = 0; i < 2; ++i) {   
    hMuonDtToF[i]->GetXaxis()->SetTitle("Time of Flight of Hits (ns)");
    hMuonDtToF[i]->GetYaxis()->SetTitle("Count");
    histMap_[hMuonDtToF[i]->GetName()] = hMuonDtToF[i];
  } 
  
  sprintf(hname,"hMuonDtR");
  histName_.push_back(hname);
  sprintf(htitle,"Muon DT hits, R/cm");
  hMuonDtR = new TH1F(hname,htitle,100,0.,1500.); 
  hMuonDtR->GetXaxis()->SetTitle("R of Hits (cm)");
  hMuonDtR->GetYaxis()->SetTitle("Count");
  histMap_[hMuonDtR->GetName()] = hMuonDtR;

  sprintf(hname,"hMuonRpcFToF1");
  histName_.push_back(hname);
  sprintf(htitle,"Muon RPC forward hits, ToF/ns");
  hMuonRpcFToF[0] = new TH1F(hname,htitle,100,0.,250.);
  sprintf(hname,"hMuonRpcFToF2");
  histName_.push_back(hname);
  hMuonRpcFToF[1] = new TH1F(hname,htitle,100,0.,50.);
  for (Int_t i = 0; i < 2; ++i) {   
    hMuonRpcFToF[i]->GetXaxis()->SetTitle("Time of Flight of Hits (ns)");
    hMuonRpcFToF[i]->GetYaxis()->SetTitle("Count");
    histMap_[hMuonRpcFToF[i]->GetName()] = hMuonRpcFToF[i];
  }  
  
  sprintf(hname,"hMuonRpcFZ");
  histName_.push_back(hname);
  sprintf(htitle,"Muon RPC forward hits, Z/cm");
  hMuonRpcFZ = new TH1F(hname,htitle,201,-1500.,1500.);
  hMuonRpcFZ->GetXaxis()->SetTitle("Z of Hits (cm)");
  hMuonRpcFZ->GetYaxis()->SetTitle("Count");
  histMap_[hMuonRpcFZ->GetName()] = hMuonRpcFZ;

  sprintf(hname,"hMuonRpcBToF1");
  histName_.push_back(hname);
  sprintf(htitle,"Muon RPC barrel hits, ToF/ns");
  hMuonRpcBToF[0] = new TH1F(hname,htitle,100,0.,250.);
  sprintf(hname,"hMuonRpcBToF2");
  histName_.push_back(hname);
  hMuonRpcBToF[1] = new TH1F(hname,htitle,100,0.,50.);
  for (Int_t i = 0; i < 2; ++i) {   
    hMuonRpcBToF[i]->GetXaxis()->SetTitle("Time of Flight of Hits (ns)");
    hMuonRpcBToF[i]->GetYaxis()->SetTitle("Count");
    histMap_[hMuonRpcBToF[i]->GetName()] = hMuonRpcBToF[i];
  }
  
  sprintf(hname,"hMuonRpcBR");
  histName_.push_back(hname);
  sprintf(htitle,"Muon RPC barrel hits, R/cm");
  hMuonRpcBR = new TH1F(hname,htitle,100,0.,1500.);
  hMuonRpcBR->GetXaxis()->SetTitle("R of Hits (cm)");
  hMuonRpcBR->GetYaxis()->SetTitle("Count"); 
  histMap_[hMuonRpcBR->GetName()] = hMuonRpcBR;

  // create persistent objects
  for (std::size_t i = 0; i < histName_.size(); ++i) {
    produces<TH1F, edm::InRun>(histName_[i]).setBranchAlias(histName_[i]);
  }
}

GlobalHitsProdHist::~GlobalHitsProdHist() 
{
}

void GlobalHitsProdHist::beginJob( void )
{
  return;
}

void GlobalHitsProdHist::endJob()
{
  std::string MsgLoggerCat = "GlobalHitsProdHist_endJob";
  if (verbosity >= 0)
    edm::LogInfo(MsgLoggerCat) 
      << "Terminating having processed " << count << " events.";
  return;
}

void GlobalHitsProdHist::produce(edm::Event& iEvent, 
                                 const edm::EventSetup& iSetup)
{
  std::string MsgLoggerCat = "GlobalHitsProdHist_produce";

  // keep track of number of events processed
  ++count;

  // get event id information
  int nrun = iEvent.id().run();
  int nevt = iEvent.id().event();

  if (verbosity > 0) {
    edm::LogInfo(MsgLoggerCat)
      << "Processing run " << nrun << ", event " << nevt
      << " (" << count << " events total)";
  } else if (verbosity == 0) {
    if (nevt%frequency == 0 || nevt == 1) {
      edm::LogInfo(MsgLoggerCat)
        << "Processing run " << nrun << ", event " << nevt
        << " (" << count << " events total)";
    }
  }

  // look at information available in the event
  if (getAllProvenances) {

    std::vector<const edm::Provenance*> AllProv;
    iEvent.getAllProvenance(AllProv);

    if (verbosity >= 0)
      edm::LogInfo(MsgLoggerCat)
        << "Number of Provenances = " << AllProv.size();

    if (printProvenanceInfo && (verbosity >= 0)) {
      TString eventout("\nProvenance info:\n");      

      for (unsigned int i = 0; i < AllProv.size(); ++i) {
        eventout += "\n       ******************************";
        eventout += "\n       Module       : ";
        eventout += AllProv[i]->moduleLabel();
        eventout += "\n       ProductID    : ";
        eventout += AllProv[i]->productID().id();
        eventout += "\n       ClassName    : ";
        eventout += AllProv[i]->className();
        eventout += "\n       InstanceName : ";
        eventout += AllProv[i]->productInstanceName();
        eventout += "\n       BranchName   : ";
        eventout += AllProv[i]->branchName();
      }
      eventout += "\n       ******************************\n";
      edm::LogInfo(MsgLoggerCat) << eventout << "\n";
      printProvenanceInfo = false;
    }
    getAllProvenances = false;
  }

  // call fill functions
  //gather G4MC information from event
  fillG4MC(iEvent);
  // gather Tracker information from event
  fillTrk(iEvent,iSetup);
  // gather muon information from event
  fillMuon(iEvent, iSetup);
  // gather Ecal information from event
  fillECal(iEvent, iSetup);
  // gather Hcal information from event
  fillHCal(iEvent, iSetup);

  if (verbosity > 0)
    edm::LogInfo (MsgLoggerCat)
      << "Done gathering data from event.";

  return;
}

void GlobalHitsProdHist::endRun(edm::Run& iRun, const edm::EventSetup& iSetup)
{

  std::string MsgLoggerCat = "GlobalHitsProdHist_endRun";

  TString eventout;
  TString eventoutw;
  bool warning = false;

  if (verbosity > 0)
    edm::LogInfo (MsgLoggerCat)
      << "\nStoring histograms.";

  // store persistent objects
  std::map<std::string, TH1F*>::iterator iter;
  for (std::size_t i = 0; i < histName_.size(); ++i) {
    iter = histMap_.find(histName_[i]);
    if (iter != histMap_.end()) {
      std::auto_ptr<TH1F> hist1D(iter->second);
      eventout += "\n Storing histogram " + histName_[i];
      iRun.put(hist1D, histName_[i]);
    } else {
      warning = true;
      eventoutw += "\n Unable to find histogram with name " + histName_[i];
    }
  }

  if (verbosity > 0) {
    edm::LogInfo(MsgLoggerCat) << eventout << "\n";
    if (warning)
      edm::LogWarning(MsgLoggerCat) << eventoutw << "\n";
  }
  return;
}

//==================fill and store functions================================
void GlobalHitsProdHist::fillG4MC(edm::Event& iEvent)
{

  std::string MsgLoggerCat = "GlobalHitsProdHist_fillG4MC";
 
  TString eventout;
  if (verbosity > 0)
    eventout = "\nGathering info:";

  // get MC information
  edm::Handle<edm::HepMCProduct> HepMCEvt;
  std::vector<edm::Handle<edm::HepMCProduct> > AllHepMCEvt;
  iEvent.getManyByType(AllHepMCEvt);

  // loop through products and extract VtxSmearing if available. Any of them
  // should have the information needed
  for (unsigned int i = 0; i < AllHepMCEvt.size(); ++i) {
    HepMCEvt = AllHepMCEvt[i];
    if ((HepMCEvt.provenance()->product()).moduleLabel() == "VtxSmeared")
      break;
  }

  if (!HepMCEvt.isValid()) {
    edm::LogWarning(MsgLoggerCat)
      << "Unable to find HepMCProduct in event!";
    return;
  } else {
    eventout += "\n          Using HepMCProduct: ";
    eventout += (HepMCEvt.provenance()->product()).moduleLabel();
  }
  const HepMC::GenEvent* MCEvt = HepMCEvt->GetEvent();
  nRawGenPart = MCEvt->particles_size();

  if (verbosity > 1) {
    eventout += "\n          Number of Raw Particles collected:......... ";
    eventout += nRawGenPart;
  }  

  if (hMCRGP[0]) hMCRGP[0]->Fill((float)nRawGenPart);
  if (hMCRGP[1]) hMCRGP[1]->Fill((float)nRawGenPart);  

  // get G4Vertex information
  // convert unit stored in SimVertex to mm
  float unit = 0.;
  if (vtxunit == 0) unit = 1.;  // already in mm
  if (vtxunit == 1) unit = 10.; // stored in cm, convert to mm

  edm::Handle<edm::SimVertexContainer> G4VtxContainer;
  iEvent.getByType(G4VtxContainer);
  if (!G4VtxContainer.isValid()) {
    edm::LogWarning(MsgLoggerCat)
      << "Unable to find SimVertex in event!";
    return;
  }
  int i = 0;
  edm::SimVertexContainer::const_iterator itVtx;
  for (itVtx = G4VtxContainer->begin(); itVtx != G4VtxContainer->end(); 
       ++itVtx) {
    
    ++i;

    const math::XYZTLorentzVector G4Vtx1(itVtx->position().x(),
                                         itVtx->position().y(),
                                         itVtx->position().z(),
                                         itVtx->position().e());

    double G4Vtx[4];
    G4Vtx1.GetCoordinates(G4Vtx);

    if (hGeantVtxX[0]) hGeantVtxX[0]->Fill((G4Vtx[0]*unit)/micrometer);
    if (hGeantVtxX[1]) hGeantVtxX[1]->Fill((G4Vtx[0]*unit)/micrometer);
    
    if (hGeantVtxY[0]) hGeantVtxY[0]->Fill((G4Vtx[1]*unit)/micrometer);
    if (hGeantVtxY[1]) hGeantVtxY[1]->Fill((G4Vtx[1]*unit)/micrometer);
    
    if (hGeantVtxZ[0]) hGeantVtxZ[0]->Fill((G4Vtx[2]*unit)/millimeter);
    if (hGeantVtxZ[1]) hGeantVtxZ[1]->Fill((G4Vtx[2]*unit)/millimeter); 
    
  }

  if (verbosity > 1) {
    eventout += "\n          Number of G4Vertices collected:............ ";
    eventout += i;
  }  

  if (hMCG4Vtx[0]) hMCG4Vtx[0]->Fill((float)i);
  if (hMCG4Vtx[1]) hMCG4Vtx[1]->Fill((float)i);  

  // get G4Track information
  edm::Handle<edm::SimTrackContainer> G4TrkContainer;
  iEvent.getByType(G4TrkContainer);
  if (!G4TrkContainer.isValid()) {
    edm::LogWarning(MsgLoggerCat)
      << "Unable to find SimTrack in event!";
    return;
  }
  i = 0;
  edm::SimTrackContainer::const_iterator itTrk;
  for (itTrk = G4TrkContainer->begin(); itTrk != G4TrkContainer->end(); 
       ++itTrk) {

    ++i;

    const math::XYZTLorentzVector G4Trk1(itTrk->momentum().x(),
                                         itTrk->momentum().y(),
                                         itTrk->momentum().z(),
                                         itTrk->momentum().e());
    double G4Trk[4];
    G4Trk1.GetCoordinates(G4Trk);

    if (hGeantTrkPt) hGeantTrkPt->
                        Fill(sqrt(G4Trk[0]*G4Trk[0]+G4Trk[1]*G4Trk[1]));
    if (hGeantTrkE) hGeantTrkE->Fill(G4Trk[3]);
  } 

  if (verbosity > 1) {
    eventout += "\n          Number of G4Tracks collected:.............. ";
    eventout += i;
  }  

  if (hMCG4Trk[0]) hMCG4Trk[0]->Fill((float)i);
  if (hMCG4Trk[1]) hMCG4Trk[1]->Fill((float)i); 

  if (verbosity > 0)
    edm::LogInfo(MsgLoggerCat) << eventout << "\n";
    
  return;
}

void GlobalHitsProdHist::fillTrk(edm::Event& iEvent, 
                                 const edm::EventSetup& iSetup)
{

  nPxlHits = 0;
  std::string MsgLoggerCat = "GlobalHitsProdHist_fillTrk";

  TString eventout;
  if (verbosity > 0)
    eventout = "\nGathering info:";  
  
  // access the tracker geometry
  edm::ESHandle<TrackerGeometry> theTrackerGeometry;
  iSetup.get<TrackerDigiGeometryRecord>().get(theTrackerGeometry);
  if (!theTrackerGeometry.isValid()) {
    edm::LogWarning(MsgLoggerCat)
      << "Unable to find TrackerDigiGeometryRecord in event!";
    return;
  }
  const TrackerGeometry& theTracker(*theTrackerGeometry);
    
  // iterator to access containers
  edm::PSimHitContainer::const_iterator itHit;

  // get Pixel Barrel information
  edm::PSimHitContainer thePxlBrlHits;
  // extract low container
  edm::Handle<edm::PSimHitContainer> PxlBrlLowContainer;
  iEvent.getByLabel(PxlBrlLowSrc_,PxlBrlLowContainer);
  if (!PxlBrlLowContainer.isValid()) {
    edm::LogWarning(MsgLoggerCat)
      << "Unable to find TrackerHitsPixelBarrelLowTof in event!";
    return;
  }
  // extract high container
  edm::Handle<edm::PSimHitContainer> PxlBrlHighContainer;
  iEvent.getByLabel(PxlBrlHighSrc_,PxlBrlHighContainer);
  if (!PxlBrlHighContainer.isValid()) {
    edm::LogWarning(MsgLoggerCat)
      << "Unable to find TrackerHitsPixelBarrelHighTof in event!";
    return;
  }
  // place both containers into new container
  thePxlBrlHits.insert(thePxlBrlHits.end(),PxlBrlLowContainer->begin(),
                       PxlBrlLowContainer->end());
  thePxlBrlHits.insert(thePxlBrlHits.end(),PxlBrlHighContainer->begin(),
                       PxlBrlHighContainer->end());

  // cycle through new container
  int i = 0, j = 0;
  for (itHit = thePxlBrlHits.begin(); itHit != thePxlBrlHits.end(); ++itHit) {

    ++i;

    // create a DetId from the detUnitId
    DetId theDetUnitId(itHit->detUnitId());
    int detector = theDetUnitId.det();
    int subdetector = theDetUnitId.subdetId();

    // check that expected detector is returned
    if ((detector == dTrk) && (subdetector == sdPxlBrl)) {

      // get the GeomDetUnit from the geometry using theDetUnitID
      const GeomDetUnit *theDet = theTracker.idToDetUnit(theDetUnitId);

      if (!theDet) {
        edm::LogWarning(MsgLoggerCat)
          << "Unable to get GeomDetUnit from PxlBrlHits for Hit " << i;
        continue;
      }

      ++j;

      // get the Surface of the hit (knows how to go from local <-> global)
      const BoundPlane& bSurface = theDet->surface();

      if(hTrackerPxBToF) hTrackerPxBToF->Fill(itHit->tof());
      if(hTrackerPxBR) 
        hTrackerPxBR->Fill(bSurface.toGlobal(itHit->localPosition()).perp());
      if(hTrackerPxPhi) 
        hTrackerPxPhi->Fill(bSurface.toGlobal(itHit->localPosition()).phi());
      if(hTrackerPxEta) 
        hTrackerPxEta->Fill(bSurface.toGlobal(itHit->localPosition()).eta());

    } else {
      edm::LogWarning(MsgLoggerCat)
        << "PxlBrl PSimHit " << i 
        << " is expected to be (det,subdet) = (" 
        << dTrk << "," << sdPxlBrl
        << "); value returned is: ("
        << detector << "," << subdetector << ")";
      continue;
    } // end detector type check
  } // end loop through PxlBrl Hits

  if (verbosity > 1) {
    eventout += "\n          Number of Pixel Barrel Hits collected:..... ";
    eventout += j;
  }  
  
  nPxlHits += j;

  // get Pixel Forward information
  edm::PSimHitContainer thePxlFwdHits;
  // extract low container
  edm::Handle<edm::PSimHitContainer> PxlFwdLowContainer;
  iEvent.getByLabel(PxlFwdLowSrc_,PxlFwdLowContainer);
  if (!PxlFwdLowContainer.isValid()) {
    edm::LogWarning(MsgLoggerCat)
      << "Unable to find TrackerHitsPixelEndcapLowTof in event!";
    return;
  }
  // extract high container
  edm::Handle<edm::PSimHitContainer> PxlFwdHighContainer;
  iEvent.getByLabel(PxlFwdHighSrc_,PxlFwdHighContainer);
  if (!PxlFwdHighContainer.isValid()) {
    edm::LogWarning("GlobalHitsProdHist_fillTrk")
      << "Unable to find TrackerHitsPixelEndcapHighTof in event!";
    return;
  }
  // place both containers into new container
  thePxlFwdHits.insert(thePxlFwdHits.end(),PxlFwdLowContainer->begin(),
                       PxlFwdLowContainer->end());
  thePxlFwdHits.insert(thePxlFwdHits.end(),PxlFwdHighContainer->begin(),
                       PxlFwdHighContainer->end());

  // cycle through new container
  i = 0; j = 0;
  for (itHit = thePxlFwdHits.begin(); itHit != thePxlFwdHits.end(); ++itHit) {

    ++i;

    // create a DetId from the detUnitId
    DetId theDetUnitId(itHit->detUnitId());
    int detector = theDetUnitId.det();
    int subdetector = theDetUnitId.subdetId();

    // check that expected detector is returned
    if ((detector == dTrk) && (subdetector == sdPxlFwd)) {

      // get the GeomDetUnit from the geometry using theDetUnitID
      const GeomDetUnit *theDet = theTracker.idToDetUnit(theDetUnitId);

      if (!theDet) {
        edm::LogWarning(MsgLoggerCat)
          << "Unable to get GeomDetUnit from PxlFwdHits for Hit " << i;;
        continue;
      }

      ++j;

      // get the Surface of the hit (knows how to go from local <-> global)
      const BoundPlane& bSurface = theDet->surface();

      if(hTrackerPxFToF) hTrackerPxFToF->Fill(itHit->tof());
      if(hTrackerPxFZ) 
        hTrackerPxFZ->Fill(bSurface.toGlobal(itHit->localPosition()).z());
      if(hTrackerPxPhi) 
        hTrackerPxPhi->Fill(bSurface.toGlobal(itHit->localPosition()).phi());
      if(hTrackerPxEta) 
        hTrackerPxEta->Fill(bSurface.toGlobal(itHit->localPosition()).eta());

    } else {
      edm::LogWarning(MsgLoggerCat)
        << "PxlFwd PSimHit " << i 
        << " is expected to be (det,subdet) = (" 
        << dTrk << "," << sdPxlFwd
        << "); value returned is: ("
        << detector << "," << subdetector << ")";
      continue;
    } // end detector type check
  } // end loop through PxlFwd Hits

  if (verbosity > 1) {
    eventout += "\n          Number of Pixel Forward Hits collected:.... ";
    eventout += j;
  }  

  nPxlHits += j;

  if (hTrackerPx[0]) hTrackerPx[0]->Fill((float)nPxlHits);
  if (hTrackerPx[1]) hTrackerPx[1]->Fill((float)nPxlHits); 

  // get Silicon Barrel information
  nSiHits = 0;
  edm::PSimHitContainer theSiBrlHits;
  // extract TIB low container
  edm::Handle<edm::PSimHitContainer> SiTIBLowContainer;
  iEvent.getByLabel(SiTIBLowSrc_,SiTIBLowContainer);
  if (!SiTIBLowContainer.isValid()) {
    edm::LogWarning(MsgLoggerCat)
      << "Unable to find TrackerHitsTIBLowTof in event!";
    return;
  }
  // extract TIB high container
  edm::Handle<edm::PSimHitContainer> SiTIBHighContainer;
  iEvent.getByLabel(SiTIBHighSrc_,SiTIBHighContainer);
  if (!SiTIBHighContainer.isValid()) {
    edm::LogWarning(MsgLoggerCat)
      << "Unable to find TrackerHitsTIBHighTof in event!";
    return;
  }
  // extract TOB low container
  edm::Handle<edm::PSimHitContainer> SiTOBLowContainer;
  iEvent.getByLabel(SiTOBLowSrc_,SiTOBLowContainer);
  if (!SiTOBLowContainer.isValid()) {
    edm::LogWarning(MsgLoggerCat)
      << "Unable to find TrackerHitsTOBLowTof in event!";
    return;
  }
  // extract TOB high container
  edm::Handle<edm::PSimHitContainer> SiTOBHighContainer;
  iEvent.getByLabel(SiTOBHighSrc_,SiTOBHighContainer);
  if (!SiTOBHighContainer.isValid()) {
    edm::LogWarning(MsgLoggerCat)
      << "Unable to find TrackerHitsTOBHighTof in event!";
    return;
  }
  // place all containers into new container
  theSiBrlHits.insert(theSiBrlHits.end(),SiTIBLowContainer->begin(),
                       SiTIBLowContainer->end());
  theSiBrlHits.insert(theSiBrlHits.end(),SiTIBHighContainer->begin(),
                       SiTIBHighContainer->end());
  theSiBrlHits.insert(theSiBrlHits.end(),SiTOBLowContainer->begin(),
                       SiTOBLowContainer->end());
  theSiBrlHits.insert(theSiBrlHits.end(),SiTOBHighContainer->begin(),
                       SiTOBHighContainer->end());

  // cycle through new container
  i = 0; j = 0;
  for (itHit = theSiBrlHits.begin(); itHit != theSiBrlHits.end(); ++itHit) {

    ++i;

    // create a DetId from the detUnitId
    DetId theDetUnitId(itHit->detUnitId());
    int detector = theDetUnitId.det();
    int subdetector = theDetUnitId.subdetId();

    // check that expected detector is returned
    if ((detector == dTrk) && 
        ((subdetector == sdSiTIB) ||
         (subdetector == sdSiTOB))) {

      // get the GeomDetUnit from the geometry using theDetUnitID
      const GeomDetUnit *theDet = theTracker.idToDetUnit(theDetUnitId);

      if (!theDet) {
        edm::LogWarning(MsgLoggerCat)
          << "Unable to get GeomDetUnit from SiBrlHits for Hit " << i;
        continue;
      }

      ++j;

      // get the Surface of the hit (knows how to go from local <-> global)
      const BoundPlane& bSurface = theDet->surface();

      if(hTrackerSiBToF) hTrackerSiBToF->Fill(itHit->tof());
      if(hTrackerSiBR) 
        hTrackerSiBR->Fill(bSurface.toGlobal(itHit->localPosition()).perp());
      if(hTrackerSiPhi) 
        hTrackerSiPhi->Fill(bSurface.toGlobal(itHit->localPosition()).phi());
      if(hTrackerSiEta) 
        hTrackerSiEta->Fill(bSurface.toGlobal(itHit->localPosition()).eta());

    } else {
      edm::LogWarning(MsgLoggerCat)
        << "SiBrl PSimHit " << i 
        << " is expected to be (det,subdet) = (" 
        << dTrk << "," << sdSiTIB
        << " || " << sdSiTOB << "); value returned is: ("
        << detector << "," << subdetector << ")";
      continue;
    } // end detector type check
  } // end loop through SiBrl Hits

  if (verbosity > 1) {
    eventout += "\n          Number of Silicon Barrel Hits collected:... ";
    eventout += j;
  }  

  nSiHits += j;

  // get Silicon Forward information
  edm::PSimHitContainer theSiFwdHits;
  // extract TID low container
  edm::Handle<edm::PSimHitContainer> SiTIDLowContainer;
  iEvent.getByLabel(SiTIDLowSrc_,SiTIDLowContainer);
  if (!SiTIDLowContainer.isValid()) {
    edm::LogWarning(MsgLoggerCat)
      << "Unable to find TrackerHitsTIDLowTof in event!";
    return;
  }
  // extract TID high container
  edm::Handle<edm::PSimHitContainer> SiTIDHighContainer;
  iEvent.getByLabel(SiTIDHighSrc_,SiTIDHighContainer);
  if (!SiTIDHighContainer.isValid()) {
    edm::LogWarning("GlobalHitsProdHist_fillTrk")
      << "Unable to find TrackerHitsTIDHighTof in event!";
    return;
  }
  // extract TEC low container
  edm::Handle<edm::PSimHitContainer> SiTECLowContainer;
  iEvent.getByLabel(SiTECLowSrc_,SiTECLowContainer);
  if (!SiTECLowContainer.isValid()) {
    edm::LogWarning(MsgLoggerCat)
      << "Unable to find TrackerHitsTECLowTof in event!";
    return;
  }
  // extract TEC high container
  edm::Handle<edm::PSimHitContainer> SiTECHighContainer;
  iEvent.getByLabel(SiTECHighSrc_,SiTECHighContainer);
  if (!SiTECHighContainer.isValid()) {
    edm::LogWarning(MsgLoggerCat)
      << "Unable to find TrackerHitsTECHighTof in event!";
    return;
  }
  // place all containers into new container
  theSiFwdHits.insert(theSiFwdHits.end(),SiTIDLowContainer->begin(),
                       SiTIDLowContainer->end());
  theSiFwdHits.insert(theSiFwdHits.end(),SiTIDHighContainer->begin(),
                       SiTIDHighContainer->end());
  theSiFwdHits.insert(theSiFwdHits.end(),SiTECLowContainer->begin(),
                       SiTECLowContainer->end());
  theSiFwdHits.insert(theSiFwdHits.end(),SiTECHighContainer->begin(),
                       SiTECHighContainer->end());

  // cycle through container
  i = 0; j = 0;
  for (itHit = theSiFwdHits.begin(); itHit != theSiFwdHits.end(); ++itHit) {

    ++i;

    // create a DetId from the detUnitId
    DetId theDetUnitId(itHit->detUnitId());
    int detector = theDetUnitId.det();
    int subdetector = theDetUnitId.subdetId();

    // check that expected detector is returned 
    if ((detector == dTrk) && 
        ((subdetector == sdSiTID) ||
         (subdetector == sdSiTEC))) {
      
      // get the GeomDetUnit from the geometry using theDetUnitID
      const GeomDetUnit *theDet = theTracker.idToDetUnit(theDetUnitId);
      
      if (!theDet) {
        edm::LogWarning(MsgLoggerCat)
          << "Unable to get GeomDetUnit from SiFwdHits Hit " << i;
        return;
      }
      
      ++j;

      // get the Surface of the hit (knows how to go from local <-> global)
      const BoundPlane& bSurface = theDet->surface();
      
      if(hTrackerSiFToF) hTrackerSiFToF->Fill(itHit->tof());
      if(hTrackerSiFZ) 
        hTrackerSiFZ->Fill(bSurface.toGlobal(itHit->localPosition()).z());
      if(hTrackerSiPhi) 
        hTrackerSiPhi->Fill(bSurface.toGlobal(itHit->localPosition()).phi());
      if(hTrackerSiEta) 
        hTrackerSiEta->Fill(bSurface.toGlobal(itHit->localPosition()).eta());

    } else {
      edm::LogWarning(MsgLoggerCat)
        << "SiFwd PSimHit " << i 
        << " is expected to be (det,subdet) = (" 
        << dTrk << "," << sdSiTOB
        << " || " << sdSiTEC << "); value returned is: ("
        << detector << "," << subdetector << ")";
      continue;
    } // end check detector type
  } // end loop through SiFwd Hits

  if (verbosity > 1) {
    eventout += "\n          Number of Silicon Forward Hits collected:.. ";
    eventout += j;
  }  

  nSiHits +=j;

  if (hTrackerSi[0]) hTrackerSi[0]->Fill((float)nSiHits);
  if (hTrackerSi[1]) hTrackerSi[1]->Fill((float)nSiHits); 

  if (verbosity > 0)
    edm::LogInfo(MsgLoggerCat) << eventout << "\n";

  return;
}

void GlobalHitsProdHist::fillMuon(edm::Event& iEvent, 
                                  const edm::EventSetup& iSetup)
{
  nMuonHits = 0;
  std::string MsgLoggerCat = "GlobalHitsProdHist_fillMuon";

  TString eventout;
  if (verbosity > 0)
    eventout = "\nGathering info:";  

  // iterator to access containers
  edm::PSimHitContainer::const_iterator itHit;

  // access the CSC Muon
  // access the CSC Muon geometry
  edm::ESHandle<CSCGeometry> theCSCGeometry;
  iSetup.get<MuonGeometryRecord>().get(theCSCGeometry);
  if (!theCSCGeometry.isValid()) {
    edm::LogWarning(MsgLoggerCat)
      << "Unable to find MuonGeometryRecord for the CSCGeometry in event!";
    return;
  }
  const CSCGeometry& theCSCMuon(*theCSCGeometry);

  // get Muon CSC information
  edm::Handle<edm::PSimHitContainer> MuonCSCContainer;
  iEvent.getByLabel(MuonCscSrc_,MuonCSCContainer);
  if (!MuonCSCContainer.isValid()) {
    edm::LogWarning(MsgLoggerCat)
      << "Unable to find MuonCSCHits in event!";
    return;
  }

  // cycle through container
  int i = 0, j = 0;
  for (itHit = MuonCSCContainer->begin(); itHit != MuonCSCContainer->end(); 
       ++itHit) {

    ++i;

    // create a DetId from the detUnitId
    DetId theDetUnitId(itHit->detUnitId());
    int detector = theDetUnitId.det();
    int subdetector = theDetUnitId.subdetId();

    // check that expected detector is returned
    if ((detector == dMuon) && 
        (subdetector == sdMuonCSC)) {

      // get the GeomDetUnit from the geometry using theDetUnitID
      const GeomDetUnit *theDet = theCSCMuon.idToDetUnit(theDetUnitId);
    
      if (!theDet) {
        edm::LogWarning(MsgLoggerCat)
          << "Unable to get GeomDetUnit from theCSCMuon for hit " << i;
        continue;
      }
     
      ++j;

      // get the Surface of the hit (knows how to go from local <-> global)
      const BoundPlane& bSurface = theDet->surface();
    
      if (hMuonCscToF[0]) hMuonCscToF[0]->Fill(itHit->tof());
      if (hMuonCscToF[1]) hMuonCscToF[1]->Fill(itHit->tof());
      if (hMuonCscZ) 
        hMuonCscZ->Fill(bSurface.toGlobal(itHit->localPosition()).z());
      if (hMuonPhi)
        hMuonPhi->Fill(bSurface.toGlobal(itHit->localPosition()).phi());
      if (hMuonEta)
        hMuonEta->Fill(bSurface.toGlobal(itHit->localPosition()).eta());

    } else {
      edm::LogWarning(MsgLoggerCat)
        << "MuonCsc PSimHit " << i 
        << " is expected to be (det,subdet) = (" 
        << dMuon << "," << sdMuonCSC
        << "); value returned is: ("
        << detector << "," << subdetector << ")";
      continue;
    } // end detector type check
  } // end loop through CSC Hits

  if (verbosity > 1) {
    eventout += "\n          Number of CSC muon Hits collected:......... ";
    eventout += j;
  }  

  nMuonHits += j;

  // access the DT Muon
  // access the DT Muon geometry
  edm::ESHandle<DTGeometry> theDTGeometry;
  iSetup.get<MuonGeometryRecord>().get(theDTGeometry);
  if (!theDTGeometry.isValid()) {
    edm::LogWarning(MsgLoggerCat)
      << "Unable to find MuonGeometryRecord for the DTGeometry in event!";
    return;
  }
  const DTGeometry& theDTMuon(*theDTGeometry);

  // get Muon DT information
  edm::Handle<edm::PSimHitContainer> MuonDtContainer;
  iEvent.getByLabel(MuonDtSrc_,MuonDtContainer);
  if (!MuonDtContainer.isValid()) {
    edm::LogWarning(MsgLoggerCat)
      << "Unable to find MuonDTHits in event!";
    return;
  }

  // cycle through container
  i = 0, j = 0;
  for (itHit = MuonDtContainer->begin(); itHit != MuonDtContainer->end(); 
       ++itHit) {

    ++i;

    // create a DetId from the detUnitId
    DetId theDetUnitId(itHit->detUnitId());
    int detector = theDetUnitId.det();
    int subdetector = theDetUnitId.subdetId();

    // check that expected detector is returned
    if ((detector == dMuon) && 
        (subdetector == sdMuonDT)) {

      // CSC uses wires and layers rather than the full detID
      // get the wireId
      DTWireId wireId(itHit->detUnitId());

      // get the DTLayer from the geometry using the wireID
      const DTLayer *theDet = theDTMuon.layer(wireId.layerId());

      if (!theDet) {
        edm::LogWarning(MsgLoggerCat)
          << "Unable to get GeomDetUnit from theDtMuon for hit " << i;
        continue;
      }
     
      ++j;

      // get the Surface of the hit (knows how to go from local <-> global)
      const BoundPlane& bSurface = theDet->surface();
    
      if (hMuonDtToF[0]) hMuonDtToF[0]->Fill(itHit->tof());
      if (hMuonDtToF[1]) hMuonDtToF[1]->Fill(itHit->tof());
      if (hMuonDtR) 
        hMuonDtR->Fill(bSurface.toGlobal(itHit->localPosition()).perp());
      if (hMuonPhi)
        hMuonPhi->Fill(bSurface.toGlobal(itHit->localPosition()).phi());
      if (hMuonEta)
        hMuonEta->Fill(bSurface.toGlobal(itHit->localPosition()).eta());

    } else {
      edm::LogWarning(MsgLoggerCat)
        << "MuonDt PSimHit " << i 
        << " is expected to be (det,subdet) = (" 
        << dMuon << "," << sdMuonDT
        << "); value returned is: ("
        << detector << "," << subdetector << ")";
      continue;
    } // end detector type check
  } // end loop through DT Hits

  if (verbosity > 1) {
    eventout += "\n          Number of DT muon Hits collected:.......... ";
    eventout += j;
  } 

  nMuonHits += j;

  //int RPCBrl = 0, RPCFwd = 0;
  // access the RPC Muon
  // access the RPC Muon geometry
  edm::ESHandle<RPCGeometry> theRPCGeometry;
  iSetup.get<MuonGeometryRecord>().get(theRPCGeometry);
  if (!theRPCGeometry.isValid()) {
    edm::LogWarning(MsgLoggerCat)
      << "Unable to find MuonGeometryRecord for the RPCGeometry in event!";
    return;
  }
  const RPCGeometry& theRPCMuon(*theRPCGeometry);

  // get Muon RPC information
  edm::Handle<edm::PSimHitContainer> MuonRPCContainer;
  iEvent.getByLabel(MuonRpcSrc_,MuonRPCContainer);
  if (!MuonRPCContainer.isValid()) {
    edm::LogWarning(MsgLoggerCat)
      << "Unable to find MuonRPCHits in event!";
    return;
  }

  // cycle through container
  i = 0, j = 0;
  int RPCBrl =0, RPCFwd = 0;
  for (itHit = MuonRPCContainer->begin(); itHit != MuonRPCContainer->end(); 
       ++itHit) {

    ++i;

    // create a DetID from the detUnitId
    DetId theDetUnitId(itHit->detUnitId());
    int detector = theDetUnitId.det();
    int subdetector = theDetUnitId.subdetId();

    // check that expected detector is returned
    if ((detector == dMuon) && 
        (subdetector == sdMuonRPC)) {
      
      // get an RPCDetID from the detUnitID
      RPCDetId RPCId(itHit->detUnitId());      

      // find the region of the RPC hit
      int region = RPCId.region();

      // get the GeomDetUnit from the geometry using the RPCDetId
      const GeomDetUnit *theDet = theRPCMuon.idToDetUnit(theDetUnitId);

      if (!theDet) {
        edm::LogWarning(MsgLoggerCat)
          << "Unable to get GeomDetUnit from theRPCMuon for hit " << i;
        continue;
      }

      ++j;

      // get the Surface of the hit (knows how to go from local <-> global)
      const BoundPlane& bSurface = theDet->surface();
    
      // gather necessary information
      if ((region == sdMuonRPCRgnFwdp) || (region == sdMuonRPCRgnFwdn)) {
        ++RPCFwd;

        if (hMuonRpcFToF[0]) hMuonRpcFToF[0]->Fill(itHit->tof());
        if (hMuonRpcFToF[1]) hMuonRpcFToF[1]->Fill(itHit->tof());
        if (hMuonRpcFZ) 
          hMuonRpcFZ->Fill(bSurface.toGlobal(itHit->localPosition()).z());
        if (hMuonPhi)
          hMuonPhi->Fill(bSurface.toGlobal(itHit->localPosition()).phi());
        if (hMuonEta)
          hMuonEta->Fill(bSurface.toGlobal(itHit->localPosition()).eta());

      } else if (region == sdMuonRPCRgnBrl) {
        ++RPCBrl;

        if (hMuonRpcBToF[0]) hMuonRpcBToF[0]->Fill(itHit->tof());
        if (hMuonRpcBToF[1]) hMuonRpcBToF[1]->Fill(itHit->tof());
        if (hMuonRpcBR) 
          hMuonRpcBR->Fill(bSurface.toGlobal(itHit->localPosition()).perp());
        if (hMuonPhi)
          hMuonPhi->Fill(bSurface.toGlobal(itHit->localPosition()).phi());
        if (hMuonEta)
          hMuonEta->Fill(bSurface.toGlobal(itHit->localPosition()).eta());
        
      } else {
        edm::LogWarning(MsgLoggerCat)
          << "Invalid region for RPC Muon hit" << i;
        continue;
      } // end check of region
    } else {
      edm::LogWarning(MsgLoggerCat)
        << "MuonRpc PSimHit " << i 
        << " is expected to be (det,subdet) = (" 
        << dMuon << "," << sdMuonRPC
        << "); value returned is: ("
        << detector << "," << subdetector << ")";
      continue;
    } // end detector type check
  } // end loop through RPC Hits

  if (verbosity > 1) {
    eventout += "\n          Number of RPC muon Hits collected:......... ";
    eventout += j;
    eventout += "\n                    RPC Barrel muon Hits:............ ";
    eventout += RPCBrl;
    eventout += "\n                    RPC Forward muon Hits:........... ";
    eventout += RPCFwd;    
  }  

  nMuonHits += j;

  if (hMuon[0]) hMuon[0]->Fill((float)nMuonHits);
  if (hMuon[1]) hMuon[1]->Fill((float)nMuonHits); 

  if (verbosity > 0)
    edm::LogInfo(MsgLoggerCat) << eventout << "\n";

  return;
}

void GlobalHitsProdHist::fillECal(edm::Event& iEvent, 
                                  const edm::EventSetup& iSetup)
{
  std::string MsgLoggerCat = "GlobalHitsProdHist_fillECal";

  TString eventout;
  if (verbosity > 0)
    eventout = "\nGathering info:";  
  
  // access the calorimeter geometry
  edm::ESHandle<CaloGeometry> theCaloGeometry;
  iSetup.get<CaloGeometryRecord>().get(theCaloGeometry);
  if (!theCaloGeometry.isValid()) {
    edm::LogWarning(MsgLoggerCat)
      << "Unable to find CaloGeometryRecord in event!";
    return;
  }
  const CaloGeometry& theCalo(*theCaloGeometry);
    
  // iterator to access containers
  edm::PCaloHitContainer::const_iterator itHit;

  // get  ECal information
  edm::PCaloHitContainer theECalHits;
  // extract EB container
  edm::Handle<edm::PCaloHitContainer> EBContainer;
  iEvent.getByLabel(ECalEBSrc_,EBContainer);
  if (!EBContainer.isValid()) {
    edm::LogWarning(MsgLoggerCat)
      << "Unable to find EcalHitsEB in event!";
    return;
  }
  // extract EE container
  edm::Handle<edm::PCaloHitContainer> EEContainer;
  iEvent.getByLabel(ECalEESrc_,EEContainer);
  if (!EEContainer.isValid()) {
    edm::LogWarning(MsgLoggerCat)
      << "Unable to find EcalHitsEE in event!";
    return;
  }
  // place both containers into new container
  theECalHits.insert(theECalHits.end(),EBContainer->begin(),
                       EBContainer->end());
  theECalHits.insert(theECalHits.end(),EEContainer->begin(),
                       EEContainer->end());

  // cycle through new container
  int i = 0, j = 0;
  for (itHit = theECalHits.begin(); itHit != theECalHits.end(); ++itHit) {

    ++i;

    // create a DetId from the detUnitId
    DetId theDetUnitId(itHit->id());
    int detector = theDetUnitId.det();
    int subdetector = theDetUnitId.subdetId();

    // check that expected detector is returned
    if ((detector == dEcal) && 
        ((subdetector == sdEcalBrl) ||
         (subdetector == sdEcalFwd))) {

      // get the Cell geometry
      const CaloCellGeometry *theDet = theCalo.
        getSubdetectorGeometry(theDetUnitId)->getGeometry(theDetUnitId);

      if (!theDet) {
        edm::LogWarning(MsgLoggerCat)
          << "Unable to get CaloCellGeometry from ECalHits for Hit " << i;
        continue;
      }

      ++j;

      // get the global position of the cell
      const GlobalPoint& globalposition = theDet->getPosition();

      if (hCaloEcalE[0]) hCaloEcalE[0]->Fill(itHit->energy());
      if (hCaloEcalE[1]) hCaloEcalE[1]->Fill(itHit->energy());
      if (hCaloEcalToF[0]) hCaloEcalToF[0]->Fill(itHit->time());
      if (hCaloEcalToF[1]) hCaloEcalToF[1]->Fill(itHit->time());
      if (hCaloEcalPhi) hCaloEcalPhi->Fill(globalposition.phi());
      if (hCaloEcalEta) hCaloEcalEta->Fill(globalposition.eta());

    } else {
      edm::LogWarning(MsgLoggerCat)
        << "ECal PCaloHit " << i 
        << " is expected to be (det,subdet) = (" 
        << dEcal << "," << sdEcalBrl
        << " || " << sdEcalFwd << "); value returned is: ("
        << detector << "," << subdetector << ")";
      continue;
    } // end detector type check
  } // end loop through ECal Hits

  if (verbosity > 1) {
    eventout += "\n          Number of ECal Hits collected:............. ";
    eventout += j;
  }  

  if (hCaloEcal[0]) hCaloEcal[0]->Fill((float)j);
  if (hCaloEcal[1]) hCaloEcal[1]->Fill((float)j); 

  // Get Preshower information
  // extract PreShower container
  edm::Handle<edm::PCaloHitContainer> PreShContainer;
  iEvent.getByLabel(ECalESSrc_,PreShContainer);
  if (!PreShContainer.isValid()) {
    edm::LogWarning(MsgLoggerCat)
      << "Unable to find EcalHitsES in event!";
    return;
  }

  // cycle through container
  i = 0, j = 0;
  for (itHit = PreShContainer->begin(); 
       itHit != PreShContainer->end(); ++itHit) {

    ++i;

    // create a DetId from the detUnitId
    DetId theDetUnitId(itHit->id());
    int detector = theDetUnitId.det();
    int subdetector = theDetUnitId.subdetId();

    // check that expected detector is returned
    if ((detector == dEcal) && 
        (subdetector == sdEcalPS)) {

      // get the Cell geometry
      const CaloCellGeometry *theDet = theCalo.
        getSubdetectorGeometry(theDetUnitId)->getGeometry(theDetUnitId);

      if (!theDet) {
        edm::LogWarning(MsgLoggerCat)
          << "Unable to get CaloCellGeometry from PreShContainer for Hit " 
          << i;
        continue;
      }

      ++j;

      // get the global position of the cell
      const GlobalPoint& globalposition = theDet->getPosition();

      if (hCaloPreShE[0]) hCaloPreShE[0]->Fill(itHit->energy());
      if (hCaloPreShE[1]) hCaloPreShE[1]->Fill(itHit->energy());
      if (hCaloPreShToF[0]) hCaloPreShToF[0]->Fill(itHit->time());
      if (hCaloPreShToF[1]) hCaloPreShToF[1]->Fill(itHit->time());
      if (hCaloPreShPhi) hCaloPreShPhi->Fill(globalposition.phi());
      if (hCaloPreShEta) hCaloPreShEta->Fill(globalposition.eta());

    } else {
      edm::LogWarning(MsgLoggerCat)
        << "PreSh PCaloHit " << i 
        << " is expected to be (det,subdet) = (" 
        << dEcal << "," << sdEcalPS
        << "); value returned is: ("
        << detector << "," << subdetector << ")";
      continue;
    } // end detector type check
  } // end loop through PreShower Hits

  if (verbosity > 1) {
    eventout += "\n          Number of PreSh Hits collected:............ ";
    eventout += j;
  }  

  if (hCaloPreSh[0]) hCaloPreSh[0]->Fill((float)j);
  if (hCaloPreSh[1]) hCaloPreSh[1]->Fill((float)j); 

  if (verbosity > 0)
    edm::LogInfo(MsgLoggerCat) << eventout << "\n";

  return;
}

void GlobalHitsProdHist::fillHCal(edm::Event& iEvent, 
                                  const edm::EventSetup& iSetup)
{
  std::string MsgLoggerCat = "GlobalHitsProdHist_fillHCal";

  TString eventout;
  if (verbosity > 0)
    eventout = "\nGathering info:";  
  
  // access the calorimeter geometry
  edm::ESHandle<CaloGeometry> theCaloGeometry;
  iSetup.get<CaloGeometryRecord>().get(theCaloGeometry);
  if (!theCaloGeometry.isValid()) {
    edm::LogWarning(MsgLoggerCat)
      << "Unable to find CaloGeometryRecord in event!";
    return;
  }
  const CaloGeometry& theCalo(*theCaloGeometry);
    
  // iterator to access containers
  edm::PCaloHitContainer::const_iterator itHit;

  // get  HCal information
  // extract HCal container
  edm::Handle<edm::PCaloHitContainer> HCalContainer;
  iEvent.getByLabel(HCalSrc_,HCalContainer);
  if (!HCalContainer.isValid()) {
    edm::LogWarning(MsgLoggerCat)
      << "Unable to find HCalHits in event!";
    return;
  }

  // cycle through container
  int i = 0, j = 0;
  for (itHit = HCalContainer->begin(); 
       itHit != HCalContainer->end(); ++itHit) {

    ++i;

    // create a DetId from the detUnitId
    DetId theDetUnitId(itHit->id());
    int detector = theDetUnitId.det();
    int subdetector = theDetUnitId.subdetId();

    // check that expected detector is returned
    if ((detector == dHcal) && 
        ((subdetector == sdHcalBrl) ||
         (subdetector == sdHcalEC) ||
         (subdetector == sdHcalOut) ||
         (subdetector == sdHcalFwd))) {

      // get the Cell geometry
      const CaloCellGeometry *theDet = theCalo.
        getSubdetectorGeometry(theDetUnitId)->getGeometry(theDetUnitId);

      if (!theDet) {
        edm::LogWarning(MsgLoggerCat)
          << "Unable to get CaloCellGeometry from HCalContainer for Hit " << i;
        continue;
      }

      ++j;

      // get the global position of the cell
      const GlobalPoint& globalposition = theDet->getPosition();

      if (hCaloHcalE[0]) hCaloHcalE[0]->Fill(itHit->energy());
      if (hCaloHcalE[1]) hCaloHcalE[1]->Fill(itHit->energy());
      if (hCaloHcalToF[0]) hCaloHcalToF[0]->Fill(itHit->time());
      if (hCaloHcalToF[1]) hCaloHcalToF[1]->Fill(itHit->time());
      if (hCaloHcalPhi) hCaloHcalPhi->Fill(globalposition.phi());
      if (hCaloHcalEta) hCaloHcalEta->Fill(globalposition.eta());

    } else {
      edm::LogWarning(MsgLoggerCat)
        << "HCal PCaloHit " << i 
        << " is expected to be (det,subdet) = (" 
        << dHcal << "," << sdHcalBrl
        << " || " << sdHcalEC << " || " << sdHcalOut << " || " << sdHcalFwd
        << "); value returned is: ("
        << detector << "," << subdetector << ")";
      continue;
    } // end detector type check
  } // end loop through HCal Hits

  if (verbosity > 1) {
    eventout += "\n          Number of HCal Hits collected:............. ";
    eventout += j;
  }  

  if (hCaloHcal[0]) hCaloHcal[0]->Fill((float)j);
  if (hCaloHcal[1]) hCaloHcal[1]->Fill((float)j); 

  if (verbosity > 0)
    edm::LogInfo(MsgLoggerCat) << eventout << "\n";

  return;
}