Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:32:22

0001 /** \file GlobalHitsAnalyzer.cc
0002  *
0003  *  See header file for description of class
0004  *
0005  *  \author M. Strang SUNY-Buffalo
0006  */
0007 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0008 
0009 #include "DQMServices/Core/interface/DQMStore.h"
0010 #include "Validation/GlobalHits/interface/GlobalHitsAnalyzer.h"
0011 
0012 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0013 #include "Geometry/HcalCommonData/interface/HcalDDDRecConstants.h"
0014 #include "Geometry/HcalCommonData/interface/HcalHitRelabeller.h"
0015 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0016 #include "Geometry/Records/interface/HcalRecNumberingRecord.h"
0017 #include "FWCore/Framework/interface/GetterOfProducts.h"
0018 #include "FWCore/Framework/interface/ProcessMatch.h"
0019 
0020 GlobalHitsAnalyzer::GlobalHitsAnalyzer(const edm::ParameterSet &iPSet)
0021     : fName(""),
0022       verbosity(0),
0023       frequency(0),
0024       vtxunit(0),
0025       label(""),
0026       getAllProvenances(false),
0027       printProvenanceInfo(false),
0028       testNumber(false),
0029       G4VtxSrc_Token_(consumes<edm::SimVertexContainer>((iPSet.getParameter<edm::InputTag>("G4VtxSrc")))),
0030       G4TrkSrc_Token_(consumes<edm::SimTrackContainer>(iPSet.getParameter<edm::InputTag>("G4TrkSrc"))),
0031       tGeomToken_(esConsumes()),
0032       cscGeomToken_(esConsumes()),
0033       dtGeomToken_(esConsumes()),
0034       rpcGeomToken_(esConsumes()),
0035       caloGeomToken_(esConsumes()),
0036       hcaldddRecToken_(esConsumes()),
0037       count(0) {
0038   getterOfProducts_ = edm::GetterOfProducts<edm::HepMCProduct>(edm::ProcessMatch("*"), this);
0039   callWhenNewProductsRegistered(getterOfProducts_);
0040 
0041   std::string MsgLoggerCat = "GlobalHitsAnalyzer_GlobalHitsAnalyzer";
0042 
0043   // get information from parameter set
0044   fName = iPSet.getUntrackedParameter<std::string>("Name");
0045   verbosity = iPSet.getUntrackedParameter<int>("Verbosity");
0046   frequency = iPSet.getUntrackedParameter<int>("Frequency");
0047   vtxunit = iPSet.getUntrackedParameter<int>("VtxUnit");
0048   edm::ParameterSet m_Prov = iPSet.getParameter<edm::ParameterSet>("ProvenanceLookup");
0049   getAllProvenances = m_Prov.getUntrackedParameter<bool>("GetAllProvenances");
0050   printProvenanceInfo = m_Prov.getUntrackedParameter<bool>("PrintProvenanceInfo");
0051   testNumber = iPSet.getUntrackedParameter<bool>("testNumber");
0052 
0053   // get Labels to use to extract information
0054   PxlBrlLowSrc_ = iPSet.getParameter<edm::InputTag>("PxlBrlLowSrc");
0055   PxlBrlHighSrc_ = iPSet.getParameter<edm::InputTag>("PxlBrlHighSrc");
0056   PxlFwdLowSrc_ = iPSet.getParameter<edm::InputTag>("PxlFwdLowSrc");
0057   PxlFwdHighSrc_ = iPSet.getParameter<edm::InputTag>("PxlFwdHighSrc");
0058 
0059   SiTIBLowSrc_ = iPSet.getParameter<edm::InputTag>("SiTIBLowSrc");
0060   SiTIBHighSrc_ = iPSet.getParameter<edm::InputTag>("SiTIBHighSrc");
0061   SiTOBLowSrc_ = iPSet.getParameter<edm::InputTag>("SiTOBLowSrc");
0062   SiTOBHighSrc_ = iPSet.getParameter<edm::InputTag>("SiTOBHighSrc");
0063   SiTIDLowSrc_ = iPSet.getParameter<edm::InputTag>("SiTIDLowSrc");
0064   SiTIDHighSrc_ = iPSet.getParameter<edm::InputTag>("SiTIDHighSrc");
0065   SiTECLowSrc_ = iPSet.getParameter<edm::InputTag>("SiTECLowSrc");
0066   SiTECHighSrc_ = iPSet.getParameter<edm::InputTag>("SiTECHighSrc");
0067 
0068   MuonCscSrc_ = iPSet.getParameter<edm::InputTag>("MuonCscSrc");
0069   MuonDtSrc_ = iPSet.getParameter<edm::InputTag>("MuonDtSrc");
0070   MuonRpcSrc_ = iPSet.getParameter<edm::InputTag>("MuonRpcSrc");
0071 
0072   ECalEBSrc_ = iPSet.getParameter<edm::InputTag>("ECalEBSrc");
0073   ECalEESrc_ = iPSet.getParameter<edm::InputTag>("ECalEESrc");
0074   ECalESSrc_ = iPSet.getParameter<edm::InputTag>("ECalESSrc");
0075 
0076   HCalSrc_ = iPSet.getParameter<edm::InputTag>("HCalSrc");
0077 
0078   // fix for consumes
0079   PxlBrlLowSrc_Token_ = consumes<edm::PSimHitContainer>(iPSet.getParameter<edm::InputTag>("PxlBrlLowSrc"));
0080   PxlBrlHighSrc_Token_ = consumes<edm::PSimHitContainer>(iPSet.getParameter<edm::InputTag>("PxlBrlHighSrc"));
0081   PxlFwdLowSrc_Token_ = consumes<edm::PSimHitContainer>(iPSet.getParameter<edm::InputTag>("PxlFwdLowSrc"));
0082   PxlFwdHighSrc_Token_ = consumes<edm::PSimHitContainer>(iPSet.getParameter<edm::InputTag>("PxlFwdHighSrc"));
0083 
0084   SiTIBLowSrc_Token_ = consumes<edm::PSimHitContainer>(iPSet.getParameter<edm::InputTag>("SiTIBLowSrc"));
0085   SiTIBHighSrc_Token_ = consumes<edm::PSimHitContainer>(iPSet.getParameter<edm::InputTag>("SiTIBHighSrc"));
0086   SiTOBLowSrc_Token_ = consumes<edm::PSimHitContainer>(iPSet.getParameter<edm::InputTag>("SiTOBLowSrc"));
0087   SiTOBHighSrc_Token_ = consumes<edm::PSimHitContainer>(iPSet.getParameter<edm::InputTag>("SiTOBHighSrc"));
0088   SiTIDLowSrc_Token_ = consumes<edm::PSimHitContainer>(iPSet.getParameter<edm::InputTag>("SiTIDLowSrc"));
0089   SiTIDHighSrc_Token_ = consumes<edm::PSimHitContainer>(iPSet.getParameter<edm::InputTag>("SiTIDHighSrc"));
0090   SiTECLowSrc_Token_ = consumes<edm::PSimHitContainer>(iPSet.getParameter<edm::InputTag>("SiTECLowSrc"));
0091   SiTECHighSrc_Token_ = consumes<edm::PSimHitContainer>(iPSet.getParameter<edm::InputTag>("SiTECHighSrc"));
0092 
0093   MuonCscSrc_Token_ = consumes<edm::PSimHitContainer>(iPSet.getParameter<edm::InputTag>("MuonCscSrc"));
0094   MuonDtSrc_Token_ = consumes<edm::PSimHitContainer>(iPSet.getParameter<edm::InputTag>("MuonDtSrc"));
0095   MuonRpcSrc_Token_ = consumes<edm::PSimHitContainer>(iPSet.getParameter<edm::InputTag>("MuonRpcSrc"));
0096 
0097   ECalEBSrc_Token_ = consumes<edm::PCaloHitContainer>(iPSet.getParameter<edm::InputTag>("ECalEBSrc"));
0098   ECalEESrc_Token_ = consumes<edm::PCaloHitContainer>(iPSet.getParameter<edm::InputTag>("ECalEESrc"));
0099   ECalESSrc_Token_ = consumes<edm::PCaloHitContainer>(iPSet.getParameter<edm::InputTag>("ECalESSrc"));
0100   HCalSrc_Token_ = consumes<edm::PCaloHitContainer>(iPSet.getParameter<edm::InputTag>("HCalSrc"));
0101 
0102   // determine whether to process subdetector or not
0103   validHepMCevt = iPSet.getUntrackedParameter<bool>("validHepMCevt");
0104   validG4VtxContainer = iPSet.getUntrackedParameter<bool>("validG4VtxContainer");
0105   validG4trkContainer = iPSet.getUntrackedParameter<bool>("validG4trkContainer");
0106   validPxlBrlLow = iPSet.getUntrackedParameter<bool>("validPxlBrlLow", true);
0107   validPxlBrlHigh = iPSet.getUntrackedParameter<bool>("validPxlBrlHigh", true);
0108   validPxlFwdLow = iPSet.getUntrackedParameter<bool>("validPxlFwdLow", true);
0109   validPxlFwdHigh = iPSet.getUntrackedParameter<bool>("validPxlFwdHigh", true);
0110   validSiTIBLow = iPSet.getUntrackedParameter<bool>("validSiTIBLow", true);
0111   validSiTIBHigh = iPSet.getUntrackedParameter<bool>("validSiTIBHigh", true);
0112   validSiTOBLow = iPSet.getUntrackedParameter<bool>("validSiTOBLow", true);
0113   validSiTOBHigh = iPSet.getUntrackedParameter<bool>("validSiTOBHigh", true);
0114   validSiTIDLow = iPSet.getUntrackedParameter<bool>("validSiTIDLow", true);
0115   validSiTIDHigh = iPSet.getUntrackedParameter<bool>("validSiTIDHigh", true);
0116   validSiTECLow = iPSet.getUntrackedParameter<bool>("validSiTECLow", true);
0117   validSiTECHigh = iPSet.getUntrackedParameter<bool>("validSiTECHigh", true);
0118   validMuonCSC = iPSet.getUntrackedParameter<bool>("validMuonCSC", true);
0119   validMuonDt = iPSet.getUntrackedParameter<bool>("validMuonDt", true);
0120   validMuonRPC = iPSet.getUntrackedParameter<bool>("validMuonRPC", true);
0121   validEB = iPSet.getUntrackedParameter<bool>("validEB", true);
0122   validEE = iPSet.getUntrackedParameter<bool>("validEE", true);
0123   validPresh = iPSet.getUntrackedParameter<bool>("validPresh", true);
0124   validHcal = iPSet.getUntrackedParameter<bool>("validHcal", true);
0125 
0126   // use value of first digit to determine default output level (inclusive)
0127   // 0 is none, 1 is basic, 2 is fill output, 3 is gather output
0128   verbosity %= 10;
0129 
0130   // print out Parameter Set information being used
0131   if (verbosity >= 0) {
0132     edm::LogInfo(MsgLoggerCat)
0133         << "\n===============================\n"
0134         << "Initialized as EDAnalyzer with parameter values:\n"
0135         << "    Name                  = " << fName << "\n"
0136         << "    Verbosity             = " << verbosity << "\n"
0137         << "    Frequency             = " << frequency << "\n"
0138         << "    VtxUnit               = " << vtxunit << "\n"
0139         << "    GetProv               = " << getAllProvenances << "\n"
0140         << "    PrintProv             = " << printProvenanceInfo << "\n"
0141         << "    PxlBrlLowSrc          = " << PxlBrlLowSrc_.label() << ":" << PxlBrlLowSrc_.instance() << "\n"
0142         << "    PxlBrlHighSrc         = " << PxlBrlHighSrc_.label() << ":" << PxlBrlHighSrc_.instance() << "\n"
0143         << "    PxlFwdLowSrc          = " << PxlFwdLowSrc_.label() << ":" << PxlBrlLowSrc_.instance() << "\n"
0144         << "    PxlFwdHighSrc         = " << PxlFwdHighSrc_.label() << ":" << PxlBrlHighSrc_.instance() << "\n"
0145         << "    SiTIBLowSrc           = " << SiTIBLowSrc_.label() << ":" << SiTIBLowSrc_.instance() << "\n"
0146         << "    SiTIBHighSrc          = " << SiTIBHighSrc_.label() << ":" << SiTIBHighSrc_.instance() << "\n"
0147         << "    SiTOBLowSrc           = " << SiTOBLowSrc_.label() << ":" << SiTOBLowSrc_.instance() << "\n"
0148         << "    SiTOBHighSrc          = " << SiTOBHighSrc_.label() << ":" << SiTOBHighSrc_.instance() << "\n"
0149         << "    SiTIDLowSrc           = " << SiTIDLowSrc_.label() << ":" << SiTIDLowSrc_.instance() << "\n"
0150         << "    SiTIDHighSrc          = " << SiTIDHighSrc_.label() << ":" << SiTIDHighSrc_.instance() << "\n"
0151         << "    SiTECLowSrc           = " << SiTECLowSrc_.label() << ":" << SiTECLowSrc_.instance() << "\n"
0152         << "    SiTECHighSrc          = " << SiTECHighSrc_.label() << ":" << SiTECHighSrc_.instance() << "\n"
0153         << "    MuonCscSrc            = " << MuonCscSrc_.label() << ":" << MuonCscSrc_.instance() << "\n"
0154         << "    MuonDtSrc             = " << MuonDtSrc_.label() << ":" << MuonDtSrc_.instance() << "\n"
0155         << "    MuonRpcSrc            = " << MuonRpcSrc_.label() << ":" << MuonRpcSrc_.instance() << "\n"
0156         << "    ECalEBSrc             = " << ECalEBSrc_.label() << ":" << ECalEBSrc_.instance() << "\n"
0157         << "    ECalEESrc             = " << ECalEESrc_.label() << ":" << ECalEESrc_.instance() << "\n"
0158         << "    ECalESSrc             = " << ECalESSrc_.label() << ":" << ECalESSrc_.instance() << "\n"
0159         << "    HCalSrc               = " << HCalSrc_.label() << ":" << HCalSrc_.instance() << "\n"
0160         << "\n"
0161         << "    validHepMCevt         = "
0162         << ":" << validHepMCevt << "\n"
0163         << "    validG4VtxContainer   = "
0164         << ":" << validG4VtxContainer << "\n"
0165         << "    validG4trkContainer   = "
0166         << ":" << validG4trkContainer << "\n"
0167         << "    validPxlBrlLow        = "
0168         << ":" << validPxlBrlLow << "\n"
0169         << "    validPxlBrlHigh       = "
0170         << ":" << validPxlBrlHigh << "\n"
0171         << "    validPxlFwdLow        = "
0172         << ":" << validPxlFwdLow << "\n"
0173         << "    validPxlFwdHigh       = "
0174         << ":" << validPxlFwdHigh << "\n"
0175         << "    validSiTIBLow         = "
0176         << ":" << validSiTIBLow << "\n"
0177         << "    validSiTIBHigh        = "
0178         << ":" << validSiTIBHigh << "\n"
0179         << "    validSiTOBLow         = "
0180         << ":" << validSiTOBLow << "\n"
0181         << "    validSiTOBHigh        = "
0182         << ":" << validSiTOBHigh << "\n"
0183         << "    validSiTIDLow         = "
0184         << ":" << validSiTIDLow << "\n"
0185         << "    validSiTIDHigh        = "
0186         << ":" << validSiTIDHigh << "\n"
0187         << "    validSiTECLow         = "
0188         << ":" << validSiTECLow << "\n"
0189         << "    validSiTECHigh        = "
0190         << ":" << validSiTECHigh << "\n"
0191         << "    validMuonCSC          = "
0192         << ":" << validMuonCSC << "\n"
0193         << "    validMuonDt           = "
0194         << ":" << validMuonDt << "\n"
0195         << "    validMuonRPC          = "
0196         << ":" << validMuonRPC << "\n"
0197         << "    validEB               = "
0198         << ":" << validEB << "\n"
0199         << "    validEE               = "
0200         << ":" << validEE << "\n"
0201         << "    validPresh            = "
0202         << ":" << validPresh << "\n"
0203         << "    validHcal             = "
0204         << ":" << validHcal << "\n"
0205         << "===============================\n";
0206   }
0207 
0208   // initialize monitor elements
0209   for (Int_t i = 0; i < 2; ++i) {
0210     meMCRGP[i] = nullptr;
0211     meMCG4Vtx[i] = nullptr;
0212     meGeantVtxX[i] = nullptr;
0213     meGeantVtxY[i] = nullptr;
0214     meGeantVtxZ[i] = nullptr;
0215     meMCG4Trk[i] = nullptr;
0216     meCaloEcal[i] = nullptr;
0217     meCaloEcalE[i] = nullptr;
0218     meCaloEcalToF[i] = nullptr;
0219     meCaloPreSh[i] = nullptr;
0220     meCaloPreShE[i] = nullptr;
0221     meCaloPreShToF[i] = nullptr;
0222     meCaloHcal[i] = nullptr;
0223     meCaloHcalE[i] = nullptr;
0224     meCaloHcalToF[i] = nullptr;
0225     meTrackerPx[i] = nullptr;
0226     meTrackerSi[i] = nullptr;
0227     meMuon[i] = nullptr;
0228     meMuonDtToF[i] = nullptr;
0229     meMuonCscToF[i] = nullptr;
0230     meMuonRpcFToF[i] = nullptr;
0231     meMuonRpcBToF[i] = nullptr;
0232     meGeantVtxRad[i] = nullptr;
0233   }
0234   meGeantTrkPt = nullptr;
0235   meGeantTrkE = nullptr;
0236   meGeantVtxEta = nullptr;
0237   meGeantVtxPhi = nullptr;
0238   meGeantVtxMulti = nullptr;
0239   meCaloEcalPhi = nullptr;
0240   meCaloEcalEta = nullptr;
0241   meCaloPreShPhi = nullptr;
0242   meCaloPreShEta = nullptr;
0243   meCaloHcalPhi = nullptr;
0244   meCaloHcalEta = nullptr;
0245   meTrackerPxPhi = nullptr;
0246   meTrackerPxEta = nullptr;
0247   meTrackerPxBToF = nullptr;
0248   meTrackerPxBR = nullptr;
0249   meTrackerPxFToF = nullptr;
0250   meTrackerPxFZ = nullptr;
0251   meTrackerSiPhi = nullptr;
0252   meTrackerSiEta = nullptr;
0253   meTrackerSiBToF = nullptr;
0254   meTrackerSiBR = nullptr;
0255   meTrackerSiFToF = nullptr;
0256   meTrackerSiFZ = nullptr;
0257   meMuonPhi = nullptr;
0258   meMuonEta = nullptr;
0259   meMuonDtR = nullptr;
0260   meMuonCscZ = nullptr;
0261   meMuonRpcBR = nullptr;
0262   meMuonRpcFZ = nullptr;
0263 }
0264 
0265 GlobalHitsAnalyzer::~GlobalHitsAnalyzer() {}
0266 
0267 void GlobalHitsAnalyzer::bookHistograms(DQMStore::IBooker &iBooker, edm::Run const &run, edm::EventSetup const &es) {
0268   // book histograms
0269   Char_t hname[200];
0270   Char_t htitle[200];
0271 
0272   // MCGeant
0273   iBooker.setCurrentFolder("GlobalHitsV/MCGeant");
0274   sprintf(hname, "hMCRGP1");
0275   sprintf(htitle, "RawGenParticles");
0276   meMCRGP[0] = iBooker.book1D(hname, htitle, 100, 0., 5000.);
0277   sprintf(hname, "hMCRGP2");
0278   meMCRGP[1] = iBooker.book1D(hname, htitle, 100, 0., 500.);
0279   for (Int_t i = 0; i < 2; ++i) {
0280     meMCRGP[i]->setAxisTitle("Number of Raw Generated Particles", 1);
0281     meMCRGP[i]->setAxisTitle("Count", 2);
0282   }
0283 
0284   sprintf(hname, "hMCG4Vtx1");
0285   sprintf(htitle, "G4 Vertices");
0286   meMCG4Vtx[0] = iBooker.book1D(hname, htitle, 150, 0., 15000.);
0287   sprintf(hname, "hMCG4Vtx2");
0288   meMCG4Vtx[1] = iBooker.book1D(hname, htitle, 100, -0.5, 99.5);
0289   for (Int_t i = 0; i < 2; ++i) {
0290     meMCG4Vtx[i]->setAxisTitle("Number of Vertices", 1);
0291     meMCG4Vtx[i]->setAxisTitle("Count", 2);
0292   }
0293 
0294   sprintf(hname, "hMCG4Trk1");
0295   sprintf(htitle, "G4 Tracks");
0296   meMCG4Trk[0] = iBooker.book1D(hname, htitle, 150, 0., 15000.);
0297   sprintf(hname, "hMCG4Trk2");
0298   meMCG4Trk[1] = iBooker.book1D(hname, htitle, 150, -0.5, 99.5);
0299   for (Int_t i = 0; i < 2; ++i) {
0300     meMCG4Trk[i]->setAxisTitle("Number of Tracks", 1);
0301     meMCG4Trk[i]->setAxisTitle("Count", 2);
0302   }
0303 
0304   sprintf(hname, "hGeantVtxX1");
0305   sprintf(htitle, "Geant vertex x/micrometer");
0306   meGeantVtxX[0] = iBooker.book1D(hname, htitle, 100, -8000000., 8000000.);
0307   sprintf(hname, "hGeantVtxX2");
0308   meGeantVtxX[1] = iBooker.book1D(hname, htitle, 100, -50., 50.);
0309   for (Int_t i = 0; i < 2; ++i) {
0310     meGeantVtxX[i]->setAxisTitle("x of Vertex (um)", 1);
0311     meGeantVtxX[i]->setAxisTitle("Count", 2);
0312   }
0313 
0314   sprintf(hname, "hGeantVtxY1");
0315   sprintf(htitle, "Geant vertex y/micrometer");
0316   meGeantVtxY[0] = iBooker.book1D(hname, htitle, 100, -8000000, 8000000.);
0317   sprintf(hname, "hGeantVtxY2");
0318   meGeantVtxY[1] = iBooker.book1D(hname, htitle, 100, -50., 50.);
0319   for (Int_t i = 0; i < 2; ++i) {
0320     meGeantVtxY[i]->setAxisTitle("y of Vertex (um)", 1);
0321     meGeantVtxY[i]->setAxisTitle("Count", 2);
0322   }
0323 
0324   sprintf(hname, "hGeantVtxZ1");
0325   sprintf(htitle, "Geant vertex z/millimeter");
0326   meGeantVtxZ[0] = iBooker.book1D(hname, htitle, 100, -11000., 11000.);
0327   sprintf(hname, "hGeantVtxZ2");
0328   meGeantVtxZ[1] = iBooker.book1D(hname, htitle, 200, -500., 500.);
0329   // meGeantVtxZ[1] = iBooker.book1D(hname,htitle,100,-250.,250.);
0330   for (Int_t i = 0; i < 2; ++i) {
0331     meGeantVtxZ[i]->setAxisTitle("z of Vertex (mm)", 1);
0332     meGeantVtxZ[i]->setAxisTitle("Count", 2);
0333   }
0334 
0335   sprintf(hname, "hGeantTrkPt");
0336   sprintf(htitle, "Log10 Geant track pt/GeV");
0337   meGeantTrkPt = iBooker.book1D(hname, htitle, 80, -4., 4.);
0338   meGeantTrkPt->setAxisTitle("Log10 pT of Track (GeV)", 1);
0339   meGeantTrkPt->setAxisTitle("Count", 2);
0340 
0341   sprintf(hname, "hGeantTrkE");
0342   sprintf(htitle, "Log10 Geant track E/GeV");
0343   meGeantTrkE = iBooker.book1D(hname, htitle, 80, -4., 4.);
0344   meGeantTrkE->setAxisTitle("Log10 E of Track (GeV)", 1);
0345   meGeantTrkE->setAxisTitle("Count", 2);
0346 
0347   sprintf(hname, "hGeantVtxEta");
0348   sprintf(htitle, "Geant vertices eta");
0349   meGeantVtxEta = iBooker.book1D(hname, htitle, 220, -5.5, 5.5);
0350   meGeantVtxEta->setAxisTitle("eta of SimVertex", 1);
0351   meGeantVtxEta->setAxisTitle("Count", 2);
0352 
0353   sprintf(hname, "hGeantVtxPhi");
0354   sprintf(htitle, "Geant vertices phi/rad");
0355   meGeantVtxPhi = iBooker.book1D(hname, htitle, 100, -3.2, 3.2);
0356   meGeantVtxPhi->setAxisTitle("phi of SimVertex (rad)", 1);
0357   meGeantVtxPhi->setAxisTitle("Count", 2);
0358 
0359   sprintf(hname, "hGeantVtxRad1");
0360   sprintf(htitle, "Geant vertices radius/cm");
0361   meGeantVtxRad[0] = iBooker.book1D(hname, htitle, 130, 0., 130.);
0362   sprintf(hname, "hGeantVtxRad2");
0363   meGeantVtxRad[1] = iBooker.book1D(hname, htitle, 100, 0., 1000.);
0364   for (Int_t i = 0; i < 2; ++i) {
0365     meGeantVtxRad[i]->setAxisTitle("radius of SimVertex (cm)", 1);
0366     meGeantVtxRad[i]->setAxisTitle("Count", 2);
0367   }
0368 
0369   sprintf(hname, "hGeantVtxMulti");
0370   sprintf(htitle, "Geant vertices outgoing multiplicity");
0371   meGeantVtxMulti = iBooker.book1D(hname, htitle, 20, 0., 20);
0372   meGeantVtxMulti->setAxisTitle("multiplicity of particles attached to a SimVertex", 1);
0373   meGeantVtxMulti->setAxisTitle("Count", 2);
0374 
0375   // ECal
0376   iBooker.setCurrentFolder("GlobalHitsV/ECals");
0377   sprintf(hname, "hCaloEcal1");
0378   sprintf(htitle, "Ecal hits");
0379   meCaloEcal[0] = iBooker.book1D(hname, htitle, 100, 0., 10000.);
0380   sprintf(hname, "hCaloEcal2");
0381   meCaloEcal[1] = iBooker.book1D(hname, htitle, 100, -0.5, 99.5);
0382 
0383   sprintf(hname, "hCaloEcalE1");
0384   sprintf(htitle, "Ecal hits, energy/GeV");
0385   meCaloEcalE[0] = iBooker.book1D(hname, htitle, 100, 0., 10.);
0386   sprintf(hname, "hCaloEcalE2");
0387   meCaloEcalE[1] = iBooker.book1D(hname, htitle, 100, 0., 0.1);
0388   sprintf(hname, "hCaloEcalToF1");
0389   sprintf(htitle, "Ecal hits, ToF/ns");
0390   meCaloEcalToF[0] = iBooker.book1D(hname, htitle, 100, 0., 1000.);
0391   sprintf(hname, "hCaloEcalToF2");
0392   meCaloEcalToF[1] = iBooker.book1D(hname, htitle, 100, 0., 100.);
0393 
0394   for (Int_t i = 0; i < 2; ++i) {
0395     meCaloEcal[i]->setAxisTitle("Number of Hits", 1);
0396     meCaloEcal[i]->setAxisTitle("Count", 2);
0397     meCaloEcalE[i]->setAxisTitle("Energy of Hits (GeV)", 1);
0398     meCaloEcalE[i]->setAxisTitle("Count", 2);
0399     meCaloEcalToF[i]->setAxisTitle("Time of Flight of Hits (ns)", 1);
0400     meCaloEcalToF[i]->setAxisTitle("Count", 2);
0401   }
0402 
0403   sprintf(hname, "hCaloEcalPhi");
0404   sprintf(htitle, "Ecal hits, phi/rad");
0405   meCaloEcalPhi = iBooker.book1D(hname, htitle, 100, -3.2, 3.2);
0406   meCaloEcalPhi->setAxisTitle("Phi of Hits (rad)", 1);
0407   meCaloEcalPhi->setAxisTitle("Count", 2);
0408 
0409   sprintf(hname, "hCaloEcalEta");
0410   sprintf(htitle, "Ecal hits, eta");
0411   meCaloEcalEta = iBooker.book1D(hname, htitle, 100, -5.5, 5.5);
0412   meCaloEcalEta->setAxisTitle("Eta of Hits", 1);
0413   meCaloEcalEta->setAxisTitle("Count", 2);
0414 
0415   sprintf(hname, "hCaloPreSh1");
0416   sprintf(htitle, "PreSh hits");
0417   meCaloPreSh[0] = iBooker.book1D(hname, htitle, 100, 0., 10000.);
0418   sprintf(hname, "hCaloPreSh2");
0419   meCaloPreSh[1] = iBooker.book1D(hname, htitle, 100, -0.5, 99.5);
0420 
0421   sprintf(hname, "hCaloPreShE1");
0422   sprintf(htitle, "PreSh hits, energy/GeV");
0423   meCaloPreShE[0] = iBooker.book1D(hname, htitle, 100, 0., 10.);
0424   sprintf(hname, "hCaloPreShE2");
0425   meCaloPreShE[1] = iBooker.book1D(hname, htitle, 100, 0., 0.1);
0426 
0427   sprintf(hname, "hCaloPreShToF1");
0428   sprintf(htitle, "PreSh hits, ToF/ns");
0429   meCaloPreShToF[0] = iBooker.book1D(hname, htitle, 100, 0., 1000.);
0430   sprintf(hname, "hCaloPreShToF2");
0431   meCaloPreShToF[1] = iBooker.book1D(hname, htitle, 100, 0., 100.);
0432 
0433   for (Int_t i = 0; i < 2; ++i) {
0434     meCaloPreSh[i]->setAxisTitle("Number of Hits", 1);
0435     meCaloPreSh[i]->setAxisTitle("Count", 2);
0436     meCaloPreShE[i]->setAxisTitle("Energy of Hits (GeV)", 1);
0437     meCaloPreShE[i]->setAxisTitle("Count", 2);
0438     meCaloPreShToF[i]->setAxisTitle("Time of Flight of Hits (ns)", 1);
0439     meCaloPreShToF[i]->setAxisTitle("Count", 2);
0440   }
0441 
0442   sprintf(hname, "hCaloPreShPhi");
0443   sprintf(htitle, "PreSh hits, phi/rad");
0444   meCaloPreShPhi = iBooker.book1D(hname, htitle, 100, -3.2, 3.2);
0445   meCaloPreShPhi->setAxisTitle("Phi of Hits (rad)", 1);
0446   meCaloPreShPhi->setAxisTitle("Count", 2);
0447 
0448   sprintf(hname, "hCaloPreShEta");
0449   sprintf(htitle, "PreSh hits, eta");
0450   meCaloPreShEta = iBooker.book1D(hname, htitle, 100, -5.5, 5.5);
0451   meCaloPreShEta->setAxisTitle("Eta of Hits", 1);
0452   meCaloPreShEta->setAxisTitle("Count", 2);
0453 
0454   // Hcal
0455   iBooker.setCurrentFolder("GlobalHitsV/HCals");
0456   sprintf(hname, "hCaloHcal1");
0457   sprintf(htitle, "Hcal hits");
0458   meCaloHcal[0] = iBooker.book1D(hname, htitle, 100, 0., 10000.);
0459   sprintf(hname, "hCaloHcal2");
0460   meCaloHcal[1] = iBooker.book1D(hname, htitle, 100, -0.5, 99.5);
0461 
0462   sprintf(hname, "hCaloHcalE1");
0463   sprintf(htitle, "Hcal hits, energy/GeV");
0464   meCaloHcalE[0] = iBooker.book1D(hname, htitle, 100, 0., 10.);
0465   sprintf(hname, "hCaloHcalE2");
0466   meCaloHcalE[1] = iBooker.book1D(hname, htitle, 100, 0., 0.1);
0467 
0468   sprintf(hname, "hCaloHcalToF1");
0469   sprintf(htitle, "Hcal hits, ToF/ns");
0470   meCaloHcalToF[0] = iBooker.book1D(hname, htitle, 100, 0., 1000.);
0471   sprintf(hname, "hCaloHcalToF2");
0472   meCaloHcalToF[1] = iBooker.book1D(hname, htitle, 100, 0., 100.);
0473 
0474   for (Int_t i = 0; i < 2; ++i) {
0475     meCaloHcal[i]->setAxisTitle("Number of Hits", 1);
0476     meCaloHcal[i]->setAxisTitle("Count", 2);
0477     meCaloHcalE[i]->setAxisTitle("Energy of Hits (GeV)", 1);
0478     meCaloHcalE[i]->setAxisTitle("Count", 2);
0479     meCaloHcalToF[i]->setAxisTitle("Time of Flight of Hits (ns)", 1);
0480     meCaloHcalToF[i]->setAxisTitle("Count", 2);
0481   }
0482 
0483   sprintf(hname, "hCaloHcalPhi");
0484   sprintf(htitle, "Hcal hits, phi/rad");
0485   meCaloHcalPhi = iBooker.book1D(hname, htitle, 100, -3.2, 3.2);
0486   meCaloHcalPhi->setAxisTitle("Phi of Hits (rad)", 1);
0487   meCaloHcalPhi->setAxisTitle("Count", 2);
0488 
0489   sprintf(hname, "hCaloHcalEta");
0490   sprintf(htitle, "Hcal hits, eta");
0491   meCaloHcalEta = iBooker.book1D(hname, htitle, 100, -5.5, 5.5);
0492   meCaloHcalEta->setAxisTitle("Eta of Hits", 1);
0493   meCaloHcalEta->setAxisTitle("Count", 2);
0494 
0495   // SiPixels
0496   iBooker.setCurrentFolder("GlobalHitsV/SiPixels");
0497   sprintf(hname, "hTrackerPx1");
0498   sprintf(htitle, "Pixel hits");
0499   meTrackerPx[0] = iBooker.book1D(hname, htitle, 100, 0., 10000.);
0500   sprintf(hname, "hTrackerPx2");
0501   meTrackerPx[1] = iBooker.book1D(hname, htitle, 100, -0.5, 99.5);
0502   for (Int_t i = 0; i < 2; ++i) {
0503     meTrackerPx[i]->setAxisTitle("Number of Pixel Hits", 1);
0504     meTrackerPx[i]->setAxisTitle("Count", 2);
0505   }
0506 
0507   sprintf(hname, "hTrackerPxPhi");
0508   sprintf(htitle, "Pixel hits phi/rad");
0509   meTrackerPxPhi = iBooker.book1D(hname, htitle, 100, -3.2, 3.2);
0510   meTrackerPxPhi->setAxisTitle("Phi of Hits (rad)", 1);
0511   meTrackerPxPhi->setAxisTitle("Count", 2);
0512 
0513   sprintf(hname, "hTrackerPxEta");
0514   sprintf(htitle, "Pixel hits eta");
0515   meTrackerPxEta = iBooker.book1D(hname, htitle, 100, -3.5, 3.5);
0516   meTrackerPxEta->setAxisTitle("Eta of Hits", 1);
0517   meTrackerPxEta->setAxisTitle("Count", 2);
0518 
0519   sprintf(hname, "hTrackerPxBToF");
0520   sprintf(htitle, "Pixel barrel hits, ToF/ns");
0521   meTrackerPxBToF = iBooker.book1D(hname, htitle, 100, 0., 40.);
0522   meTrackerPxBToF->setAxisTitle("Time of Flight of Hits (ns)", 1);
0523   meTrackerPxBToF->setAxisTitle("Count", 2);
0524 
0525   sprintf(hname, "hTrackerPxBR");
0526   sprintf(htitle, "Pixel barrel hits, R/cm");
0527   meTrackerPxBR = iBooker.book1D(hname, htitle, 100, 0., 50.);
0528   meTrackerPxBR->setAxisTitle("R of Hits (cm)", 1);
0529   meTrackerPxBR->setAxisTitle("Count", 2);
0530 
0531   sprintf(hname, "hTrackerPxFToF");
0532   sprintf(htitle, "Pixel forward hits, ToF/ns");
0533   meTrackerPxFToF = iBooker.book1D(hname, htitle, 100, 0., 50.);
0534   meTrackerPxFToF->setAxisTitle("Time of Flight of Hits (ns)", 1);
0535   meTrackerPxFToF->setAxisTitle("Count", 2);
0536 
0537   sprintf(hname, "hTrackerPxFZ");
0538   sprintf(htitle, "Pixel forward hits, Z/cm");
0539   meTrackerPxFZ = iBooker.book1D(hname, htitle, 200, -100., 100.);
0540   meTrackerPxFZ->setAxisTitle("Z of Hits (cm)", 1);
0541   meTrackerPxFZ->setAxisTitle("Count", 2);
0542 
0543   // SiStrips
0544   iBooker.setCurrentFolder("GlobalHitsV/SiStrips");
0545   sprintf(hname, "hTrackerSi1");
0546   sprintf(htitle, "Silicon hits");
0547   meTrackerSi[0] = iBooker.book1D(hname, htitle, 100, 0., 10000.);
0548   sprintf(hname, "hTrackerSi2");
0549   meTrackerSi[1] = iBooker.book1D(hname, htitle, 100, -0.5, 99.5);
0550   for (Int_t i = 0; i < 2; ++i) {
0551     meTrackerSi[i]->setAxisTitle("Number of Silicon Hits", 1);
0552     meTrackerSi[i]->setAxisTitle("Count", 2);
0553   }
0554 
0555   sprintf(hname, "hTrackerSiPhi");
0556   sprintf(htitle, "Silicon hits phi/rad");
0557   meTrackerSiPhi = iBooker.book1D(hname, htitle, 100, -3.2, 3.2);
0558   meTrackerSiPhi->setAxisTitle("Phi of Hits (rad)", 1);
0559   meTrackerSiPhi->setAxisTitle("Count", 2);
0560 
0561   sprintf(hname, "hTrackerSiEta");
0562   sprintf(htitle, "Silicon hits eta");
0563   meTrackerSiEta = iBooker.book1D(hname, htitle, 100, -3.5, 3.5);
0564   meTrackerSiEta->setAxisTitle("Eta of Hits", 1);
0565   meTrackerSiEta->setAxisTitle("Count", 2);
0566 
0567   sprintf(hname, "hTrackerSiBToF");
0568   sprintf(htitle, "Silicon barrel hits, ToF/ns");
0569   meTrackerSiBToF = iBooker.book1D(hname, htitle, 100, 0., 50.);
0570   meTrackerSiBToF->setAxisTitle("Time of Flight of Hits (ns)", 1);
0571   meTrackerSiBToF->setAxisTitle("Count", 2);
0572 
0573   sprintf(hname, "hTrackerSiBR");
0574   sprintf(htitle, "Silicon barrel hits, R/cm");
0575   meTrackerSiBR = iBooker.book1D(hname, htitle, 100, 0., 200.);
0576   meTrackerSiBR->setAxisTitle("R of Hits (cm)", 1);
0577   meTrackerSiBR->setAxisTitle("Count", 2);
0578 
0579   sprintf(hname, "hTrackerSiFToF");
0580   sprintf(htitle, "Silicon forward hits, ToF/ns");
0581   meTrackerSiFToF = iBooker.book1D(hname, htitle, 100, 0., 75.);
0582   meTrackerSiFToF->setAxisTitle("Time of Flight of Hits (ns)", 1);
0583   meTrackerSiFToF->setAxisTitle("Count", 2);
0584 
0585   sprintf(hname, "hTrackerSiFZ");
0586   sprintf(htitle, "Silicon forward hits, Z/cm");
0587   meTrackerSiFZ = iBooker.book1D(hname, htitle, 200, -300., 300.);
0588   meTrackerSiFZ->setAxisTitle("Z of Hits (cm)", 1);
0589   meTrackerSiFZ->setAxisTitle("Count", 2);
0590 
0591   // Muon
0592   iBooker.setCurrentFolder("GlobalHitsV/Muons");
0593   sprintf(hname, "hMuon1");
0594   sprintf(htitle, "Muon hits");
0595   meMuon[0] = iBooker.book1D(hname, htitle, 100, 0., 10000.);
0596   sprintf(hname, "hMuon2");
0597   meMuon[1] = iBooker.book1D(hname, htitle, 100, -0.5, 99.5);
0598   for (Int_t i = 0; i < 2; ++i) {
0599     meMuon[i]->setAxisTitle("Number of Muon Hits", 1);
0600     meMuon[i]->setAxisTitle("Count", 2);
0601   }
0602 
0603   sprintf(hname, "hMuonPhi");
0604   sprintf(htitle, "Muon hits phi/rad");
0605   meMuonPhi = iBooker.book1D(hname, htitle, 100, -3.2, 3.2);
0606   meMuonPhi->setAxisTitle("Phi of Hits (rad)", 1);
0607   meMuonPhi->setAxisTitle("Count", 2);
0608 
0609   sprintf(hname, "hMuonEta");
0610   sprintf(htitle, "Muon hits eta");
0611   meMuonEta = iBooker.book1D(hname, htitle, 100, -3.5, 3.5);
0612   meMuonEta->setAxisTitle("Eta of Hits", 1);
0613   meMuonEta->setAxisTitle("Count", 2);
0614 
0615   sprintf(hname, "hMuonCscToF1");
0616   sprintf(htitle, "Muon CSC hits, ToF/ns");
0617   meMuonCscToF[0] = iBooker.book1D(hname, htitle, 100, 0., 250.);
0618   sprintf(hname, "hMuonCscToF2");
0619   meMuonCscToF[1] = iBooker.book1D(hname, htitle, 100, 0., 50.);
0620   for (Int_t i = 0; i < 2; ++i) {
0621     meMuonCscToF[i]->setAxisTitle("Time of Flight of Hits (ns)", 1);
0622     meMuonCscToF[i]->setAxisTitle("Count", 2);
0623   }
0624 
0625   sprintf(hname, "hMuonCscZ");
0626   sprintf(htitle, "Muon CSC hits, Z/cm");
0627   meMuonCscZ = iBooker.book1D(hname, htitle, 200, -1500., 1500.);
0628   meMuonCscZ->setAxisTitle("Z of Hits (cm)", 1);
0629   meMuonCscZ->setAxisTitle("Count", 2);
0630 
0631   sprintf(hname, "hMuonDtToF1");
0632   sprintf(htitle, "Muon DT hits, ToF/ns");
0633   meMuonDtToF[0] = iBooker.book1D(hname, htitle, 100, 0., 250.);
0634   sprintf(hname, "hMuonDtToF2");
0635   meMuonDtToF[1] = iBooker.book1D(hname, htitle, 100, 0., 50.);
0636   for (Int_t i = 0; i < 2; ++i) {
0637     meMuonDtToF[i]->setAxisTitle("Time of Flight of Hits (ns)", 1);
0638     meMuonDtToF[i]->setAxisTitle("Count", 2);
0639   }
0640 
0641   sprintf(hname, "hMuonDtR");
0642   sprintf(htitle, "Muon DT hits, R/cm");
0643   meMuonDtR = iBooker.book1D(hname, htitle, 100, 0., 1500.);
0644   meMuonDtR->setAxisTitle("R of Hits (cm)", 1);
0645   meMuonDtR->setAxisTitle("Count", 2);
0646 
0647   sprintf(hname, "hMuonRpcFToF1");
0648   sprintf(htitle, "Muon RPC forward hits, ToF/ns");
0649   meMuonRpcFToF[0] = iBooker.book1D(hname, htitle, 100, 0., 250.);
0650   sprintf(hname, "hMuonRpcFToF2");
0651   meMuonRpcFToF[1] = iBooker.book1D(hname, htitle, 100, 0., 50.);
0652   for (Int_t i = 0; i < 2; ++i) {
0653     meMuonRpcFToF[i]->setAxisTitle("Time of Flight of Hits (ns)", 1);
0654     meMuonRpcFToF[i]->setAxisTitle("Count", 2);
0655   }
0656   sprintf(hname, "hMuonRpcFZ");
0657   sprintf(htitle, "Muon RPC forward hits, Z/cm");
0658   meMuonRpcFZ = iBooker.book1D(hname, htitle, 201, -1500., 1500.);
0659   meMuonRpcFZ->setAxisTitle("Z of Hits (cm)", 1);
0660   meMuonRpcFZ->setAxisTitle("Count", 2);
0661 
0662   sprintf(hname, "hMuonRpcBToF1");
0663   sprintf(htitle, "Muon RPC barrel hits, ToF/ns");
0664   meMuonRpcBToF[0] = iBooker.book1D(hname, htitle, 100, 0., 250.);
0665   sprintf(hname, "hMuonRpcBToF2");
0666   meMuonRpcBToF[1] = iBooker.book1D(hname, htitle, 100, 0., 50.);
0667   for (Int_t i = 0; i < 2; ++i) {
0668     meMuonRpcBToF[i]->setAxisTitle("Time of Flight of Hits (ns)", 1);
0669     meMuonRpcBToF[i]->setAxisTitle("Count", 2);
0670   }
0671 
0672   sprintf(hname, "hMuonRpcBR");
0673   sprintf(htitle, "Muon RPC barrel hits, R/cm");
0674   meMuonRpcBR = iBooker.book1D(hname, htitle, 100, 0., 1500.);
0675   meMuonRpcBR->setAxisTitle("R of Hits (cm)", 1);
0676   meMuonRpcBR->setAxisTitle("Count", 2);
0677 }
0678 
0679 void GlobalHitsAnalyzer::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0680   std::string MsgLoggerCat = "GlobalHitsAnalyzer_analyze";
0681 
0682   // keep track of number of events processed
0683   ++count;
0684 
0685   // get event id information
0686   edm::RunNumber_t nrun = iEvent.id().run();
0687   edm::EventNumber_t nevt = iEvent.id().event();
0688 
0689   if (verbosity > 0) {
0690     edm::LogInfo(MsgLoggerCat) << "Processing run " << nrun << ", event " << nevt << " (" << count << " events total)";
0691   } else if (verbosity == 0) {
0692     if (nevt % frequency == 0 || nevt == 1) {
0693       edm::LogInfo(MsgLoggerCat) << "Processing run " << nrun << ", event " << nevt << " (" << count
0694                                  << " events total)";
0695     }
0696   }
0697 
0698   // look at information available in the event
0699   if (getAllProvenances) {
0700     std::vector<const edm::StableProvenance *> AllProv;
0701     iEvent.getAllStableProvenance(AllProv);
0702 
0703     if (verbosity >= 0)
0704       edm::LogInfo(MsgLoggerCat) << "Number of Provenances = " << AllProv.size();
0705 
0706     if (printProvenanceInfo && (verbosity >= 0)) {
0707       TString eventout("\nProvenance info:\n");
0708 
0709       for (unsigned int i = 0; i < AllProv.size(); ++i) {
0710         eventout += "\n       ******************************";
0711         eventout += "\n       Module       : ";
0712         eventout += AllProv[i]->moduleLabel();
0713         eventout += "\n       ProductID    : ";
0714         eventout += AllProv[i]->productID().id();
0715         eventout += "\n       ClassName    : ";
0716         eventout += AllProv[i]->className();
0717         eventout += "\n       InstanceName : ";
0718         eventout += AllProv[i]->productInstanceName();
0719         eventout += "\n       BranchName   : ";
0720         eventout += AllProv[i]->branchName();
0721       }
0722       eventout += "\n       ******************************\n";
0723       edm::LogInfo(MsgLoggerCat) << eventout << "\n";
0724       printProvenanceInfo = false;
0725     }
0726     getAllProvenances = false;
0727   }
0728 
0729   // call fill functions
0730   // gather G4MC information from event
0731   fillG4MC(iEvent);
0732   // gather Tracker information from event
0733   fillTrk(iEvent, iSetup);
0734   // gather muon information from event
0735   fillMuon(iEvent, iSetup);
0736   // gather Ecal information from event
0737   fillECal(iEvent, iSetup);
0738   // gather Hcal information from event
0739   fillHCal(iEvent, iSetup);
0740 
0741   if (verbosity > 0)
0742     edm::LogInfo(MsgLoggerCat) << "Done gathering data from event.";
0743 
0744   return;
0745 }
0746 
0747 //==================fill and store functions================================
0748 void GlobalHitsAnalyzer::fillG4MC(const edm::Event &iEvent) {
0749   std::string MsgLoggerCat = "GlobalHitsAnalyzer_fillG4MC";
0750 
0751   TString eventout;
0752   if (verbosity > 0)
0753     eventout = "\nGathering info:";
0754 
0755   //////////////////////
0756   // get MC information
0757   /////////////////////
0758   edm::Handle<edm::HepMCProduct> HepMCEvt;
0759   std::vector<edm::Handle<edm::HepMCProduct>> AllHepMCEvt;
0760   getterOfProducts_.fillHandles(iEvent, AllHepMCEvt);
0761   // loop through products and extract VtxSmearing if available. Any of them
0762   // should have the information needed
0763   for (unsigned int i = 0; i < AllHepMCEvt.size(); ++i) {
0764     HepMCEvt = AllHepMCEvt[i];
0765     if ((HepMCEvt.provenance()->branchDescription()).moduleLabel() == "generatorSmeared")
0766       break;
0767   }
0768 
0769   if (!HepMCEvt.isValid()) {
0770     LogDebug(MsgLoggerCat) << "Unable to find HepMCProduct in event!";
0771     validHepMCevt = false;
0772   } else {
0773     eventout += "\n          Using HepMCProduct: ";
0774     eventout += (HepMCEvt.provenance()->branchDescription()).moduleLabel();
0775   }
0776   if (validHepMCevt) {
0777     const HepMC::GenEvent *MCEvt = HepMCEvt->GetEvent();
0778     nRawGenPart = MCEvt->particles_size();
0779 
0780     if (verbosity > 1) {
0781       eventout += "\n          Number of Raw Particles collected:......... ";
0782       eventout += nRawGenPart;
0783     }
0784 
0785     if (meMCRGP[0])
0786       meMCRGP[0]->Fill((float)nRawGenPart);
0787     if (meMCRGP[1])
0788       meMCRGP[1]->Fill((float)nRawGenPart);
0789   }
0790 
0791   ////////////////////////////
0792   // get G4Vertex information
0793   ////////////////////////////
0794   // convert unit stored in SimVertex to mm
0795   float unit = 0.;
0796   if (vtxunit == 0)
0797     unit = 1.;  // already in mm
0798   if (vtxunit == 1)
0799     unit = 10.;  // stored in cm, convert to mm
0800 
0801   edm::Handle<edm::SimVertexContainer> G4VtxContainer;
0802   iEvent.getByToken(G4VtxSrc_Token_, G4VtxContainer);
0803 
0804   // needed here by vertex multiplicity
0805   edm::Handle<edm::SimTrackContainer> G4TrkContainer;
0806   iEvent.getByToken(G4TrkSrc_Token_, G4TrkContainer);
0807 
0808   if (!G4VtxContainer.isValid()) {
0809     LogDebug(MsgLoggerCat) << "Unable to find SimVertex in event!";
0810     validG4VtxContainer = false;
0811   }
0812   if (validG4VtxContainer) {
0813     int i = 0;
0814     edm::SimVertexContainer::const_iterator itVtx;
0815     for (itVtx = G4VtxContainer->begin(); itVtx != G4VtxContainer->end(); ++itVtx) {
0816       ++i;
0817 
0818       const math::XYZTLorentzVector G4Vtx1(
0819           itVtx->position().x(), itVtx->position().y(), itVtx->position().z(), itVtx->position().e());
0820 
0821       double G4Vtx[4];
0822       G4Vtx1.GetCoordinates(G4Vtx);
0823 
0824       if (meGeantVtxX[0])
0825         meGeantVtxX[0]->Fill((G4Vtx[0] * unit) / micrometer);
0826       if (meGeantVtxX[1])
0827         meGeantVtxX[1]->Fill((G4Vtx[0] * unit) / micrometer);
0828 
0829       if (meGeantVtxY[0])
0830         meGeantVtxY[0]->Fill((G4Vtx[1] * unit) / micrometer);
0831       if (meGeantVtxY[1])
0832         meGeantVtxY[1]->Fill((G4Vtx[1] * unit) / micrometer);
0833 
0834       if (meGeantVtxZ[0])
0835         meGeantVtxZ[0]->Fill((G4Vtx[2] * unit) / millimeter);
0836       if (meGeantVtxZ[1])
0837         meGeantVtxZ[1]->Fill((G4Vtx[2] * unit) / millimeter);
0838 
0839       if (meGeantVtxEta)
0840         meGeantVtxEta->Fill(G4Vtx1.eta());
0841       if (meGeantVtxPhi)
0842         meGeantVtxPhi->Fill(G4Vtx1.phi());
0843       if (meGeantVtxRad[0])
0844         meGeantVtxRad[0]->Fill(G4Vtx1.rho());
0845       if (meGeantVtxRad[1])
0846         meGeantVtxRad[1]->Fill(G4Vtx1.rho());
0847 
0848       if (meGeantVtxMulti) {
0849         int multi = 0;
0850         if (G4TrkContainer.isValid()) {
0851           edm::SimTrackContainer::const_iterator itTrk;
0852           for (itTrk = G4TrkContainer->begin(); itTrk != G4TrkContainer->end(); ++itTrk) {
0853             if ((*itTrk).vertIndex() == i) {
0854               multi++;
0855             }
0856           }
0857         }
0858         meGeantVtxMulti->Fill(((double)multi + 0.5));
0859       }
0860     }
0861 
0862     if (verbosity > 1) {
0863       eventout += "\n          Number of G4Vertices collected:............ ";
0864       eventout += i;
0865     }
0866 
0867     if (meMCG4Vtx[0])
0868       meMCG4Vtx[0]->Fill((float)i);
0869     if (meMCG4Vtx[1])
0870       meMCG4Vtx[1]->Fill((float)i);
0871   }
0872 
0873   ///////////////////////////
0874   // get G4Track information
0875   ///////////////////////////
0876   if (!G4TrkContainer.isValid()) {
0877     LogDebug(MsgLoggerCat) << "Unable to find SimTrack in event!";
0878     validG4trkContainer = false;
0879   }
0880   if (validG4trkContainer) {
0881     int i = 0;
0882     edm::SimTrackContainer::const_iterator itTrk;
0883     for (itTrk = G4TrkContainer->begin(); itTrk != G4TrkContainer->end(); ++itTrk) {
0884       ++i;
0885 
0886       const math::XYZTLorentzVector G4Trk1(
0887           itTrk->momentum().x(), itTrk->momentum().y(), itTrk->momentum().z(), itTrk->momentum().e());
0888       double G4Trk[4];
0889       G4Trk1.GetCoordinates(G4Trk);
0890 
0891       if (meGeantTrkPt)
0892         meGeantTrkPt->Fill(std::log10(std::max(sqrt(G4Trk[0] * G4Trk[0] + G4Trk[1] * G4Trk[1]), -9.)));
0893       if (meGeantTrkE)
0894         meGeantTrkE->Fill(std::log10(std::max(G4Trk[3], -9.)));
0895     }
0896 
0897     if (verbosity > 1) {
0898       eventout += "\n          Number of G4Tracks collected:.............. ";
0899       eventout += i;
0900     }
0901 
0902     if (meMCG4Trk[0])
0903       meMCG4Trk[0]->Fill((float)i);
0904     if (meMCG4Trk[1])
0905       meMCG4Trk[1]->Fill((float)i);
0906   }
0907 
0908   if (verbosity > 0)
0909     edm::LogInfo(MsgLoggerCat) << eventout << "\n";
0910 
0911   return;
0912 }
0913 
0914 void GlobalHitsAnalyzer::fillTrk(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0915   nPxlHits = 0;
0916   std::string MsgLoggerCat = "GlobalHitsAnalyzer_fillTrk";
0917 
0918   TString eventout;
0919   if (verbosity > 0)
0920     eventout = "\nGathering info:";
0921 
0922   // access the tracker geometry
0923   const auto &theTrackerGeometry = iSetup.getHandle(tGeomToken_);
0924   if (!theTrackerGeometry.isValid()) {
0925     edm::LogWarning(MsgLoggerCat) << "Unable to find TrackerDigiGeometryRecord in event!";
0926     return;
0927   }
0928   const TrackerGeometry &theTracker(*theTrackerGeometry);
0929 
0930   // iterator to access containers
0931   edm::PSimHitContainer::const_iterator itHit;
0932 
0933   ///////////////////////////////
0934   // get Pixel Barrel information
0935   ///////////////////////////////
0936   edm::PSimHitContainer thePxlBrlHits;
0937   // extract low container
0938   edm::Handle<edm::PSimHitContainer> PxlBrlLowContainer;
0939   iEvent.getByToken(PxlBrlLowSrc_Token_, PxlBrlLowContainer);
0940   if (!PxlBrlLowContainer.isValid()) {
0941     LogDebug(MsgLoggerCat) << "Unable to find TrackerHitsPixelBarrelLowTof in event!";
0942     validPxlBrlLow = false;
0943   }
0944   // extract high container
0945   edm::Handle<edm::PSimHitContainer> PxlBrlHighContainer;
0946   iEvent.getByToken(PxlBrlHighSrc_Token_, PxlBrlHighContainer);
0947   if (!PxlBrlHighContainer.isValid()) {
0948     LogDebug(MsgLoggerCat) << "Unable to find TrackerHitsPixelBarrelHighTof in event!";
0949     validPxlBrlHigh = false;
0950   }
0951   // place both containers into new container
0952   if (validPxlBrlLow)
0953     thePxlBrlHits.insert(thePxlBrlHits.end(), PxlBrlLowContainer->begin(), PxlBrlLowContainer->end());
0954   if (validPxlBrlHigh)
0955     thePxlBrlHits.insert(thePxlBrlHits.end(), PxlBrlHighContainer->begin(), PxlBrlHighContainer->end());
0956 
0957   // cycle through new container
0958   int i = 0, j = 0;
0959   for (itHit = thePxlBrlHits.begin(); itHit != thePxlBrlHits.end(); ++itHit) {
0960     ++i;
0961 
0962     // create a DetId from the detUnitId
0963     DetId theDetUnitId(itHit->detUnitId());
0964     int detector = theDetUnitId.det();
0965     int subdetector = theDetUnitId.subdetId();
0966 
0967     // check that expected detector is returned
0968     if ((detector == dTrk) && (subdetector == sdPxlBrl)) {
0969       // get the GeomDetUnit from the geometry using theDetUnitID
0970       const GeomDetUnit *theDet = theTracker.idToDetUnit(theDetUnitId);
0971 
0972       if (!theDet) {
0973         edm::LogWarning(MsgLoggerCat) << "Unable to get GeomDetUnit from PxlBrlHits for Hit " << i;
0974         continue;
0975       }
0976 
0977       ++j;
0978 
0979       // get the Surface of the hit (knows how to go from local <-> global)
0980       const BoundPlane &bSurface = theDet->surface();
0981 
0982       if (meTrackerPxBToF)
0983         meTrackerPxBToF->Fill(itHit->tof());
0984       if (meTrackerPxBR)
0985         meTrackerPxBR->Fill(bSurface.toGlobal(itHit->localPosition()).perp());
0986       if (meTrackerPxPhi)
0987         meTrackerPxPhi->Fill(bSurface.toGlobal(itHit->localPosition()).phi());
0988       if (meTrackerPxEta)
0989         meTrackerPxEta->Fill(bSurface.toGlobal(itHit->localPosition()).eta());
0990 
0991     } else {
0992       edm::LogWarning(MsgLoggerCat) << "PxlBrl PSimHit " << i << " is expected to be (det,subdet) = (" << dTrk << ","
0993                                     << sdPxlBrl << "); value returned is: (" << detector << "," << subdetector << ")";
0994       continue;
0995     }  // end detector type check
0996   }    // end loop through PxlBrl Hits
0997 
0998   if (verbosity > 1) {
0999     eventout += "\n          Number of Pixel Barrel Hits collected:..... ";
1000     eventout += j;
1001   }
1002 
1003   nPxlHits += j;
1004 
1005   /////////////////////////////////
1006   // get Pixel Forward information
1007   ////////////////////////////////
1008   edm::PSimHitContainer thePxlFwdHits;
1009   // extract low container
1010   edm::Handle<edm::PSimHitContainer> PxlFwdLowContainer;
1011   iEvent.getByToken(PxlFwdLowSrc_Token_, PxlFwdLowContainer);
1012   if (!PxlFwdLowContainer.isValid()) {
1013     LogDebug(MsgLoggerCat) << "Unable to find TrackerHitsPixelEndcapLowTof in event!";
1014     validPxlFwdLow = false;
1015   }
1016   // extract high container
1017   edm::Handle<edm::PSimHitContainer> PxlFwdHighContainer;
1018   iEvent.getByToken(PxlFwdHighSrc_Token_, PxlFwdHighContainer);
1019   if (!PxlFwdHighContainer.isValid()) {
1020     LogDebug("GlobalHitsAnalyzer_fillTrk") << "Unable to find TrackerHitsPixelEndcapHighTof in event!";
1021     validPxlFwdHigh = false;
1022   }
1023   // place both containers into new container
1024   if (validPxlFwdLow)
1025     thePxlFwdHits.insert(thePxlFwdHits.end(), PxlFwdLowContainer->begin(), PxlFwdLowContainer->end());
1026   if (validPxlFwdHigh)
1027     thePxlFwdHits.insert(thePxlFwdHits.end(), PxlFwdHighContainer->begin(), PxlFwdHighContainer->end());
1028 
1029   // cycle through new container
1030   i = 0;
1031   j = 0;
1032   for (itHit = thePxlFwdHits.begin(); itHit != thePxlFwdHits.end(); ++itHit) {
1033     ++i;
1034 
1035     // create a DetId from the detUnitId
1036     DetId theDetUnitId(itHit->detUnitId());
1037     int detector = theDetUnitId.det();
1038     int subdetector = theDetUnitId.subdetId();
1039 
1040     // check that expected detector is returned
1041     if ((detector == dTrk) && (subdetector == sdPxlFwd)) {
1042       // get the GeomDetUnit from the geometry using theDetUnitID
1043       const GeomDetUnit *theDet = theTracker.idToDetUnit(theDetUnitId);
1044 
1045       if (!theDet) {
1046         edm::LogWarning(MsgLoggerCat) << "Unable to get GeomDetUnit from PxlFwdHits for Hit " << i;
1047         ;
1048         continue;
1049       }
1050 
1051       ++j;
1052 
1053       // get the Surface of the hit (knows how to go from local <-> global)
1054       const BoundPlane &bSurface = theDet->surface();
1055 
1056       if (meTrackerPxFToF)
1057         meTrackerPxFToF->Fill(itHit->tof());
1058       if (meTrackerPxFZ)
1059         meTrackerPxFZ->Fill(bSurface.toGlobal(itHit->localPosition()).z());
1060       if (meTrackerPxPhi)
1061         meTrackerPxPhi->Fill(bSurface.toGlobal(itHit->localPosition()).phi());
1062       if (meTrackerPxEta)
1063         meTrackerPxEta->Fill(bSurface.toGlobal(itHit->localPosition()).eta());
1064 
1065     } else {
1066       edm::LogWarning(MsgLoggerCat) << "PxlFwd PSimHit " << i << " is expected to be (det,subdet) = (" << dTrk << ","
1067                                     << sdPxlFwd << "); value returned is: (" << detector << "," << subdetector << ")";
1068       continue;
1069     }  // end detector type check
1070   }    // end loop through PxlFwd Hits
1071 
1072   if (verbosity > 1) {
1073     eventout += "\n          Number of Pixel Forward Hits collected:.... ";
1074     eventout += j;
1075   }
1076 
1077   nPxlHits += j;
1078 
1079   if (meTrackerPx[0])
1080     meTrackerPx[0]->Fill((float)nPxlHits);
1081   if (meTrackerPx[1])
1082     meTrackerPx[1]->Fill((float)nPxlHits);
1083 
1084   ///////////////////////////////////
1085   // get Silicon Barrel information
1086   //////////////////////////////////
1087   nSiHits = 0;
1088   edm::PSimHitContainer theSiBrlHits;
1089   // extract TIB low container
1090   edm::Handle<edm::PSimHitContainer> SiTIBLowContainer;
1091   iEvent.getByToken(SiTIBLowSrc_Token_, SiTIBLowContainer);
1092   if (!SiTIBLowContainer.isValid()) {
1093     LogDebug(MsgLoggerCat) << "Unable to find TrackerHitsTIBLowTof in event!";
1094     validSiTIBLow = false;
1095   }
1096   // extract TIB high container
1097   edm::Handle<edm::PSimHitContainer> SiTIBHighContainer;
1098   iEvent.getByToken(SiTIBHighSrc_Token_, SiTIBHighContainer);
1099   if (!SiTIBHighContainer.isValid()) {
1100     LogDebug(MsgLoggerCat) << "Unable to find TrackerHitsTIBHighTof in event!";
1101     validSiTIBHigh = false;
1102   }
1103   // extract TOB low container
1104   edm::Handle<edm::PSimHitContainer> SiTOBLowContainer;
1105   iEvent.getByToken(SiTOBLowSrc_Token_, SiTOBLowContainer);
1106   if (!SiTOBLowContainer.isValid()) {
1107     LogDebug(MsgLoggerCat) << "Unable to find TrackerHitsTOBLowTof in event!";
1108     validSiTOBLow = false;
1109   }
1110   // extract TOB high container
1111   edm::Handle<edm::PSimHitContainer> SiTOBHighContainer;
1112   iEvent.getByToken(SiTOBHighSrc_Token_, SiTOBHighContainer);
1113   if (!SiTOBHighContainer.isValid()) {
1114     LogDebug(MsgLoggerCat) << "Unable to find TrackerHitsTOBHighTof in event!";
1115     validSiTOBHigh = false;
1116   }
1117   // place all containers into new container
1118   if (validSiTIBLow)
1119     theSiBrlHits.insert(theSiBrlHits.end(), SiTIBLowContainer->begin(), SiTIBLowContainer->end());
1120   if (validSiTIBHigh)
1121     theSiBrlHits.insert(theSiBrlHits.end(), SiTIBHighContainer->begin(), SiTIBHighContainer->end());
1122   if (validSiTOBLow)
1123     theSiBrlHits.insert(theSiBrlHits.end(), SiTOBLowContainer->begin(), SiTOBLowContainer->end());
1124   if (validSiTOBHigh)
1125     theSiBrlHits.insert(theSiBrlHits.end(), SiTOBHighContainer->begin(), SiTOBHighContainer->end());
1126 
1127   // cycle through new container
1128   i = 0;
1129   j = 0;
1130   for (itHit = theSiBrlHits.begin(); itHit != theSiBrlHits.end(); ++itHit) {
1131     ++i;
1132 
1133     // create a DetId from the detUnitId
1134     DetId theDetUnitId(itHit->detUnitId());
1135     int detector = theDetUnitId.det();
1136     int subdetector = theDetUnitId.subdetId();
1137 
1138     // check that expected detector is returned
1139     if ((detector == dTrk) && ((subdetector == sdSiTIB) || (subdetector == sdSiTOB))) {
1140       // get the GeomDetUnit from the geometry using theDetUnitID
1141       const GeomDetUnit *theDet = theTracker.idToDetUnit(theDetUnitId);
1142 
1143       if (!theDet) {
1144         edm::LogWarning(MsgLoggerCat) << "Unable to get GeomDetUnit from SiBrlHits for Hit " << i;
1145         continue;
1146       }
1147 
1148       ++j;
1149 
1150       // get the Surface of the hit (knows how to go from local <-> global)
1151       const BoundPlane &bSurface = theDet->surface();
1152 
1153       if (meTrackerSiBToF)
1154         meTrackerSiBToF->Fill(itHit->tof());
1155       if (meTrackerSiBR)
1156         meTrackerSiBR->Fill(bSurface.toGlobal(itHit->localPosition()).perp());
1157       if (meTrackerSiPhi)
1158         meTrackerSiPhi->Fill(bSurface.toGlobal(itHit->localPosition()).phi());
1159       if (meTrackerSiEta)
1160         meTrackerSiEta->Fill(bSurface.toGlobal(itHit->localPosition()).eta());
1161 
1162     } else {
1163       edm::LogWarning(MsgLoggerCat) << "SiBrl PSimHit " << i << " is expected to be (det,subdet) = (" << dTrk << ","
1164                                     << sdSiTIB << " || " << sdSiTOB << "); value returned is: (" << detector << ","
1165                                     << subdetector << ")";
1166       continue;
1167     }  // end detector type check
1168   }    // end loop through SiBrl Hits
1169 
1170   if (verbosity > 1) {
1171     eventout += "\n          Number of Silicon Barrel Hits collected:... ";
1172     eventout += j;
1173   }
1174 
1175   nSiHits += j;
1176 
1177   ///////////////////////////////////
1178   // get Silicon Forward information
1179   ///////////////////////////////////
1180   edm::PSimHitContainer theSiFwdHits;
1181   // extract TID low container
1182   edm::Handle<edm::PSimHitContainer> SiTIDLowContainer;
1183   iEvent.getByToken(SiTIDLowSrc_Token_, SiTIDLowContainer);
1184   if (!SiTIDLowContainer.isValid()) {
1185     LogDebug(MsgLoggerCat) << "Unable to find TrackerHitsTIDLowTof in event!";
1186     validSiTIDLow = false;
1187   }
1188   // extract TID high container
1189   edm::Handle<edm::PSimHitContainer> SiTIDHighContainer;
1190   iEvent.getByToken(SiTIDHighSrc_Token_, SiTIDHighContainer);
1191   if (!SiTIDHighContainer.isValid()) {
1192     LogDebug("GlobalHitsAnalyzer_fillTrk") << "Unable to find TrackerHitsTIDHighTof in event!";
1193     validSiTIDHigh = false;
1194   }
1195   // extract TEC low container
1196   edm::Handle<edm::PSimHitContainer> SiTECLowContainer;
1197   iEvent.getByToken(SiTECLowSrc_Token_, SiTECLowContainer);
1198   if (!SiTECLowContainer.isValid()) {
1199     LogDebug(MsgLoggerCat) << "Unable to find TrackerHitsTECLowTof in event!";
1200     validSiTECLow = false;
1201   }
1202   // extract TEC high container
1203   edm::Handle<edm::PSimHitContainer> SiTECHighContainer;
1204   iEvent.getByToken(SiTECHighSrc_Token_, SiTECHighContainer);
1205   if (!SiTECHighContainer.isValid()) {
1206     LogDebug(MsgLoggerCat) << "Unable to find TrackerHitsTECHighTof in event!";
1207     validSiTECHigh = false;
1208   }
1209   // place all containers into new container
1210   if (validSiTIDLow)
1211     theSiFwdHits.insert(theSiFwdHits.end(), SiTIDLowContainer->begin(), SiTIDLowContainer->end());
1212   if (validSiTIDHigh)
1213     theSiFwdHits.insert(theSiFwdHits.end(), SiTIDHighContainer->begin(), SiTIDHighContainer->end());
1214   if (validSiTECLow)
1215     theSiFwdHits.insert(theSiFwdHits.end(), SiTECLowContainer->begin(), SiTECLowContainer->end());
1216   if (validSiTECHigh)
1217     theSiFwdHits.insert(theSiFwdHits.end(), SiTECHighContainer->begin(), SiTECHighContainer->end());
1218 
1219   // cycle through container
1220   i = 0;
1221   j = 0;
1222   for (itHit = theSiFwdHits.begin(); itHit != theSiFwdHits.end(); ++itHit) {
1223     ++i;
1224 
1225     // create a DetId from the detUnitId
1226     DetId theDetUnitId(itHit->detUnitId());
1227     int detector = theDetUnitId.det();
1228     int subdetector = theDetUnitId.subdetId();
1229 
1230     // check that expected detector is returned
1231     if ((detector == dTrk) && ((subdetector == sdSiTID) || (subdetector == sdSiTEC))) {
1232       // get the GeomDetUnit from the geometry using theDetUnitID
1233       const GeomDetUnit *theDet = theTracker.idToDetUnit(theDetUnitId);
1234 
1235       if (!theDet) {
1236         edm::LogWarning(MsgLoggerCat) << "Unable to get GeomDetUnit from SiFwdHits Hit " << i;
1237         return;
1238       }
1239 
1240       ++j;
1241 
1242       // get the Surface of the hit (knows how to go from local <-> global)
1243       const BoundPlane &bSurface = theDet->surface();
1244 
1245       if (meTrackerSiFToF)
1246         meTrackerSiFToF->Fill(itHit->tof());
1247       if (meTrackerSiFZ)
1248         meTrackerSiFZ->Fill(bSurface.toGlobal(itHit->localPosition()).z());
1249       if (meTrackerSiPhi)
1250         meTrackerSiPhi->Fill(bSurface.toGlobal(itHit->localPosition()).phi());
1251       if (meTrackerSiEta)
1252         meTrackerSiEta->Fill(bSurface.toGlobal(itHit->localPosition()).eta());
1253 
1254     } else {
1255       edm::LogWarning(MsgLoggerCat) << "SiFwd PSimHit " << i << " is expected to be (det,subdet) = (" << dTrk << ","
1256                                     << sdSiTOB << " || " << sdSiTEC << "); value returned is: (" << detector << ","
1257                                     << subdetector << ")";
1258       continue;
1259     }  // end check detector type
1260   }    // end loop through SiFwd Hits
1261 
1262   if (verbosity > 1) {
1263     eventout += "\n          Number of Silicon Forward Hits collected:.. ";
1264     eventout += j;
1265   }
1266 
1267   nSiHits += j;
1268 
1269   if (meTrackerSi[0])
1270     meTrackerSi[0]->Fill((float)nSiHits);
1271   if (meTrackerSi[1])
1272     meTrackerSi[1]->Fill((float)nSiHits);
1273 
1274   if (verbosity > 0)
1275     edm::LogInfo(MsgLoggerCat) << eventout << "\n";
1276 
1277   return;
1278 }
1279 
1280 void GlobalHitsAnalyzer::fillMuon(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
1281   nMuonHits = 0;
1282   std::string MsgLoggerCat = "GlobalHitsAnalyzer_fillMuon";
1283 
1284   TString eventout;
1285   if (verbosity > 0)
1286     eventout = "\nGathering info:";
1287 
1288   // iterator to access containers
1289   edm::PSimHitContainer::const_iterator itHit;
1290 
1291   ///////////////////////
1292   // access the CSC Muon
1293   ///////////////////////
1294   // access the CSC Muon geometry
1295   const auto &theCSCGeometry = iSetup.getHandle(cscGeomToken_);
1296   if (!theCSCGeometry.isValid()) {
1297     edm::LogWarning(MsgLoggerCat) << "Unable to find MuonGeometryRecord for the CSCGeometry in event!";
1298     return;
1299   }
1300   const CSCGeometry &theCSCMuon(*theCSCGeometry);
1301 
1302   // get Muon CSC information
1303   edm::Handle<edm::PSimHitContainer> MuonCSCContainer;
1304   iEvent.getByToken(MuonCscSrc_Token_, MuonCSCContainer);
1305   if (!MuonCSCContainer.isValid()) {
1306     LogDebug(MsgLoggerCat) << "Unable to find MuonCSCHits in event!";
1307     validMuonCSC = false;
1308   }
1309 
1310   if (validMuonCSC) {
1311     // cycle through container
1312     int i = 0, j = 0;
1313     for (itHit = MuonCSCContainer->begin(); itHit != MuonCSCContainer->end(); ++itHit) {
1314       ++i;
1315 
1316       // create a DetId from the detUnitId
1317       DetId theDetUnitId(itHit->detUnitId());
1318       int detector = theDetUnitId.det();
1319       int subdetector = theDetUnitId.subdetId();
1320 
1321       // check that expected detector is returned
1322       if ((detector == dMuon) && (subdetector == sdMuonCSC)) {
1323         // get the GeomDetUnit from the geometry using theDetUnitID
1324         const GeomDetUnit *theDet = theCSCMuon.idToDetUnit(theDetUnitId);
1325 
1326         if (!theDet) {
1327           edm::LogWarning(MsgLoggerCat) << "Unable to get GeomDetUnit from theCSCMuon for hit " << i;
1328           continue;
1329         }
1330 
1331         ++j;
1332 
1333         // get the Surface of the hit (knows how to go from local <-> global)
1334         const BoundPlane &bSurface = theDet->surface();
1335 
1336         if (meMuonCscToF[0])
1337           meMuonCscToF[0]->Fill(itHit->tof());
1338         if (meMuonCscToF[1])
1339           meMuonCscToF[1]->Fill(itHit->tof());
1340         if (meMuonCscZ)
1341           meMuonCscZ->Fill(bSurface.toGlobal(itHit->localPosition()).z());
1342         if (meMuonPhi)
1343           meMuonPhi->Fill(bSurface.toGlobal(itHit->localPosition()).phi());
1344         if (meMuonEta)
1345           meMuonEta->Fill(bSurface.toGlobal(itHit->localPosition()).eta());
1346 
1347       } else {
1348         edm::LogWarning(MsgLoggerCat) << "MuonCsc PSimHit " << i << " is expected to be (det,subdet) = (" << dMuon
1349                                       << "," << sdMuonCSC << "); value returned is: (" << detector << "," << subdetector
1350                                       << ")";
1351         continue;
1352       }  // end detector type check
1353     }    // end loop through CSC Hits
1354 
1355     if (verbosity > 1) {
1356       eventout += "\n          Number of CSC muon Hits collected:......... ";
1357       eventout += j;
1358     }
1359 
1360     nMuonHits += j;
1361   }
1362 
1363   /////////////////////
1364   // access the DT Muon
1365   /////////////////////
1366   // access the DT Muon geometry
1367   const auto &theDTGeometry = iSetup.getHandle(dtGeomToken_);
1368   if (!theDTGeometry.isValid()) {
1369     edm::LogWarning(MsgLoggerCat) << "Unable to find MuonGeometryRecord for the DTGeometry in event!";
1370     return;
1371   }
1372   const DTGeometry &theDTMuon(*theDTGeometry);
1373 
1374   // get Muon DT information
1375   edm::Handle<edm::PSimHitContainer> MuonDtContainer;
1376   iEvent.getByToken(MuonDtSrc_Token_, MuonDtContainer);
1377   if (!MuonDtContainer.isValid()) {
1378     LogDebug(MsgLoggerCat) << "Unable to find MuonDTHits in event!";
1379     validMuonDt = false;
1380   }
1381 
1382   if (validMuonDt) {
1383     // cycle through container
1384     int i = 0, j = 0;
1385     for (itHit = MuonDtContainer->begin(); itHit != MuonDtContainer->end(); ++itHit) {
1386       ++i;
1387 
1388       // create a DetId from the detUnitId
1389       DetId theDetUnitId(itHit->detUnitId());
1390       int detector = theDetUnitId.det();
1391       int subdetector = theDetUnitId.subdetId();
1392 
1393       // check that expected detector is returned
1394       if ((detector == dMuon) && (subdetector == sdMuonDT)) {
1395         // CSC uses wires and layers rather than the full detID
1396         // get the wireId
1397         DTWireId wireId(itHit->detUnitId());
1398 
1399         // get the DTLayer from the geometry using the wireID
1400         const DTLayer *theDet = theDTMuon.layer(wireId.layerId());
1401 
1402         if (!theDet) {
1403           edm::LogWarning(MsgLoggerCat) << "Unable to get GeomDetUnit from theDtMuon for hit " << i;
1404           continue;
1405         }
1406 
1407         ++j;
1408 
1409         // get the Surface of the hit (knows how to go from local <-> global)
1410         const BoundPlane &bSurface = theDet->surface();
1411 
1412         if (meMuonDtToF[0])
1413           meMuonDtToF[0]->Fill(itHit->tof());
1414         if (meMuonDtToF[1])
1415           meMuonDtToF[1]->Fill(itHit->tof());
1416         if (meMuonDtR)
1417           meMuonDtR->Fill(bSurface.toGlobal(itHit->localPosition()).perp());
1418         if (meMuonPhi)
1419           meMuonPhi->Fill(bSurface.toGlobal(itHit->localPosition()).phi());
1420         if (meMuonEta)
1421           meMuonEta->Fill(bSurface.toGlobal(itHit->localPosition()).eta());
1422 
1423       } else {
1424         edm::LogWarning(MsgLoggerCat) << "MuonDt PSimHit " << i << " is expected to be (det,subdet) = (" << dMuon << ","
1425                                       << sdMuonDT << "); value returned is: (" << detector << "," << subdetector << ")";
1426         continue;
1427       }  // end detector type check
1428     }    // end loop through DT Hits
1429 
1430     if (verbosity > 1) {
1431       eventout += "\n          Number of DT muon Hits collected:.......... ";
1432       eventout += j;
1433     }
1434 
1435     nMuonHits += j;
1436   }
1437 
1438   // int RPCBrl = 0, RPCFwd = 0;
1439   ///////////////////////
1440   // access the RPC Muon
1441   ///////////////////////
1442   // access the RPC Muon geometry
1443   const auto &theRPCGeometry = iSetup.getHandle(rpcGeomToken_);
1444   if (!theRPCGeometry.isValid()) {
1445     edm::LogWarning(MsgLoggerCat) << "Unable to find MuonGeometryRecord for the RPCGeometry in event!";
1446     return;
1447   }
1448   const RPCGeometry &theRPCMuon(*theRPCGeometry);
1449 
1450   // get Muon RPC information
1451   edm::Handle<edm::PSimHitContainer> MuonRPCContainer;
1452   iEvent.getByToken(MuonRpcSrc_Token_, MuonRPCContainer);
1453   if (!MuonRPCContainer.isValid()) {
1454     LogDebug(MsgLoggerCat) << "Unable to find MuonRPCHits in event!";
1455     validMuonRPC = false;
1456   }
1457 
1458   if (validMuonRPC) {
1459     // cycle through container
1460     int i = 0, j = 0;
1461     int RPCBrl = 0, RPCFwd = 0;
1462     for (itHit = MuonRPCContainer->begin(); itHit != MuonRPCContainer->end(); ++itHit) {
1463       ++i;
1464 
1465       // create a DetID from the detUnitId
1466       DetId theDetUnitId(itHit->detUnitId());
1467       int detector = theDetUnitId.det();
1468       int subdetector = theDetUnitId.subdetId();
1469 
1470       // check that expected detector is returned
1471       if ((detector == dMuon) && (subdetector == sdMuonRPC)) {
1472         // get an RPCDetID from the detUnitID
1473         RPCDetId RPCId(itHit->detUnitId());
1474 
1475         // find the region of the RPC hit
1476         int region = RPCId.region();
1477 
1478         // get the GeomDetUnit from the geometry using the RPCDetId
1479         const GeomDetUnit *theDet = theRPCMuon.idToDetUnit(theDetUnitId);
1480 
1481         if (!theDet) {
1482           edm::LogWarning(MsgLoggerCat) << "Unable to get GeomDetUnit from theRPCMuon for hit " << i;
1483           continue;
1484         }
1485 
1486         ++j;
1487 
1488         // get the Surface of the hit (knows how to go from local <-> global)
1489         const BoundPlane &bSurface = theDet->surface();
1490 
1491         // gather necessary information
1492         if ((region == sdMuonRPCRgnFwdp) || (region == sdMuonRPCRgnFwdn)) {
1493           ++RPCFwd;
1494 
1495           if (meMuonRpcFToF[0])
1496             meMuonRpcFToF[0]->Fill(itHit->tof());
1497           if (meMuonRpcFToF[1])
1498             meMuonRpcFToF[1]->Fill(itHit->tof());
1499           if (meMuonRpcFZ)
1500             meMuonRpcFZ->Fill(bSurface.toGlobal(itHit->localPosition()).z());
1501           if (meMuonPhi)
1502             meMuonPhi->Fill(bSurface.toGlobal(itHit->localPosition()).phi());
1503           if (meMuonEta)
1504             meMuonEta->Fill(bSurface.toGlobal(itHit->localPosition()).eta());
1505 
1506         } else if (region == sdMuonRPCRgnBrl) {
1507           ++RPCBrl;
1508 
1509           if (meMuonRpcBToF[0])
1510             meMuonRpcBToF[0]->Fill(itHit->tof());
1511           if (meMuonRpcBToF[1])
1512             meMuonRpcBToF[1]->Fill(itHit->tof());
1513           if (meMuonRpcBR)
1514             meMuonRpcBR->Fill(bSurface.toGlobal(itHit->localPosition()).perp());
1515           if (meMuonPhi)
1516             meMuonPhi->Fill(bSurface.toGlobal(itHit->localPosition()).phi());
1517           if (meMuonEta)
1518             meMuonEta->Fill(bSurface.toGlobal(itHit->localPosition()).eta());
1519 
1520         } else {
1521           edm::LogWarning(MsgLoggerCat) << "Invalid region for RPC Muon hit" << i;
1522           continue;
1523         }  // end check of region
1524       } else {
1525         edm::LogWarning(MsgLoggerCat) << "MuonRpc PSimHit " << i << " is expected to be (det,subdet) = (" << dMuon
1526                                       << "," << sdMuonRPC << "); value returned is: (" << detector << "," << subdetector
1527                                       << ")";
1528         continue;
1529       }  // end detector type check
1530     }    // end loop through RPC Hits
1531 
1532     if (verbosity > 1) {
1533       eventout += "\n          Number of RPC muon Hits collected:......... ";
1534       eventout += j;
1535       eventout += "\n                    RPC Barrel muon Hits:............ ";
1536       eventout += RPCBrl;
1537       eventout += "\n                    RPC Forward muon Hits:........... ";
1538       eventout += RPCFwd;
1539     }
1540 
1541     nMuonHits += j;
1542   }
1543 
1544   if (meMuon[0])
1545     meMuon[0]->Fill((float)nMuonHits);
1546   if (meMuon[1])
1547     meMuon[1]->Fill((float)nMuonHits);
1548 
1549   if (verbosity > 0)
1550     edm::LogInfo(MsgLoggerCat) << eventout << "\n";
1551 
1552   return;
1553 }
1554 
1555 void GlobalHitsAnalyzer::fillECal(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
1556   std::string MsgLoggerCat = "GlobalHitsAnalyzer_fillECal";
1557 
1558   TString eventout;
1559   if (verbosity > 0)
1560     eventout = "\nGathering info:";
1561 
1562   // access the calorimeter geometry
1563   const auto &theCaloGeometry = iSetup.getHandle(caloGeomToken_);
1564   if (!theCaloGeometry.isValid()) {
1565     edm::LogWarning(MsgLoggerCat) << "Unable to find CaloGeometryRecord in event!";
1566     return;
1567   }
1568   const CaloGeometry &theCalo(*theCaloGeometry);
1569 
1570   // iterator to access containers
1571   edm::PCaloHitContainer::const_iterator itHit;
1572 
1573   ///////////////////////////////
1574   // get  ECal information
1575   ///////////////////////////////
1576   edm::PCaloHitContainer theECalHits;
1577   // extract EB container
1578   edm::Handle<edm::PCaloHitContainer> EBContainer;
1579   iEvent.getByToken(ECalEBSrc_Token_, EBContainer);
1580   if (!EBContainer.isValid()) {
1581     LogDebug(MsgLoggerCat) << "Unable to find EcalHitsEB in event!";
1582     validEB = false;
1583   }
1584   // extract EE container
1585   edm::Handle<edm::PCaloHitContainer> EEContainer;
1586   iEvent.getByToken(ECalEESrc_Token_, EEContainer);
1587   if (!EEContainer.isValid()) {
1588     LogDebug(MsgLoggerCat) << "Unable to find EcalHitsEE in event!";
1589     validEE = false;
1590   }
1591   // place both containers into new container
1592   if (validEB)
1593     theECalHits.insert(theECalHits.end(), EBContainer->begin(), EBContainer->end());
1594   if (validEE)
1595     theECalHits.insert(theECalHits.end(), EEContainer->begin(), EEContainer->end());
1596 
1597   // cycle through new container
1598   int i = 0, j = 0;
1599   for (itHit = theECalHits.begin(); itHit != theECalHits.end(); ++itHit) {
1600     ++i;
1601 
1602     // create a DetId from the detUnitId
1603     DetId theDetUnitId(itHit->id());
1604     int detector = theDetUnitId.det();
1605     int subdetector = theDetUnitId.subdetId();
1606 
1607     // check that expected detector is returned
1608     if ((detector == dEcal) && ((subdetector == sdEcalBrl) || (subdetector == sdEcalFwd))) {
1609       // get the Cell geometry
1610       auto theDet = (theCalo.getSubdetectorGeometry(theDetUnitId))->getGeometry(theDetUnitId);
1611 
1612       if (!theDet) {
1613         edm::LogWarning(MsgLoggerCat) << "Unable to get CaloCellGeometry from ECalHits for Hit " << i;
1614         continue;
1615       }
1616 
1617       ++j;
1618 
1619       // get the global position of the cell
1620       const GlobalPoint &globalposition = theDet->getPosition();
1621 
1622       if (meCaloEcalE[0])
1623         meCaloEcalE[0]->Fill(itHit->energy());
1624       if (meCaloEcalE[1])
1625         meCaloEcalE[1]->Fill(itHit->energy());
1626       if (meCaloEcalToF[0])
1627         meCaloEcalToF[0]->Fill(itHit->time());
1628       if (meCaloEcalToF[1])
1629         meCaloEcalToF[1]->Fill(itHit->time());
1630       if (meCaloEcalPhi)
1631         meCaloEcalPhi->Fill(globalposition.phi());
1632       if (meCaloEcalEta)
1633         meCaloEcalEta->Fill(globalposition.eta());
1634 
1635     } else {
1636       edm::LogWarning(MsgLoggerCat) << "ECal PCaloHit " << i << " is expected to be (det,subdet) = (" << dEcal << ","
1637                                     << sdEcalBrl << " || " << sdEcalFwd << "); value returned is: (" << detector << ","
1638                                     << subdetector << ")";
1639       continue;
1640     }  // end detector type check
1641   }    // end loop through ECal Hits
1642 
1643   if (verbosity > 1) {
1644     eventout += "\n          Number of ECal Hits collected:............. ";
1645     eventout += j;
1646   }
1647 
1648   if (meCaloEcal[0])
1649     meCaloEcal[0]->Fill((float)j);
1650   if (meCaloEcal[1])
1651     meCaloEcal[1]->Fill((float)j);
1652 
1653   ////////////////////////////
1654   // Get Preshower information
1655   ////////////////////////////
1656   // extract PreShower container
1657   edm::Handle<edm::PCaloHitContainer> PreShContainer;
1658   iEvent.getByToken(ECalESSrc_Token_, PreShContainer);
1659   if (!PreShContainer.isValid()) {
1660     LogDebug(MsgLoggerCat) << "Unable to find EcalHitsES in event!";
1661     validPresh = false;
1662   }
1663 
1664   if (validPresh) {
1665     // cycle through container
1666     int i = 0, j = 0;
1667     for (itHit = PreShContainer->begin(); itHit != PreShContainer->end(); ++itHit) {
1668       ++i;
1669 
1670       // create a DetId from the detUnitId
1671       DetId theDetUnitId(itHit->id());
1672       int detector = theDetUnitId.det();
1673       int subdetector = theDetUnitId.subdetId();
1674 
1675       // check that expected detector is returned
1676       if ((detector == dEcal) && (subdetector == sdEcalPS)) {
1677         // get the Cell geometry
1678         auto theDet = (theCalo.getSubdetectorGeometry(theDetUnitId))->getGeometry(theDetUnitId);
1679 
1680         if (!theDet) {
1681           edm::LogWarning(MsgLoggerCat) << "Unable to get CaloCellGeometry from PreShContainer for Hit " << i;
1682           continue;
1683         }
1684 
1685         ++j;
1686 
1687         // get the global position of the cell
1688         const GlobalPoint &globalposition = theDet->getPosition();
1689 
1690         if (meCaloPreShE[0])
1691           meCaloPreShE[0]->Fill(itHit->energy());
1692         if (meCaloPreShE[1])
1693           meCaloPreShE[1]->Fill(itHit->energy());
1694         if (meCaloPreShToF[0])
1695           meCaloPreShToF[0]->Fill(itHit->time());
1696         if (meCaloPreShToF[1])
1697           meCaloPreShToF[1]->Fill(itHit->time());
1698         if (meCaloPreShPhi)
1699           meCaloPreShPhi->Fill(globalposition.phi());
1700         if (meCaloPreShEta)
1701           meCaloPreShEta->Fill(globalposition.eta());
1702 
1703       } else {
1704         edm::LogWarning(MsgLoggerCat) << "PreSh PCaloHit " << i << " is expected to be (det,subdet) = (" << dEcal << ","
1705                                       << sdEcalPS << "); value returned is: (" << detector << "," << subdetector << ")";
1706         continue;
1707       }  // end detector type check
1708     }    // end loop through PreShower Hits
1709 
1710     if (verbosity > 1) {
1711       eventout += "\n          Number of PreSh Hits collected:............ ";
1712       eventout += j;
1713     }
1714 
1715     if (meCaloPreSh[0])
1716       meCaloPreSh[0]->Fill((float)j);
1717     if (meCaloPreSh[1])
1718       meCaloPreSh[1]->Fill((float)j);
1719   }
1720 
1721   if (verbosity > 0)
1722     edm::LogInfo(MsgLoggerCat) << eventout << "\n";
1723 
1724   return;
1725 }
1726 
1727 void GlobalHitsAnalyzer::fillHCal(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
1728   std::string MsgLoggerCat = "GlobalHitsAnalyzer_fillHCal";
1729 
1730   TString eventout;
1731   if (verbosity > 0)
1732     eventout = "\nGathering info:";
1733 
1734   // access the calorimeter geometry
1735   const auto &theCaloGeometry = iSetup.getHandle(caloGeomToken_);
1736   if (!theCaloGeometry.isValid()) {
1737     edm::LogWarning(MsgLoggerCat) << "Unable to find CaloGeometryRecord in event!";
1738     return;
1739   }
1740   const CaloGeometry &theCalo(*theCaloGeometry);
1741 
1742   const HcalDDDRecConstants *hcons = &iSetup.getData(hcaldddRecToken_);
1743 
1744   // iterator to access containers
1745   edm::PCaloHitContainer::const_iterator itHit;
1746 
1747   ///////////////////////////////
1748   // get  HCal information
1749   ///////////////////////////////
1750   // extract HCal container
1751   edm::Handle<edm::PCaloHitContainer> HCalContainer;
1752   iEvent.getByToken(HCalSrc_Token_, HCalContainer);
1753   if (!HCalContainer.isValid()) {
1754     LogDebug(MsgLoggerCat) << "Unable to find HCalHits in event!";
1755     validHcal = false;
1756   }
1757 
1758   if (validHcal) {
1759     // cycle through container
1760     int i = 0, j = 0;
1761     for (itHit = HCalContainer->begin(); itHit != HCalContainer->end(); ++itHit) {
1762       ++i;
1763 
1764       // create a DetId from the detUnitId
1765       DetId theDetUnitId;
1766       unsigned int id_ = itHit->id();
1767       if (testNumber)
1768         theDetUnitId = HcalHitRelabeller::relabel(id_, hcons);
1769       else
1770         theDetUnitId = id_;
1771 
1772       int detector = theDetUnitId.det();
1773       int subdetector = theDetUnitId.subdetId();
1774 
1775       // check that expected detector is returned
1776       if ((detector == dHcal) && ((subdetector == sdHcalBrl) || (subdetector == sdHcalEC) ||
1777                                   (subdetector == sdHcalOut) || (subdetector == sdHcalFwd))) {
1778         // get the Cell geometry
1779         const HcalGeometry *theDet = dynamic_cast<const HcalGeometry *>(theCalo.getSubdetectorGeometry(theDetUnitId));
1780 
1781         if (!theDet) {
1782           edm::LogWarning(MsgLoggerCat) << "Unable to get HcalGeometry from HCalContainer for Hit " << i;
1783           continue;
1784         }
1785 
1786         ++j;
1787 
1788         // get the global position of the cell
1789         const GlobalPoint &globalposition = theDet->getPosition(theDetUnitId);
1790 
1791         if (meCaloHcalE[0])
1792           meCaloHcalE[0]->Fill(itHit->energy());
1793         if (meCaloHcalE[1])
1794           meCaloHcalE[1]->Fill(itHit->energy());
1795         if (meCaloHcalToF[0])
1796           meCaloHcalToF[0]->Fill(itHit->time());
1797         if (meCaloHcalToF[1])
1798           meCaloHcalToF[1]->Fill(itHit->time());
1799         if (meCaloHcalPhi)
1800           meCaloHcalPhi->Fill(globalposition.phi());
1801         if (meCaloHcalEta)
1802           meCaloHcalEta->Fill(globalposition.eta());
1803 
1804       } else {
1805         edm::LogWarning(MsgLoggerCat) << "HCal PCaloHit " << i << " is expected to be (det,subdet) = (" << dHcal << ","
1806                                       << sdHcalBrl << " || " << sdHcalEC << " || " << sdHcalOut << " || " << sdHcalFwd
1807                                       << "); value returned is: (" << detector << "," << subdetector << ")";
1808         continue;
1809       }  // end detector type check
1810     }    // end loop through HCal Hits
1811 
1812     if (verbosity > 1) {
1813       eventout += "\n          Number of HCal Hits collected:............. ";
1814       eventout += j;
1815     }
1816 
1817     if (meCaloHcal[0])
1818       meCaloHcal[0]->Fill((float)j);
1819     if (meCaloHcal[1])
1820       meCaloHcal[1]->Fill((float)j);
1821   }
1822 
1823   if (verbosity > 0)
1824     edm::LogInfo(MsgLoggerCat) << eventout << "\n";
1825 
1826   return;
1827 }