Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-02-14 23:39:57 

0001 /** \file GlobalRecHitsAnalyzer.cc
0002  *  
0003  *  See header file for description of class
0004  *
0005  *  \author M. Strang SUNY-Buffalo
0006  *  Testing by Ken Smith
0007  */
0008 #include "Validation/GlobalRecHits/interface/GlobalRecHitsAnalyzer.h"
0009 #include "DQMServices/Core/interface/DQMStore.h"
0010 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0011 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0012 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0013 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0014 using namespace std;
0015 
0016 GlobalRecHitsAnalyzer::GlobalRecHitsAnalyzer(const edm::ParameterSet& iPSet)
0017     : fName(""),
0018       verbosity(0),
0019       frequency(0),
0020       label(""),
0021       getAllProvenances(false),
0022       printProvenanceInfo(false),
0023       trackerHitAssociatorConfig_(iPSet, consumesCollector()),
0024       caloGeomToken_(esConsumes()),
0025       tTopoToken_(esConsumes()),
0026       tGeomToken_(esConsumes()),
0027       dtGeomToken_(esConsumes()),
0028       cscGeomToken_(esConsumes()),
0029       rpcGeomToken_(esConsumes()),
0030       count(0) {
0031   consumesMany<edm::SortedCollection<HBHERecHit, edm::StrictWeakOrdering<HBHERecHit>>>();
0032   consumesMany<edm::SortedCollection<HFRecHit, edm::StrictWeakOrdering<HFRecHit>>>();
0033   consumesMany<edm::SortedCollection<HORecHit, edm::StrictWeakOrdering<HORecHit>>>();
0034   std::string MsgLoggerCat = "GlobalRecHitsAnalyzer_GlobalRecHitsAnalyzer";
0035 
0036   // get information from parameter set
0037   fName = iPSet.getUntrackedParameter<std::string>("Name");
0038   verbosity = iPSet.getUntrackedParameter<int>("Verbosity");
0039   frequency = iPSet.getUntrackedParameter<int>("Frequency");
0040   edm::ParameterSet m_Prov = iPSet.getParameter<edm::ParameterSet>("ProvenanceLookup");
0041   getAllProvenances = m_Prov.getUntrackedParameter<bool>("GetAllProvenances");
0042   printProvenanceInfo = m_Prov.getUntrackedParameter<bool>("PrintProvenanceInfo");
0043   hitsProducer = iPSet.getParameter<std::string>("hitsProducer");
0044 
0045   //get Labels to use to extract information
0046   ECalEBSrc_ = iPSet.getParameter<edm::InputTag>("ECalEBSrc");
0047   ECalUncalEBSrc_ = iPSet.getParameter<edm::InputTag>("ECalUncalEBSrc");
0048   ECalEESrc_ = iPSet.getParameter<edm::InputTag>("ECalEESrc");
0049   ECalUncalEESrc_ = iPSet.getParameter<edm::InputTag>("ECalUncalEESrc");
0050   ECalESSrc_ = iPSet.getParameter<edm::InputTag>("ECalESSrc");
0051   HCalSrc_ = iPSet.getParameter<edm::InputTag>("HCalSrc");
0052   SiStripSrc_ = iPSet.getParameter<edm::InputTag>("SiStripSrc");
0053   SiPxlSrc_ = iPSet.getParameter<edm::InputTag>("SiPxlSrc");
0054   MuDTSrc_ = iPSet.getParameter<edm::InputTag>("MuDTSrc");
0055   MuDTSimSrc_ = iPSet.getParameter<edm::InputTag>("MuDTSimSrc");
0056   MuCSCSrc_ = iPSet.getParameter<edm::InputTag>("MuCSCSrc");
0057   MuRPCSrc_ = iPSet.getParameter<edm::InputTag>("MuRPCSrc");
0058   MuRPCSimSrc_ = iPSet.getParameter<edm::InputTag>("MuRPCSimSrc");
0059 
0060   // fix for consumes
0061   ECalUncalEBSrc_Token_ = consumes<EBUncalibratedRecHitCollection>(iPSet.getParameter<edm::InputTag>("ECalUncalEBSrc"));
0062   ECalUncalEESrc_Token_ = consumes<EEUncalibratedRecHitCollection>(iPSet.getParameter<edm::InputTag>("ECalUncalEESrc"));
0063   ECalEBSrc_Token_ = consumes<EBRecHitCollection>(iPSet.getParameter<edm::InputTag>("ECalEBSrc"));
0064   ECalEESrc_Token_ = consumes<EERecHitCollection>(iPSet.getParameter<edm::InputTag>("ECalEESrc"));
0065   ECalESSrc_Token_ = consumes<ESRecHitCollection>(iPSet.getParameter<edm::InputTag>("ECalESSrc"));
0066   HCalSrc_Token_ = consumes<edm::PCaloHitContainer>(iPSet.getParameter<edm::InputTag>("HCalSrc"));
0067   SiStripSrc_Token_ = consumes<SiStripMatchedRecHit2DCollection>(iPSet.getParameter<edm::InputTag>("SiStripSrc"));
0068   SiPxlSrc_Token_ = consumes<SiPixelRecHitCollection>(iPSet.getParameter<edm::InputTag>("SiPxlSrc"));
0069 
0070   MuDTSrc_Token_ = consumes<DTRecHitCollection>(iPSet.getParameter<edm::InputTag>("MuDTSrc"));
0071   MuDTSimSrc_Token_ = consumes<edm::PSimHitContainer>(iPSet.getParameter<edm::InputTag>("MuDTSimSrc"));
0072 
0073   MuCSCSrc_Token_ = consumes<CSCRecHit2DCollection>(iPSet.getParameter<edm::InputTag>("MuCSCSrc"));
0074   MuCSCHits_Token_ = consumes<CrossingFrame<PSimHit>>(
0075       edm::InputTag(std::string("mix"), iPSet.getParameter<std::string>("hitsProducer") + std::string("MuonCSCHits")));
0076 
0077   MuRPCSrc_Token_ = consumes<RPCRecHitCollection>(iPSet.getParameter<edm::InputTag>("MuRPCSrc"));
0078   MuRPCSimSrc_Token_ = consumes<edm::PSimHitContainer>(iPSet.getParameter<edm::InputTag>("MuRPCSimSrc"));
0079 
0080   EBHits_Token_ = consumes<CrossingFrame<PCaloHit>>(
0081       edm::InputTag(std::string("mix"), iPSet.getParameter<std::string>("hitsProducer") + std::string("EcalHitsEB")));
0082   EEHits_Token_ = consumes<CrossingFrame<PCaloHit>>(
0083       edm::InputTag(std::string("mix"), iPSet.getParameter<std::string>("hitsProducer") + std::string("EcalHitsEE")));
0084   ESHits_Token_ = consumes<CrossingFrame<PCaloHit>>(
0085       edm::InputTag(std::string("mix"), iPSet.getParameter<std::string>("hitsProducer") + std::string("EcalHitsES")));
0086 
0087   // use value of first digit to determine default output level (inclusive)
0088   // 0 is none, 1 is basic, 2 is fill output, 3 is gather output
0089   verbosity %= 10;
0090 
0091   // create persistent object
0092   // produces<PGlobalRecHit>(label);
0093 
0094   // print out Parameter Set information being used
0095   if (verbosity >= 0) {
0096     edm::LogInfo(MsgLoggerCat) << "\n===============================\n"
0097                                << "Initialized as EDProducer with parameter values:\n"
0098                                << "    Name           = " << fName << "\n"
0099                                << "    Verbosity      = " << verbosity << "\n"
0100                                << "    Frequency      = " << frequency << "\n"
0101                                << "    GetProv        = " << getAllProvenances << "\n"
0102                                << "    PrintProv      = " << printProvenanceInfo << "\n"
0103                                << "    ECalEBSrc      = " << ECalEBSrc_.label() << ":" << ECalEBSrc_.instance() << "\n"
0104                                << "    ECalUncalEBSrc = " << ECalUncalEBSrc_.label() << ":"
0105                                << ECalUncalEBSrc_.instance() << "\n"
0106                                << "    ECalEESrc      = " << ECalEESrc_.label() << ":" << ECalUncalEESrc_.instance()
0107                                << "\n"
0108                                << "    ECalUncalEESrc = " << ECalUncalEESrc_.label() << ":" << ECalEESrc_.instance()
0109                                << "\n"
0110                                << "    ECalESSrc      = " << ECalESSrc_.label() << ":" << ECalESSrc_.instance() << "\n"
0111                                << "    HCalSrc        = " << HCalSrc_.label() << ":" << HCalSrc_.instance() << "\n"
0112                                << "    SiStripSrc     = " << SiStripSrc_.label() << ":" << SiStripSrc_.instance()
0113                                << "\n"
0114                                << "    SiPixelSrc     = " << SiPxlSrc_.label() << ":" << SiPxlSrc_.instance() << "\n"
0115                                << "    MuDTSrc        = " << MuDTSrc_.label() << ":" << MuDTSrc_.instance() << "\n"
0116                                << "    MuDTSimSrc     = " << MuDTSimSrc_.label() << ":" << MuDTSimSrc_.instance()
0117                                << "\n"
0118                                << "    MuCSCSrc       = " << MuCSCSrc_.label() << ":" << MuCSCSrc_.instance() << "\n"
0119                                << "    MuRPCSrc       = " << MuRPCSrc_.label() << ":" << MuRPCSrc_.instance() << "\n"
0120                                << "    MuRPCSimSrc    = " << MuRPCSimSrc_.label() << ":" << MuRPCSimSrc_.instance()
0121                                << "\n"
0122                                << "===============================\n";
0123   }
0124 }
0125 
0126 GlobalRecHitsAnalyzer::~GlobalRecHitsAnalyzer() {}
0127 
0128 void GlobalRecHitsAnalyzer::bookHistograms(DQMStore::IBooker& iBooker, edm::Run const&, edm::EventSetup const&) {
0129   // Si Strip
0130   string SiStripString[19] = {"TECW1",
0131                               "TECW2",
0132                               "TECW3",
0133                               "TECW4",
0134                               "TECW5",
0135                               "TECW6",
0136                               "TECW7",
0137                               "TECW8",
0138                               "TIBL1",
0139                               "TIBL2",
0140                               "TIBL3",
0141                               "TIBL4",
0142                               "TIDW1",
0143                               "TIDW2",
0144                               "TIDW3",
0145                               "TOBL1",
0146                               "TOBL2",
0147                               "TOBL3",
0148                               "TOBL4"};
0149   for (int i = 0; i < 19; ++i) {
0150     mehSiStripn[i] = nullptr;
0151     mehSiStripResX[i] = nullptr;
0152     mehSiStripResY[i] = nullptr;
0153   }
0154   string hcharname, hchartitle;
0155   iBooker.setCurrentFolder("GlobalRecHitsV/SiStrips");
0156   for (int amend = 0; amend < 19; ++amend) {
0157     hcharname = "hSiStripn_" + SiStripString[amend];
0158     hchartitle = SiStripString[amend] + "  rechits";
0159     mehSiStripn[amend] = iBooker.book1D(hcharname, hchartitle, 200, 0., 200.);
0160     mehSiStripn[amend]->setAxisTitle("Number of hits in " + SiStripString[amend], 1);
0161     mehSiStripn[amend]->setAxisTitle("Count", 2);
0162     hcharname = "hSiStripResX_" + SiStripString[amend];
0163     hchartitle = SiStripString[amend] + " rechit x resolution";
0164     mehSiStripResX[amend] = iBooker.book1D(hcharname, hchartitle, 200, -0.02, .02);
0165     mehSiStripResX[amend]->setAxisTitle("X-resolution in " + SiStripString[amend], 1);
0166     mehSiStripResX[amend]->setAxisTitle("Count", 2);
0167     hcharname = "hSiStripResY_" + SiStripString[amend];
0168     hchartitle = SiStripString[amend] + " rechit y resolution";
0169     mehSiStripResY[amend] = iBooker.book1D(hcharname, hchartitle, 200, -0.02, .02);
0170     mehSiStripResY[amend]->setAxisTitle("Y-resolution in " + SiStripString[amend], 1);
0171     mehSiStripResY[amend]->setAxisTitle("Count", 2);
0172   }
0173 
0174   //HCal
0175   //string hcharname, hchartitle;
0176   string HCalString[4] = {"HB", "HE", "HF", "HO"};
0177   float HCalnUpper[4] = {3000., 3000., 3000., 3000.};
0178   float HCalnLower[4] = {0., 0., 0., 0.};
0179   for (int j = 0; j < 4; ++j) {
0180     mehHcaln[j] = nullptr;
0181     mehHcalRes[j] = nullptr;
0182   }
0183 
0184   iBooker.setCurrentFolder("GlobalRecHitsV/HCals");
0185   for (int amend = 0; amend < 4; ++amend) {
0186     hcharname = "hHcaln_" + HCalString[amend];
0187     hchartitle = HCalString[amend] + "  rechits";
0188     mehHcaln[amend] = iBooker.book1D(hcharname, hchartitle, 1000, HCalnLower[amend], HCalnUpper[amend]);
0189     mehHcaln[amend]->setAxisTitle("Number of RecHits", 1);
0190     mehHcaln[amend]->setAxisTitle("Count", 2);
0191     hcharname = "hHcalRes_" + HCalString[amend];
0192     hchartitle = HCalString[amend] + "  rechit resolution";
0193     mehHcalRes[amend] = iBooker.book1D(hcharname, hchartitle, 25, -2., 2.);
0194     mehHcalRes[amend]->setAxisTitle("RecHit E - SimHit E", 1);
0195     mehHcalRes[amend]->setAxisTitle("Count", 2);
0196   }
0197 
0198   //Ecal
0199   string ECalString[3] = {"EB", "EE", "ES"};
0200   int ECalnBins[3] = {1000, 3000, 150};
0201   float ECalnUpper[3] = {20000., 62000., 3000.};
0202   float ECalnLower[3] = {0., 0., 0.};
0203   int ECalResBins[3] = {200, 200, 200};
0204   float ECalResUpper[3] = {1., 0.3, .0002};
0205   float ECalResLower[3] = {-1., -0.3, -.0002};
0206   for (int i = 0; i < 3; ++i) {
0207     mehEcaln[i] = nullptr;
0208     mehEcalRes[i] = nullptr;
0209   }
0210   iBooker.setCurrentFolder("GlobalRecHitsV/ECals");
0211 
0212   for (int amend = 0; amend < 3; ++amend) {
0213     hcharname = "hEcaln_" + ECalString[amend];
0214     hchartitle = ECalString[amend] + "  rechits";
0215     mehEcaln[amend] = iBooker.book1D(hcharname, hchartitle, ECalnBins[amend], ECalnLower[amend], ECalnUpper[amend]);
0216     mehEcaln[amend]->setAxisTitle("Number of RecHits", 1);
0217     mehEcaln[amend]->setAxisTitle("Count", 2);
0218     hcharname = "hEcalRes_" + ECalString[amend];
0219     hchartitle = ECalString[amend] + "  rechit resolution";
0220     mehEcalRes[amend] =
0221         iBooker.book1D(hcharname, hchartitle, ECalResBins[amend], ECalResLower[amend], ECalResUpper[amend]);
0222     mehEcalRes[amend]->setAxisTitle("RecHit E - SimHit E", 1);
0223     mehEcalRes[amend]->setAxisTitle("Count", 2);
0224   }
0225 
0226   //Si Pixels
0227   string SiPixelString[7] = {"BRL1", "BRL2", "BRL3", "FWD1n", "FWD1p", "FWD2n", "FWD2p"};
0228   for (int j = 0; j < 7; ++j) {
0229     mehSiPixeln[j] = nullptr;
0230     mehSiPixelResX[j] = nullptr;
0231     mehSiPixelResY[j] = nullptr;
0232   }
0233 
0234   iBooker.setCurrentFolder("GlobalRecHitsV/SiPixels");
0235   for (int amend = 0; amend < 7; ++amend) {
0236     hcharname = "hSiPixeln_" + SiPixelString[amend];
0237     hchartitle = SiPixelString[amend] + " rechits";
0238     mehSiPixeln[amend] = iBooker.book1D(hcharname, hchartitle, 200, 0., 200.);
0239     mehSiPixeln[amend]->setAxisTitle("Number of hits in " + SiPixelString[amend], 1);
0240     mehSiPixeln[amend]->setAxisTitle("Count", 2);
0241     hcharname = "hSiPixelResX_" + SiPixelString[amend];
0242     hchartitle = SiPixelString[amend] + " rechit x resolution";
0243     mehSiPixelResX[amend] = iBooker.book1D(hcharname, hchartitle, 200, -0.02, .02);
0244     mehSiPixelResX[amend]->setAxisTitle("X-resolution in " + SiPixelString[amend], 1);
0245     mehSiPixelResX[amend]->setAxisTitle("Count", 2);
0246     hcharname = "hSiPixelResY_" + SiPixelString[amend];
0247     hchartitle = SiPixelString[amend] + " rechit y resolution";
0248 
0249     mehSiPixelResY[amend] = iBooker.book1D(hcharname, hchartitle, 200, -0.02, .02);
0250     mehSiPixelResY[amend]->setAxisTitle("Y-resolution in " + SiPixelString[amend], 1);
0251     mehSiPixelResY[amend]->setAxisTitle("Count", 2);
0252   }
0253 
0254   //Muons
0255   iBooker.setCurrentFolder("GlobalRecHitsV/Muons");
0256 
0257   mehDtMuonn = nullptr;
0258   mehCSCn = nullptr;
0259   mehRPCn = nullptr;
0260 
0261   string n_List[3] = {"hDtMuonn", "hCSCn", "hRPCn"};
0262   string hist_string[3] = {"Dt", "CSC", "RPC"};
0263 
0264   for (int amend = 0; amend < 3; ++amend) {
0265     hchartitle = hist_string[amend] + " rechits";
0266     if (amend == 0) {
0267       mehDtMuonn = iBooker.book1D(n_List[amend], hchartitle, 50, 0., 500.);
0268       mehDtMuonn->setAxisTitle("Number of Rechits", 1);
0269       mehDtMuonn->setAxisTitle("Count", 2);
0270     }
0271     if (amend == 1) {
0272       mehCSCn = iBooker.book1D(n_List[amend], hchartitle, 50, 0., 500.);
0273       mehCSCn->setAxisTitle("Number of Rechits", 1);
0274       mehCSCn->setAxisTitle("Count", 2);
0275     }
0276     if (amend == 2) {
0277       mehRPCn = iBooker.book1D(n_List[amend], hchartitle, 50, 0., 500.);
0278       mehRPCn->setAxisTitle("Number of Rechits", 1);
0279       mehRPCn->setAxisTitle("Count", 2);
0280     }
0281   }
0282 
0283   mehDtMuonRes = nullptr;
0284   mehCSCResRDPhi = nullptr;
0285   mehRPCResX = nullptr;
0286 
0287   hcharname = "hDtMuonRes";
0288   hchartitle = "DT wire distance resolution";
0289   mehDtMuonRes = iBooker.book1D(hcharname, hchartitle, 200, -0.2, 0.2);
0290   hcharname = "CSCResRDPhi";
0291   hchartitle = "CSC perp*dphi resolution";
0292   mehCSCResRDPhi = iBooker.book1D(hcharname, hchartitle, 200, -0.2, 0.2);
0293   hcharname = "hRPCResX";
0294   hchartitle = "RPC rechits x resolution";
0295   mehRPCResX = iBooker.book1D(hcharname, hchartitle, 50, -5., 5.);
0296 }
0297 
0298 void GlobalRecHitsAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0299   std::string MsgLoggerCat = "GlobalRecHitsAnalyzer_analyze";
0300 
0301   // keep track of number of events processed
0302   ++count;
0303 
0304   // get event id information
0305   edm::RunNumber_t nrun = iEvent.id().run();
0306   edm::EventNumber_t nevt = iEvent.id().event();
0307 
0308   if (verbosity > 0) {
0309     edm::LogInfo(MsgLoggerCat) << "Processing run " << nrun << ", event " << nevt << " (" << count << " events total)";
0310   } else if (verbosity == 0) {
0311     if (nevt % frequency == 0 || nevt == 1) {
0312       edm::LogInfo(MsgLoggerCat) << "Processing run " << nrun << ", event " << nevt << " (" << count
0313                                  << " events total)";
0314     }
0315   }
0316 
0317   // look at information available in the event
0318   if (getAllProvenances) {
0319     std::vector<const edm::StableProvenance*> AllProv;
0320     iEvent.getAllStableProvenance(AllProv);
0321 
0322     if (verbosity >= 0)
0323       edm::LogInfo(MsgLoggerCat) << "Number of Provenances = " << AllProv.size();
0324 
0325     if (printProvenanceInfo && (verbosity >= 0)) {
0326       TString eventout("\nProvenance info:\n");
0327 
0328       for (unsigned int i = 0; i < AllProv.size(); ++i) {
0329         eventout += "\n       ******************************";
0330         eventout += "\n       Module       : ";
0331         eventout += AllProv[i]->moduleLabel();
0332         eventout += "\n       ProductID    : ";
0333         eventout += AllProv[i]->productID().id();
0334         eventout += "\n       ClassName    : ";
0335         eventout += AllProv[i]->className();
0336         eventout += "\n       InstanceName : ";
0337         eventout += AllProv[i]->productInstanceName();
0338         eventout += "\n       BranchName   : ";
0339         eventout += AllProv[i]->branchName();
0340       }
0341       eventout += "\n       ******************************\n";
0342       edm::LogInfo(MsgLoggerCat) << eventout << "\n";
0343       printProvenanceInfo = false;
0344     }
0345     getAllProvenances = false;
0346   }
0347 
0348   // call fill functions
0349   // gather Ecal information from event
0350   fillECal(iEvent, iSetup);
0351   // gather Hcal information from event
0352   fillHCal(iEvent, iSetup);
0353   // gather Track information from event
0354   fillTrk(iEvent, iSetup);
0355   // gather Muon information from event
0356   fillMuon(iEvent, iSetup);
0357 
0358   if (verbosity > 0)
0359     edm::LogInfo(MsgLoggerCat) << "Done gathering data from event.";
0360 
0361   return;
0362 }
0363 
0364 void GlobalRecHitsAnalyzer::fillECal(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0365   std::string MsgLoggerCat = "GlobalRecHitsAnalyzer_fillECal";
0366 
0367   TString eventout;
0368   if (verbosity > 0)
0369     eventout = "\nGathering info:";
0370 
0371   // extract crossing frame from event
0372   edm::Handle<CrossingFrame<PCaloHit>> crossingFrame;
0373 
0374   ////////////////////////
0375   //extract EB information
0376   ////////////////////////
0377   edm::Handle<EBUncalibratedRecHitCollection> EcalUncalibRecHitEB;
0378   iEvent.getByToken(ECalUncalEBSrc_Token_, EcalUncalibRecHitEB);
0379   bool validUncalibRecHitEB = EcalUncalibRecHitEB.isValid();
0380   if (!validUncalibRecHitEB) {
0381     LogDebug(MsgLoggerCat) << "Unable to find EcalUncalRecHitEB in event!";
0382   }
0383 
0384   edm::Handle<EBRecHitCollection> EcalRecHitEB;
0385   iEvent.getByToken(ECalEBSrc_Token_, EcalRecHitEB);
0386   bool validRecHitEB = EcalRecHitEB.isValid();
0387   if (!validRecHitEB) {
0388     LogDebug(MsgLoggerCat) << "Unable to find EcalRecHitEB in event!";
0389   }
0390 
0391   // loop over simhits
0392   iEvent.getByToken(EBHits_Token_, crossingFrame);
0393   bool validXFrame = crossingFrame.isValid();
0394   if (!validXFrame) {
0395     LogDebug(MsgLoggerCat) << "Unable to find cal barrel crossingFrame in event!";
0396   }
0397 
0398   MapType ebSimMap;
0399   if (validXFrame) {
0400     const MixCollection<PCaloHit> barrelHits(crossingFrame.product());
0401     // keep track of sum of simhit energy in each crystal
0402     for (auto const& iHit : barrelHits) {
0403       EBDetId ebid = EBDetId(iHit.id());
0404 
0405       uint32_t crystid = ebid.rawId();
0406       ebSimMap[crystid] += iHit.energy();
0407     }
0408   }
0409 
0410   int nEBRecHits = 0;
0411   // loop over RecHits
0412   if (validUncalibRecHitEB && validRecHitEB) {
0413     const EBUncalibratedRecHitCollection* EBUncalibRecHit = EcalUncalibRecHitEB.product();
0414     const EBRecHitCollection* EBRecHit = EcalRecHitEB.product();
0415 
0416     for (EcalUncalibratedRecHitCollection::const_iterator uncalibRecHit = EBUncalibRecHit->begin();
0417          uncalibRecHit != EBUncalibRecHit->end();
0418          ++uncalibRecHit) {
0419       EBDetId EBid = EBDetId(uncalibRecHit->id());
0420 
0421       EcalRecHitCollection::const_iterator myRecHit = EBRecHit->find(EBid);
0422 
0423       if (myRecHit != EBRecHit->end()) {
0424         ++nEBRecHits;
0425         mehEcalRes[1]->Fill(myRecHit->energy() - ebSimMap[EBid.rawId()]);
0426       }
0427     }
0428 
0429     if (verbosity > 1) {
0430       eventout += "\n          Number of EBRecHits collected:............ ";
0431       eventout += nEBRecHits;
0432     }
0433     mehEcaln[1]->Fill((float)nEBRecHits);
0434   }
0435 
0436   ////////////////////////
0437   //extract EE information
0438   ////////////////////////
0439   edm::Handle<EEUncalibratedRecHitCollection> EcalUncalibRecHitEE;
0440   iEvent.getByToken(ECalUncalEESrc_Token_, EcalUncalibRecHitEE);
0441   bool validuncalibRecHitEE = EcalUncalibRecHitEE.isValid();
0442   if (!validuncalibRecHitEE) {
0443     LogDebug(MsgLoggerCat) << "Unable to find EcalUncalRecHitEE in event!";
0444   }
0445 
0446   edm::Handle<EERecHitCollection> EcalRecHitEE;
0447   iEvent.getByToken(ECalEESrc_Token_, EcalRecHitEE);
0448   bool validRecHitEE = EcalRecHitEE.isValid();
0449   if (!validRecHitEE) {
0450     LogDebug(MsgLoggerCat) << "Unable to find EcalRecHitEE in event!";
0451   }
0452 
0453   // loop over simhits
0454   iEvent.getByToken(EEHits_Token_, crossingFrame);
0455   validXFrame = crossingFrame.isValid();
0456   if (!validXFrame) {
0457     LogDebug(MsgLoggerCat) << "Unable to find cal endcap crossingFrame in event!";
0458   }
0459 
0460   MapType eeSimMap;
0461   if (validXFrame) {
0462     const MixCollection<PCaloHit> endcapHits(crossingFrame.product());
0463     // keep track of sum of simhit energy in each crystal
0464     for (auto const& iHit : endcapHits) {
0465       EEDetId eeid = EEDetId(iHit.id());
0466 
0467       uint32_t crystid = eeid.rawId();
0468       eeSimMap[crystid] += iHit.energy();
0469     }
0470   }
0471 
0472   int nEERecHits = 0;
0473   if (validuncalibRecHitEE && validRecHitEE) {
0474     // loop over RecHits
0475     const EEUncalibratedRecHitCollection* EEUncalibRecHit = EcalUncalibRecHitEE.product();
0476     const EERecHitCollection* EERecHit = EcalRecHitEE.product();
0477 
0478     for (EcalUncalibratedRecHitCollection::const_iterator uncalibRecHit = EEUncalibRecHit->begin();
0479          uncalibRecHit != EEUncalibRecHit->end();
0480          ++uncalibRecHit) {
0481       EEDetId EEid = EEDetId(uncalibRecHit->id());
0482 
0483       EcalRecHitCollection::const_iterator myRecHit = EERecHit->find(EEid);
0484 
0485       if (myRecHit != EERecHit->end()) {
0486         ++nEERecHits;
0487         mehEcalRes[0]->Fill(myRecHit->energy() - eeSimMap[EEid.rawId()]);
0488       }
0489     }
0490 
0491     if (verbosity > 1) {
0492       eventout += "\n          Number of EERecHits collected:............ ";
0493       eventout += nEERecHits;
0494     }
0495     mehEcaln[0]->Fill((float)nEERecHits);
0496   }
0497 
0498   ////////////////////////
0499   //extract ES information
0500   ////////////////////////
0501   edm::Handle<ESRecHitCollection> EcalRecHitES;
0502   iEvent.getByToken(ECalESSrc_Token_, EcalRecHitES);
0503   bool validRecHitES = EcalRecHitES.isValid();
0504   if (!validRecHitES) {
0505     LogDebug(MsgLoggerCat) << "Unable to find EcalRecHitES in event!";
0506   }
0507 
0508   // loop over simhits
0509   iEvent.getByToken(ESHits_Token_, crossingFrame);
0510   validXFrame = crossingFrame.isValid();
0511   if (!validXFrame) {
0512     LogDebug(MsgLoggerCat) << "Unable to find cal preshower crossingFrame in event!";
0513   }
0514 
0515   MapType esSimMap;
0516   if (validXFrame) {
0517     const MixCollection<PCaloHit> preshowerHits(crossingFrame.product());
0518     // keep track of sum of simhit energy in each crystal
0519     for (auto const& iHit : preshowerHits) {
0520       ESDetId esid = ESDetId(iHit.id());
0521 
0522       uint32_t crystid = esid.rawId();
0523       esSimMap[crystid] += iHit.energy();
0524     }
0525   }
0526 
0527   int nESRecHits = 0;
0528   if (validRecHitES) {
0529     // loop over RecHits
0530     const ESRecHitCollection* ESRecHit = EcalRecHitES.product();
0531     for (EcalRecHitCollection::const_iterator recHit = ESRecHit->begin(); recHit != ESRecHit->end(); ++recHit) {
0532       ESDetId ESid = ESDetId(recHit->id());
0533 
0534       ++nESRecHits;
0535       mehEcalRes[2]->Fill(recHit->energy() - esSimMap[ESid.rawId()]);
0536     }
0537 
0538     if (verbosity > 1) {
0539       eventout += "\n          Number of ESRecHits collected:............ ";
0540       eventout += nESRecHits;
0541     }
0542     mehEcaln[2]->Fill(float(nESRecHits));
0543   }
0544 
0545   if (verbosity > 0)
0546     edm::LogInfo(MsgLoggerCat) << eventout << "\n";
0547 
0548   return;
0549 }
0550 
0551 void GlobalRecHitsAnalyzer::fillHCal(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0552   std::string MsgLoggerCat = "GlobalRecHitsAnalyzer_fillHCal";
0553 
0554   TString eventout;
0555   if (verbosity > 0)
0556     eventout = "\nGathering info:";
0557 
0558   // get geometry
0559   const auto& geometry = iSetup.getHandle(caloGeomToken_);
0560   if (!geometry.isValid()) {
0561     edm::LogWarning(MsgLoggerCat) << "Unable to find CaloGeometry in event!";
0562     return;
0563   }
0564 
0565   ///////////////////////
0566   // extract simhit info
0567   //////////////////////
0568   edm::Handle<edm::PCaloHitContainer> hcalHits;
0569   iEvent.getByToken(HCalSrc_Token_, hcalHits);
0570   bool validhcalHits = hcalHits.isValid();
0571   if (!validhcalHits) {
0572     LogDebug(MsgLoggerCat) << "Unable to find hcalHits in event!";
0573   }
0574 
0575   std::map<HcalDetId, float> fHBEnergySimHits;
0576   std::map<HcalDetId, float> fHEEnergySimHits;
0577   std::map<HcalDetId, float> fHOEnergySimHits;
0578   std::map<HcalDetId, float> fHFEnergySimHits;
0579   if (validhcalHits) {
0580     const edm::PCaloHitContainer* simhitResult = hcalHits.product();
0581 
0582     for (std::vector<PCaloHit>::const_iterator simhits = simhitResult->begin(); simhits != simhitResult->end();
0583          ++simhits) {
0584       HcalDetId detId(simhits->id());
0585 
0586       if (detId.subdet() == sdHcalBrl) {
0587         fHBEnergySimHits[detId] += simhits->energy();
0588       }
0589       if (detId.subdet() == sdHcalEC) {
0590         fHEEnergySimHits[detId] += simhits->energy();
0591       }
0592       if (detId.subdet() == sdHcalOut) {
0593         fHOEnergySimHits[detId] += simhits->energy();
0594       }
0595       if (detId.subdet() == sdHcalFwd) {
0596         fHFEnergySimHits[detId] += simhits->energy();
0597       }
0598     }
0599   }
0600 
0601   ////////////////////////
0602   // get HBHE information
0603   ///////////////////////
0604   std::vector<edm::Handle<HBHERecHitCollection>> hbhe;
0605   iEvent.getManyByType(hbhe);
0606 
0607   bool validHBHE = hbhe[0].isValid();
0608 
0609   if (!validHBHE) {
0610     LogDebug(MsgLoggerCat) << "Unable to find any HBHERecHitCollections in event!";
0611   }
0612 
0613   else {
0614     std::vector<edm::Handle<HBHERecHitCollection>>::iterator ihbhe;
0615     int iHB = 0;
0616     int iHE = 0;
0617     for (ihbhe = hbhe.begin(); ihbhe != hbhe.end(); ++ihbhe) {
0618       for (HBHERecHitCollection::const_iterator jhbhe = (*ihbhe)->begin(); jhbhe != (*ihbhe)->end(); ++jhbhe) {
0619         HcalDetId cell(jhbhe->id());
0620 
0621         if (cell.subdet() == sdHcalBrl) {
0622           ++iHB;
0623           mehHcalRes[0]->Fill(jhbhe->energy() - fHBEnergySimHits[cell]);
0624         }
0625 
0626         else if (cell.subdet() == sdHcalEC) {
0627           ++iHE;
0628           mehHcalRes[1]->Fill(jhbhe->energy() - fHEEnergySimHits[cell]);
0629         }
0630       }
0631     }  // end loop through collection
0632 
0633     if (verbosity > 1) {
0634       eventout += "\n          Number of HBRecHits collected:............ ";
0635       eventout += iHB;
0636     }
0637 
0638     if (verbosity > 1) {
0639       eventout += "\n          Number of HERecHits collected:............ ";
0640       eventout += iHE;
0641     }
0642     mehHcaln[0]->Fill((float)iHB);
0643     mehHcaln[1]->Fill((float)iHE);
0644   }
0645 
0646   ////////////////////////
0647   // get HF information
0648   ///////////////////////
0649   std::vector<edm::Handle<HFRecHitCollection>> hf;
0650   iEvent.getManyByType(hf);
0651   bool validHF = hf[0].isValid();
0652   if (!validHF) {
0653     LogDebug(MsgLoggerCat) << "Unable to find any HFRecHitCollections in event!";
0654   } else {
0655     std::vector<edm::Handle<HFRecHitCollection>>::iterator ihf;
0656 
0657     int iHF = 0;
0658     for (ihf = hf.begin(); ihf != hf.end(); ++ihf) {
0659       for (HFRecHitCollection::const_iterator jhf = (*ihf)->begin(); jhf != (*ihf)->end(); ++jhf) {
0660         HcalDetId cell(jhf->id());
0661 
0662         if (cell.subdet() == sdHcalFwd) {
0663           ++iHF;
0664           mehHcalRes[2]->Fill(jhf->energy() - fHFEnergySimHits[cell]);
0665         }
0666       }
0667     }  // end loop through collection
0668 
0669     if (verbosity > 1) {
0670       eventout += "\n          Number of HFDigis collected:.............. ";
0671       eventout += iHF;
0672     }
0673     mehHcaln[2]->Fill((float)iHF);
0674   }
0675 
0676   ////////////////////////
0677   // get HO information
0678   ///////////////////////
0679   std::vector<edm::Handle<HORecHitCollection>> ho;
0680   iEvent.getManyByType(ho);
0681   bool validHO = ho[0].isValid();
0682   if (!validHO) {
0683     LogDebug(MsgLoggerCat) << "Unable to find any HORecHitCollections in event!";
0684   } else {
0685     std::vector<edm::Handle<HORecHitCollection>>::iterator iho;
0686 
0687     int iHO = 0;
0688     for (iho = ho.begin(); iho != ho.end(); ++iho) {
0689       for (HORecHitCollection::const_iterator jho = (*iho)->begin(); jho != (*iho)->end(); ++jho) {
0690         HcalDetId cell(jho->id());
0691 
0692         if (cell.subdet() == sdHcalOut) {
0693           ++iHO;
0694           mehHcalRes[3]->Fill(jho->energy() - fHOEnergySimHits[cell]);
0695         }
0696       }
0697     }  // end loop through collection
0698 
0699     if (verbosity > 1) {
0700       eventout += "\n          Number of HODigis collected:.............. ";
0701       eventout += iHO;
0702     }
0703     mehHcaln[3]->Fill((float)iHO);
0704   }
0705 
0706   if (verbosity > 0)
0707     edm::LogInfo(MsgLoggerCat) << eventout << "\n";
0708 
0709   return;
0710 }
0711 
0712 void GlobalRecHitsAnalyzer::fillTrk(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0713   //Retrieve tracker topology from geometry
0714   const TrackerTopology* const tTopo = &iSetup.getData(tTopoToken_);
0715   std::string MsgLoggerCat = "GlobalRecHitsAnalyzer_fillTrk";
0716   TString eventout;
0717   if (verbosity > 0)
0718     eventout = "\nGathering info:";
0719 
0720   // get strip information
0721   edm::Handle<SiStripMatchedRecHit2DCollection> rechitsmatched;
0722   iEvent.getByToken(SiStripSrc_Token_, rechitsmatched);
0723   bool validstrip = rechitsmatched.isValid();
0724   if (!validstrip) {
0725     LogDebug(MsgLoggerCat) << "Unable to find stripmatchedrechits in event!";
0726   }
0727 
0728   TrackerHitAssociator associate(iEvent, trackerHitAssociatorConfig_);
0729 
0730   const auto& tGeomHandle = iSetup.getHandle(tGeomToken_);
0731   if (!tGeomHandle.isValid()) {
0732     edm::LogWarning(MsgLoggerCat) << "Unable to find TrackerDigiGeometry in event!";
0733     return;
0734   }
0735   const TrackerGeometry& tracker(*tGeomHandle);
0736 
0737   if (validstrip) {
0738     int nStripBrl = 0, nStripFwd = 0;
0739 
0740     // loop over det units
0741     for (TrackerGeometry::DetContainer::const_iterator it = tGeomHandle->dets().begin();
0742          it != tGeomHandle->dets().end();
0743          ++it) {
0744       uint32_t myid = ((*it)->geographicalId()).rawId();
0745       DetId detid = ((*it)->geographicalId());
0746 
0747       //loop over rechits-matched in the same subdetector
0748       SiStripMatchedRecHit2DCollection::const_iterator rechitmatchedMatch = rechitsmatched->find(detid);
0749 
0750       if (rechitmatchedMatch != rechitsmatched->end()) {
0751         SiStripMatchedRecHit2DCollection::DetSet rechitmatchedRange = *rechitmatchedMatch;
0752         SiStripMatchedRecHit2DCollection::DetSet::const_iterator rechitmatchedRangeIteratorBegin =
0753             rechitmatchedRange.begin();
0754         SiStripMatchedRecHit2DCollection::DetSet::const_iterator rechitmatchedRangeIteratorEnd =
0755             rechitmatchedRange.end();
0756         SiStripMatchedRecHit2DCollection::DetSet::const_iterator itermatched;
0757 
0758         for (itermatched = rechitmatchedRangeIteratorBegin; itermatched != rechitmatchedRangeIteratorEnd;
0759              ++itermatched) {
0760           SiStripMatchedRecHit2D const rechit = *itermatched;
0761           LocalPoint position = rechit.localPosition();
0762 
0763           float mindist = 999999.;
0764           float distx = 999999.;
0765           float disty = 999999.;
0766           float dist = 999999.;
0767           std::pair<LocalPoint, LocalVector> closestPair;
0768           matched.clear();
0769 
0770           float rechitmatchedx = position.x();
0771           float rechitmatchedy = position.y();
0772 
0773           matched = associate.associateHit(rechit);
0774 
0775           if (!matched.empty()) {
0776             //project simhit;
0777             const GluedGeomDet* gluedDet = (const GluedGeomDet*)tracker.idToDet(rechit.geographicalId());
0778             const StripGeomDetUnit* partnerstripdet = (StripGeomDetUnit*)gluedDet->stereoDet();
0779             std::pair<LocalPoint, LocalVector> hitPair;
0780 
0781             for (std::vector<PSimHit>::const_iterator m = matched.begin(); m != matched.end(); m++) {
0782               //project simhit;
0783               hitPair = projectHit((*m), partnerstripdet, gluedDet->surface());
0784               distx = fabs(rechitmatchedx - hitPair.first.x());
0785               disty = fabs(rechitmatchedy - hitPair.first.y());
0786               dist = sqrt(distx * distx + disty * disty);
0787 
0788               if (dist < mindist) {
0789                 mindist = dist;
0790                 closestPair = hitPair;
0791               }
0792             }
0793 
0794             // get TIB
0795             if (detid.subdetId() == sdSiTIB) {
0796               ++nStripBrl;
0797 
0798               if (tTopo->tibLayer(myid) == 1) {
0799                 mehSiStripResX[8]->Fill(rechitmatchedx - closestPair.first.x());
0800                 mehSiStripResY[8]->Fill(rechitmatchedy - closestPair.first.y());
0801               }
0802               if (tTopo->tibLayer(myid) == 2) {
0803                 mehSiStripResX[9]->Fill(rechitmatchedx - closestPair.first.x());
0804                 mehSiStripResY[9]->Fill(rechitmatchedy - closestPair.first.y());
0805               }
0806               if (tTopo->tibLayer(myid) == 3) {
0807                 mehSiStripResX[10]->Fill(rechitmatchedx - closestPair.first.x());
0808                 mehSiStripResY[10]->Fill(rechitmatchedy - closestPair.first.y());
0809               }
0810               if (tTopo->tibLayer(myid) == 4) {
0811                 mehSiStripResX[11]->Fill(rechitmatchedx - closestPair.first.x());
0812                 mehSiStripResY[11]->Fill(rechitmatchedy - closestPair.first.y());
0813               }
0814             }
0815 
0816             // get TOB
0817             if (detid.subdetId() == sdSiTOB) {
0818               ++nStripBrl;
0819 
0820               if (tTopo->tobLayer(myid) == 1) {
0821                 mehSiStripResX[15]->Fill(rechitmatchedx - closestPair.first.x());
0822                 mehSiStripResY[15]->Fill(rechitmatchedy - closestPair.first.y());
0823               }
0824               if (tTopo->tobLayer(myid) == 2) {
0825                 mehSiStripResX[16]->Fill(rechitmatchedx - closestPair.first.x());
0826                 mehSiStripResY[16]->Fill(rechitmatchedy - closestPair.first.y());
0827               }
0828               if (tTopo->tobLayer(myid) == 3) {
0829                 mehSiStripResX[17]->Fill(rechitmatchedx - closestPair.first.x());
0830                 mehSiStripResY[17]->Fill(rechitmatchedy - closestPair.first.y());
0831               }
0832               if (tTopo->tobLayer(myid) == 4) {
0833                 mehSiStripResX[18]->Fill(rechitmatchedx - closestPair.first.x());
0834                 mehSiStripResY[18]->Fill(rechitmatchedy - closestPair.first.y());
0835               }
0836             }
0837 
0838             // get TID
0839             if (detid.subdetId() == sdSiTID) {
0840               ++nStripFwd;
0841 
0842               if (tTopo->tidWheel(myid) == 1) {
0843                 mehSiStripResX[12]->Fill(rechitmatchedx - closestPair.first.x());
0844                 mehSiStripResY[12]->Fill(rechitmatchedy - closestPair.first.y());
0845               }
0846               if (tTopo->tidWheel(myid) == 2) {
0847                 mehSiStripResX[13]->Fill(rechitmatchedx - closestPair.first.x());
0848                 mehSiStripResY[13]->Fill(rechitmatchedy - closestPair.first.y());
0849               }
0850               if (tTopo->tidWheel(myid) == 3) {
0851                 mehSiStripResX[14]->Fill(rechitmatchedx - closestPair.first.x());
0852                 mehSiStripResY[14]->Fill(rechitmatchedy - closestPair.first.y());
0853               }
0854             }
0855 
0856             // get TEC
0857             if (detid.subdetId() == sdSiTEC) {
0858               ++nStripFwd;
0859 
0860               if (tTopo->tecWheel(myid) == 1) {
0861                 mehSiStripResX[0]->Fill(rechitmatchedx - closestPair.first.x());
0862                 mehSiStripResY[0]->Fill(rechitmatchedy - closestPair.first.y());
0863               }
0864               if (tTopo->tecWheel(myid) == 2) {
0865                 mehSiStripResX[1]->Fill(rechitmatchedx - closestPair.first.x());
0866                 mehSiStripResY[1]->Fill(rechitmatchedy - closestPair.first.y());
0867               }
0868               if (tTopo->tecWheel(myid) == 3) {
0869                 mehSiStripResX[2]->Fill(rechitmatchedx - closestPair.first.x());
0870                 mehSiStripResY[2]->Fill(rechitmatchedy - closestPair.first.y());
0871               }
0872               if (tTopo->tecWheel(myid) == 4) {
0873                 mehSiStripResX[3]->Fill(rechitmatchedx - closestPair.first.x());
0874                 mehSiStripResY[3]->Fill(rechitmatchedy - closestPair.first.y());
0875               }
0876               if (tTopo->tecWheel(myid) == 5) {
0877                 mehSiStripResX[4]->Fill(rechitmatchedx - closestPair.first.x());
0878                 mehSiStripResY[4]->Fill(rechitmatchedy - closestPair.first.y());
0879               }
0880               if (tTopo->tecWheel(myid) == 6) {
0881                 mehSiStripResX[5]->Fill(rechitmatchedx - closestPair.first.x());
0882                 mehSiStripResY[5]->Fill(rechitmatchedy - closestPair.first.y());
0883               }
0884               if (tTopo->tecWheel(myid) == 7) {
0885                 mehSiStripResX[6]->Fill(rechitmatchedx - closestPair.first.x());
0886                 mehSiStripResY[6]->Fill(rechitmatchedy - closestPair.first.y());
0887               }
0888               if (tTopo->tecWheel(myid) == 8) {
0889                 mehSiStripResX[7]->Fill(rechitmatchedx - closestPair.first.x());
0890                 mehSiStripResY[7]->Fill(rechitmatchedy - closestPair.first.y());
0891               }
0892             }
0893 
0894           }  // end if matched empty
0895         }
0896       }
0897     }  // end loop over det units
0898 
0899     if (verbosity > 1) {
0900       eventout += "\n          Number of BrlStripRecHits collected:...... ";
0901       eventout += nStripBrl;
0902     }
0903 
0904     for (int i = 8; i < 12; ++i) {
0905       mehSiStripn[i]->Fill((float)nStripBrl);
0906     }
0907     for (int i = 16; i < 19; ++i) {
0908       mehSiStripn[i]->Fill((float)nStripBrl);
0909     }
0910 
0911     if (verbosity > 1) {
0912       eventout += "\n          Number of FrwdStripRecHits collected:..... ";
0913       eventout += nStripFwd;
0914     }
0915     for (int i = 0; i < 8; ++i) {
0916       mehSiStripn[i]->Fill((float)nStripFwd);
0917     }
0918     for (int i = 12; i < 16; ++i) {
0919       mehSiStripn[i]->Fill((float)nStripFwd);
0920     }
0921   }
0922 
0923   // get pixel information
0924   //Get RecHits
0925   edm::Handle<SiPixelRecHitCollection> recHitColl;
0926   iEvent.getByToken(SiPxlSrc_Token_, recHitColl);
0927   bool validpixel = recHitColl.isValid();
0928   if (!validpixel) {
0929     LogDebug(MsgLoggerCat) << "Unable to find SiPixelRecHitCollection in event!";
0930   } else {
0931     int nPxlBrl = 0, nPxlFwd = 0;
0932     //iterate over detunits
0933     for (TrackerGeometry::DetContainer::const_iterator it = tGeomHandle->dets().begin();
0934          it != tGeomHandle->dets().end();
0935          ++it) {
0936       uint32_t myid = ((*it)->geographicalId()).rawId();
0937       DetId detId = ((*it)->geographicalId());
0938       int subid = detId.subdetId();
0939 
0940       if (!((subid == sdPxlBrl) || (subid == sdPxlFwd)))
0941         continue;
0942 
0943       SiPixelRecHitCollection::const_iterator pixeldet = recHitColl->find(detId);
0944       if (pixeldet == recHitColl->end())
0945         continue;
0946       SiPixelRecHitCollection::DetSet pixelrechitRange = *pixeldet;
0947       SiPixelRecHitCollection::DetSet::const_iterator pixelrechitRangeIteratorBegin = pixelrechitRange.begin();
0948       SiPixelRecHitCollection::DetSet::const_iterator pixelrechitRangeIteratorEnd = pixelrechitRange.end();
0949       SiPixelRecHitCollection::DetSet::const_iterator pixeliter = pixelrechitRangeIteratorBegin;
0950 
0951       std::vector<PSimHit> matched;
0952 
0953       //----Loop over rechits for this detId
0954       for (; pixeliter != pixelrechitRangeIteratorEnd; ++pixeliter) {
0955         matched.clear();
0956         matched = associate.associateHit(*pixeliter);
0957 
0958         if (!matched.empty()) {
0959           float closest = 9999.9;
0960           LocalPoint lp = pixeliter->localPosition();
0961           float rechit_x = lp.x();
0962           float rechit_y = lp.y();
0963 
0964           float sim_x = 0.;
0965           float sim_y = 0.;
0966 
0967           //loop over sim hits and fill closet
0968           for (std::vector<PSimHit>::const_iterator m = matched.begin(); m != matched.end(); ++m) {
0969             float sim_x1 = (*m).entryPoint().x();
0970             float sim_x2 = (*m).exitPoint().x();
0971             float sim_xpos = 0.5 * (sim_x1 + sim_x2);
0972 
0973             float sim_y1 = (*m).entryPoint().y();
0974             float sim_y2 = (*m).exitPoint().y();
0975             float sim_ypos = 0.5 * (sim_y1 + sim_y2);
0976 
0977             float x_res = fabs(sim_xpos - rechit_x);
0978             float y_res = fabs(sim_ypos - rechit_y);
0979 
0980             float dist = sqrt(x_res * x_res + y_res * y_res);
0981 
0982             if (dist < closest) {
0983               closest = dist;
0984               sim_x = sim_xpos;
0985               sim_y = sim_ypos;
0986             }
0987           }  // end sim hit loop
0988 
0989           // get Barrel pixels ***************Pixel STuff******************
0990           if (subid == sdPxlBrl) {
0991             ++nPxlBrl;
0992 
0993             if (tTopo->pxbLayer(myid) == 1) {
0994               mehSiPixelResX[0]->Fill(rechit_x - sim_x);
0995               mehSiPixelResY[0]->Fill(rechit_y - sim_y);
0996             }
0997             if (tTopo->pxbLayer(myid) == 2) {
0998               mehSiPixelResX[1]->Fill(rechit_x - sim_x);
0999               mehSiPixelResY[1]->Fill(rechit_y - sim_y);
1000             }
1001             if (tTopo->pxbLayer(myid) == 3) {
1002               mehSiPixelResX[2]->Fill(rechit_x - sim_x);
1003               mehSiPixelResY[2]->Fill(rechit_y - sim_y);
1004             }
1005           }
1006 
1007           // get Forward pixels
1008           if (subid == sdPxlFwd) {
1009             ++nPxlFwd;
1010 
1011             if (tTopo->pxfDisk(myid) == 1) {
1012               if (tTopo->pxfSide(myid) == 1) {
1013                 mehSiPixelResX[3]->Fill(rechit_x - sim_x);
1014                 mehSiPixelResY[3]->Fill(rechit_y - sim_y);
1015               }
1016               if (tTopo->pxfSide(myid) == 2) {
1017                 mehSiPixelResX[4]->Fill(rechit_x - sim_x);
1018                 mehSiPixelResY[4]->Fill(rechit_y - sim_y);
1019               }
1020             }
1021             if (tTopo->pxfDisk(myid) == 2) {
1022               if (tTopo->pxfSide(myid) == 1) {
1023                 mehSiPixelResX[5]->Fill(rechit_x - sim_x);
1024                 mehSiPixelResY[5]->Fill(rechit_y - sim_y);
1025               }
1026               if (tTopo->pxfSide(myid) == 2) {
1027                 mehSiPixelResX[6]->Fill(rechit_x - sim_x);
1028                 mehSiPixelResY[6]->Fill(rechit_y - sim_y);
1029               }
1030             }
1031           }
1032         }  // end matched emtpy
1033       }    // <-----end rechit loop
1034     }      // <------ end detunit loop
1035 
1036     if (verbosity > 1) {
1037       eventout += "\n          Number of BrlPixelRecHits collected:...... ";
1038       eventout += nPxlBrl;
1039     }
1040     for (int i = 0; i < 3; ++i) {
1041       mehSiPixeln[i]->Fill((float)nPxlBrl);
1042     }
1043 
1044     if (verbosity > 1) {
1045       eventout += "\n          Number of FrwdPixelRecHits collected:..... ";
1046       eventout += nPxlFwd;
1047     }
1048 
1049     for (int i = 3; i < 7; ++i) {
1050       mehSiPixeln[i]->Fill((float)nPxlFwd);
1051     }
1052   }
1053 
1054   if (verbosity > 0)
1055     edm::LogInfo(MsgLoggerCat) << eventout << "\n";
1056 
1057   return;
1058 }
1059 
1060 void GlobalRecHitsAnalyzer::fillMuon(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
1061   std::string MsgLoggerCat = "GlobalRecHitsAnalyzer_fillMuon";
1062 
1063   TString eventout;
1064   if (verbosity > 0)
1065     eventout = "\nGathering info:";
1066 
1067   // get DT information
1068   const auto& dtGeom = iSetup.getHandle(dtGeomToken_);
1069   if (!dtGeom.isValid()) {
1070     edm::LogWarning(MsgLoggerCat) << "Unable to find DTMuonGeometryRecord in event!";
1071     return;
1072   }
1073 
1074   edm::Handle<edm::PSimHitContainer> dtsimHits;
1075   iEvent.getByToken(MuDTSimSrc_Token_, dtsimHits);
1076   bool validdtsim = dtsimHits.isValid();
1077   if (!validdtsim) {
1078     LogDebug(MsgLoggerCat) << "Unable to find dtsimHits in event!";
1079   }
1080 
1081   edm::Handle<DTRecHitCollection> dtRecHits;
1082   iEvent.getByToken(MuDTSrc_Token_, dtRecHits);
1083   bool validdtrec = dtRecHits.isValid();
1084   if (!validdtrec) {
1085     LogDebug(MsgLoggerCat) << "Unable to find dtRecHits in event!";
1086   }
1087 
1088   if (validdtsim && validdtrec) {
1089     std::map<DTWireId, edm::PSimHitContainer> simHitsPerWire =
1090         DTHitQualityUtils::mapSimHitsPerWire(*(dtsimHits.product()));
1091 
1092     std::map<DTWireId, std::vector<DTRecHit1DPair>> recHitsPerWire = map1DRecHitsPerWire(dtRecHits.product());
1093 
1094     int nDt = compute(dtGeom.product(), simHitsPerWire, recHitsPerWire, 1);
1095 
1096     if (verbosity > 1) {
1097       eventout += "\n          Number of DtMuonRecHits collected:........ ";
1098       eventout += nDt;
1099     }
1100     mehDtMuonn->Fill(float(nDt));
1101   }
1102 
1103   // get CSC Strip information
1104   // get map of sim hits
1105   theMap.clear();
1106   edm::Handle<CrossingFrame<PSimHit>> cf;
1107 
1108   iEvent.getByToken(MuCSCHits_Token_, cf);
1109   bool validXFrame = cf.isValid();
1110   if (!validXFrame) {
1111     LogDebug(MsgLoggerCat) << "Unable to find muo CSC crossingFrame in event!";
1112   } else {
1113     const MixCollection<PSimHit> simHits(cf.product());
1114 
1115     // arrange the hits by detUnit
1116     for (auto const& iHit : simHits) {
1117       theMap[iHit.detUnitId()].push_back(iHit);
1118     }
1119   }
1120 
1121   // get geometry
1122   const auto& hGeom = iSetup.getHandle(cscGeomToken_);
1123   if (!hGeom.isValid()) {
1124     edm::LogWarning(MsgLoggerCat) << "Unable to find CSCMuonGeometryRecord in event!";
1125     return;
1126   }
1127   const CSCGeometry* theCSCGeometry = &*hGeom;
1128 
1129   // get rechits
1130   edm::Handle<CSCRecHit2DCollection> hRecHits;
1131   iEvent.getByToken(MuCSCSrc_Token_, hRecHits);
1132   bool validCSC = hRecHits.isValid();
1133   if (!validCSC) {
1134     LogDebug(MsgLoggerCat) << "Unable to find CSC RecHits in event!";
1135   } else {
1136     const CSCRecHit2DCollection* cscRecHits = hRecHits.product();
1137 
1138     int nCSC = 0;
1139     for (CSCRecHit2DCollection::const_iterator recHitItr = cscRecHits->begin(); recHitItr != cscRecHits->end();
1140          ++recHitItr) {
1141       int detId = (*recHitItr).cscDetId().rawId();
1142 
1143       edm::PSimHitContainer simHits;
1144       std::map<int, edm::PSimHitContainer>::const_iterator mapItr = theMap.find(detId);
1145       if (mapItr != theMap.end()) {
1146         simHits = mapItr->second;
1147       }
1148 
1149       if (simHits.size() == 1) {
1150         ++nCSC;
1151 
1152         const GeomDetUnit* detUnit = theCSCGeometry->idToDetUnit(CSCDetId(detId));
1153         const CSCLayer* layer = dynamic_cast<const CSCLayer*>(detUnit);
1154 
1155         int chamberType = layer->chamber()->specs()->chamberType();
1156         plotResolution(simHits[0], *recHitItr, layer, chamberType);
1157       }
1158     }
1159 
1160     if (verbosity > 1) {
1161       eventout += "\n          Number of CSCRecHits collected:........... ";
1162       eventout += nCSC;
1163     }
1164     mehCSCn->Fill((float)nCSC);
1165   }
1166 
1167   // get RPC information
1168   std::map<double, int> mapsim, maprec;
1169   std::map<int, double> nmapsim, nmaprec;
1170   const auto& rpcGeom = iSetup.getHandle(rpcGeomToken_);
1171   if (!rpcGeom.isValid()) {
1172     edm::LogWarning(MsgLoggerCat) << "Unable to find RPCMuonGeometryRecord in event!";
1173     return;
1174   }
1175 
1176   edm::Handle<edm::PSimHitContainer> simHit;
1177   iEvent.getByToken(MuRPCSimSrc_Token_, simHit);
1178   bool validrpcsim = simHit.isValid();
1179   if (!validrpcsim) {
1180     LogDebug(MsgLoggerCat) << "Unable to find RPCSimHit in event!";
1181   }
1182 
1183   edm::Handle<RPCRecHitCollection> recHit;
1184   iEvent.getByToken(MuRPCSrc_Token_, recHit);
1185   bool validrpc = recHit.isValid();
1186   if (!validrpc) {
1187     LogDebug(MsgLoggerCat) << "Unable to find RPCRecHit in event!";
1188   }
1189 
1190   if (validrpc) {
1191     int nRPC = 0;
1192     RPCRecHitCollection::const_iterator recIt;
1193     int nrec = 0;
1194     for (recIt = recHit->begin(); recIt != recHit->end(); ++recIt) {
1195       RPCDetId Rid = (RPCDetId)(*recIt).rpcId();
1196       const RPCRoll* roll = dynamic_cast<const RPCRoll*>(rpcGeom->roll(Rid));
1197       if (roll->isForward()) {
1198         if (verbosity > 1) {
1199           eventout += "\n          Number of RPCRecHits collected:........... ";
1200           eventout += nRPC;
1201         }
1202 
1203         if (verbosity > 0)
1204           edm::LogInfo(MsgLoggerCat) << eventout << "\n";
1205         return;
1206       }
1207       nrec = nrec + 1;
1208       LocalPoint rhitlocal = (*recIt).localPosition();
1209       double rhitlocalx = rhitlocal.x();
1210       maprec[rhitlocalx] = nrec;
1211     }
1212 
1213     int i = 0;
1214     for (std::map<double, int>::iterator iter = maprec.begin(); iter != maprec.end(); ++iter) {
1215       i = i + 1;
1216       nmaprec[i] = (*iter).first;
1217     }
1218 
1219     int nsim = 0;
1220     if (validrpcsim) {
1221       edm::PSimHitContainer::const_iterator simIt;
1222       for (simIt = simHit->begin(); simIt != simHit->end(); simIt++) {
1223         int ptype = (*simIt).particleType();
1224         if (ptype == 13 || ptype == -13) {
1225           nsim = nsim + 1;
1226           LocalPoint shitlocal = (*simIt).localPosition();
1227           double shitlocalx = shitlocal.x();
1228           mapsim[shitlocalx] = nsim;
1229         }
1230       }
1231 
1232       i = 0;
1233       for (std::map<double, int>::iterator iter = mapsim.begin(); iter != mapsim.end(); ++iter) {
1234         i = i + 1;
1235         nmapsim[i] = (*iter).first;
1236       }
1237     }
1238 
1239     if (nsim == nrec) {
1240       for (int r = 0; r < nsim; r++) {
1241         ++nRPC;
1242         mehRPCResX->Fill(nmaprec[r + 1] - nmapsim[r + 1]);
1243       }
1244     }
1245 
1246     if (verbosity > 1) {
1247       eventout += "\n          Number of RPCRecHits collected:........... ";
1248       eventout += nRPC;
1249     }
1250     mehRPCn->Fill((float)nRPC);
1251   }
1252 
1253   if (verbosity > 0)
1254     edm::LogInfo(MsgLoggerCat) << eventout << "\n";
1255 
1256   return;
1257 }
1258 
1259 //needed by to do the residual for matched hits in SiStrip
1260 std::pair<LocalPoint, LocalVector> GlobalRecHitsAnalyzer::projectHit(const PSimHit& hit,
1261                                                                      const StripGeomDetUnit* stripDet,
1262                                                                      const BoundPlane& plane) {
1263   const StripTopology& topol = stripDet->specificTopology();
1264   GlobalPoint globalpos = stripDet->surface().toGlobal(hit.localPosition());
1265   LocalPoint localHit = plane.toLocal(globalpos);
1266   //track direction
1267   LocalVector locdir = hit.localDirection();
1268   //rotate track in new frame
1269 
1270   GlobalVector globaldir = stripDet->surface().toGlobal(locdir);
1271   LocalVector dir = plane.toLocal(globaldir);
1272   float scale = -localHit.z() / dir.z();
1273 
1274   LocalPoint projectedPos = localHit + scale * dir;
1275 
1276   float selfAngle = topol.stripAngle(topol.strip(hit.localPosition()));
1277 
1278   // vector along strip in hit frame
1279   LocalVector stripDir(sin(selfAngle), cos(selfAngle), 0);
1280 
1281   LocalVector localStripDir(plane.toLocal(stripDet->surface().toGlobal(stripDir)));
1282 
1283   return std::pair<LocalPoint, LocalVector>(projectedPos, localStripDir);
1284 }
1285 
1286 // Return a map between DTRecHit1DPair and wireId
1287 std::map<DTWireId, std::vector<DTRecHit1DPair>> GlobalRecHitsAnalyzer::map1DRecHitsPerWire(
1288     const DTRecHitCollection* dt1DRecHitPairs) {
1289   std::map<DTWireId, std::vector<DTRecHit1DPair>> ret;
1290 
1291   for (DTRecHitCollection::const_iterator rechit = dt1DRecHitPairs->begin(); rechit != dt1DRecHitPairs->end();
1292        rechit++) {
1293     ret[(*rechit).wireId()].push_back(*rechit);
1294   }
1295 
1296   return ret;
1297 }
1298 
1299 // Compute SimHit distance from wire (cm)
1300 float GlobalRecHitsAnalyzer::simHitDistFromWire(const DTLayer* layer, DTWireId wireId, const PSimHit& hit) {
1301   float xwire = layer->specificTopology().wirePosition(wireId.wire());
1302   LocalPoint entryP = hit.entryPoint();
1303   LocalPoint exitP = hit.exitPoint();
1304   float xEntry = entryP.x() - xwire;
1305   float xExit = exitP.x() - xwire;
1306 
1307   //FIXME: check...
1308   return fabs(xEntry - (entryP.z() * (xExit - xEntry)) / (exitP.z() - entryP.z()));
1309 }
1310 
1311 // Find the RecHit closest to the muon SimHit
1312 template <typename type>
1313 const type* GlobalRecHitsAnalyzer::findBestRecHit(const DTLayer* layer,
1314                                                   DTWireId wireId,
1315                                                   const std::vector<type>& recHits,
1316                                                   const float simHitDist) {
1317   float res = 99999;
1318   const type* theBestRecHit = nullptr;
1319   // Loop over RecHits within the cell
1320   for (typename std::vector<type>::const_iterator recHit = recHits.begin(); recHit != recHits.end(); recHit++) {
1321     float distTmp = recHitDistFromWire(*recHit, layer);
1322     if (fabs(distTmp - simHitDist) < res) {
1323       res = fabs(distTmp - simHitDist);
1324       theBestRecHit = &(*recHit);
1325     }
1326   }  // End of loop over RecHits within the cell
1327 
1328   return theBestRecHit;
1329 }
1330 
1331 // Compute the distance from wire (cm) of a hits in a DTRecHit1DPair
1332 float GlobalRecHitsAnalyzer::recHitDistFromWire(const DTRecHit1DPair& hitPair, const DTLayer* layer) {
1333   // Compute the rechit distance from wire
1334   return fabs(hitPair.localPosition(DTEnums::Left).x() - hitPair.localPosition(DTEnums::Right).x()) / 2.;
1335 }
1336 
1337 // Compute the distance from wire (cm) of a hits in a DTRecHit1D
1338 float GlobalRecHitsAnalyzer::recHitDistFromWire(const DTRecHit1D& recHit, const DTLayer* layer) {
1339   return fabs(recHit.localPosition().x() - layer->specificTopology().wirePosition(recHit.wireId().wire()));
1340 }
1341 
1342 template <typename type>
1343 int GlobalRecHitsAnalyzer::compute(const DTGeometry* dtGeom,
1344                                    const std::map<DTWireId, std::vector<PSimHit>>& _simHitsPerWire,
1345                                    const std::map<DTWireId, std::vector<type>>& _recHitsPerWire,
1346                                    int step) {
1347   std::map<DTWireId, std::vector<PSimHit>> simHitsPerWire = _simHitsPerWire;
1348   std::map<DTWireId, std::vector<type>> recHitsPerWire = _recHitsPerWire;
1349   int nDt = 0;
1350   // Loop over cells with a muon SimHit
1351   for (std::map<DTWireId, std::vector<PSimHit>>::const_iterator wireAndSHits = simHitsPerWire.begin();
1352        wireAndSHits != simHitsPerWire.end();
1353        wireAndSHits++) {
1354     DTWireId wireId = (*wireAndSHits).first;
1355     std::vector<PSimHit> simHitsInCell = (*wireAndSHits).second;
1356 
1357     // Get the layer
1358     const DTLayer* layer = dtGeom->layer(wireId);
1359 
1360     // Look for a mu hit in the cell
1361     const PSimHit* muSimHit = DTHitQualityUtils::findMuSimHit(simHitsInCell);
1362     if (muSimHit == nullptr) {
1363       continue;  // Skip this cell
1364     }
1365 
1366     // Find the distance of the simhit from the wire
1367     float simHitWireDist = simHitDistFromWire(layer, wireId, *muSimHit);
1368     // Skip simhits out of the cell
1369     if (simHitWireDist > 2.1) {
1370       continue;  // Skip this cell
1371     }
1372 
1373     // Look for RecHits in the same cell
1374     if (recHitsPerWire.find(wireId) == recHitsPerWire.end()) {
1375       continue;  // No RecHit found in this cell
1376     } else {
1377       std::vector<type> recHits = recHitsPerWire[wireId];
1378 
1379       // Find the best RecHit
1380       const type* theBestRecHit = findBestRecHit(layer, wireId, recHits, simHitWireDist);
1381 
1382       float recHitWireDist = recHitDistFromWire(*theBestRecHit, layer);
1383 
1384       ++nDt;
1385 
1386       mehDtMuonRes->Fill(recHitWireDist - simHitWireDist);
1387 
1388     }  // find rechits
1389   }    // loop over simhits
1390 
1391   return nDt;
1392 }
1393 
1394 void GlobalRecHitsAnalyzer::plotResolution(const PSimHit& simHit,
1395                                            const CSCRecHit2D& recHit,
1396                                            const CSCLayer* layer,
1397                                            int chamberType) {
1398   GlobalPoint simHitPos = layer->toGlobal(simHit.localPosition());
1399   GlobalPoint recHitPos = layer->toGlobal(recHit.localPosition());
1400 
1401   mehCSCResRDPhi->Fill(recHitPos.phi() - simHitPos.phi());
1402 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!