Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-01-31 02:20:05 

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