Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-07-13 03:16:18

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