Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-10-08 05:12:10

0001 /** \file GlobalDigisAnalyzer.cc
0002  *
0003  *  See header file for description of class
0004  *
0005  *  \author M. Strang SUNY-Buffalo
0006  */
0007 
0008 #include "CondFormats/EcalObjects/interface/EcalMGPAGainRatio.h"
0009 #include "DQMServices/Core/interface/DQMStore.h"
0010 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0011 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0012 #include "Validation/GlobalDigis/interface/GlobalDigisAnalyzer.h"
0013 
0014 GlobalDigisAnalyzer::GlobalDigisAnalyzer(const edm::ParameterSet &iPSet)
0015     : fName(""),
0016       verbosity(0),
0017       frequency(0),
0018       label(""),
0019       getAllProvenances(false),
0020       printProvenanceInfo(false),
0021       hitsProducer(""),
0022       theCSCStripPedestalSum(0),
0023       theCSCStripPedestalCount(0),
0024       count(0) {
0025   std::string MsgLoggerCat = "GlobalDigisAnalyzer_GlobalDigisAnalyzer";
0026 
0027   // get information from parameter set
0028   fName = iPSet.getUntrackedParameter<std::string>("Name");
0029   verbosity = iPSet.getUntrackedParameter<int>("Verbosity");
0030   frequency = iPSet.getUntrackedParameter<int>("Frequency");
0031   edm::ParameterSet m_Prov = iPSet.getParameter<edm::ParameterSet>("ProvenanceLookup");
0032   getAllProvenances = m_Prov.getUntrackedParameter<bool>("GetAllProvenances");
0033   printProvenanceInfo = m_Prov.getUntrackedParameter<bool>("PrintProvenanceInfo");
0034   hitsProducer = iPSet.getParameter<std::string>("hitsProducer");
0035 
0036   // get Labels to use to extract information
0037   ECalEBSrc_ = iPSet.getParameter<edm::InputTag>("ECalEBSrc");
0038   ECalEESrc_ = iPSet.getParameter<edm::InputTag>("ECalEESrc");
0039   ECalESSrc_ = iPSet.getParameter<edm::InputTag>("ECalESSrc");
0040   HCalSrc_ = iPSet.getParameter<edm::InputTag>("HCalSrc");
0041   HCalDigi_ = iPSet.getParameter<edm::InputTag>("HCalDigi");
0042   SiStripSrc_ = iPSet.getParameter<edm::InputTag>("SiStripSrc");
0043   SiPxlSrc_ = iPSet.getParameter<edm::InputTag>("SiPxlSrc");
0044   MuDTSrc_ = iPSet.getParameter<edm::InputTag>("MuDTSrc");
0045   MuCSCStripSrc_ = iPSet.getParameter<edm::InputTag>("MuCSCStripSrc");
0046   MuCSCWireSrc_ = iPSet.getParameter<edm::InputTag>("MuCSCWireSrc");
0047   MuRPCSrc_ = iPSet.getParameter<edm::InputTag>("MuRPCSrc");
0048 
0049   ECalEBSrc_Token_ = consumes<EBDigiCollection>(iPSet.getParameter<edm::InputTag>("ECalEBSrc"));
0050   ECalEESrc_Token_ = consumes<EEDigiCollection>(iPSet.getParameter<edm::InputTag>("ECalEESrc"));
0051   ECalESSrc_Token_ = consumes<ESDigiCollection>(iPSet.getParameter<edm::InputTag>("ECalESSrc"));
0052   HCalSrc_Token_ = consumes<edm::PCaloHitContainer>(iPSet.getParameter<edm::InputTag>("HCalSrc"));
0053   HBHEDigi_Token_ = consumes<edm::SortedCollection<HBHEDataFrame>>(iPSet.getParameter<edm::InputTag>("HCalDigi"));
0054   HODigi_Token_ = consumes<edm::SortedCollection<HODataFrame>>(iPSet.getParameter<edm::InputTag>("HCalDigi"));
0055   HFDigi_Token_ = consumes<edm::SortedCollection<HFDataFrame>>(iPSet.getParameter<edm::InputTag>("HCalDigi"));
0056   SiStripSrc_Token_ = consumes<edm::DetSetVector<SiStripDigi>>(iPSet.getParameter<edm::InputTag>("SiStripSrc"));
0057   SiPxlSrc_Token_ = consumes<edm::DetSetVector<PixelDigi>>(iPSet.getParameter<edm::InputTag>("SiPxlSrc"));
0058   MuDTSrc_Token_ = consumes<DTDigiCollection>(iPSet.getParameter<edm::InputTag>("MuDTSrc"));
0059   MuCSCStripSrc_Token_ = consumes<CSCStripDigiCollection>(iPSet.getParameter<edm::InputTag>("MuCSCStripSrc"));
0060   MuCSCWireSrc_Token_ = consumes<CSCWireDigiCollection>(iPSet.getParameter<edm::InputTag>("MuCSCWireSrc"));
0061   MuRPCSrc_Token_ = consumes<RPCDigiCollection>(iPSet.getParameter<edm::InputTag>("MuRPCSrc"));
0062   //
0063   const std::string barrelHitsName(hitsProducer + "EcalHitsEB");
0064   const std::string endcapHitsName(hitsProducer + "EcalHitsEE");
0065   const std::string preshowerHitsName(hitsProducer + "EcalHitsES");
0066   EBHits_Token_ = consumes<CrossingFrame<PCaloHit>>(edm::InputTag(std::string("mix"), std::string("barrelHitsName")));
0067   EEHits_Token_ = consumes<CrossingFrame<PCaloHit>>(edm::InputTag(std::string("mix"), std::string("endcapHitsName")));
0068   ESHits_Token_ =
0069       consumes<CrossingFrame<PCaloHit>>(edm::InputTag(std::string("mix"), std::string("preshowerHitsName")));
0070 
0071   RPCSimHit_Token_ =
0072       consumes<edm::PSimHitContainer>(edm::InputTag(std::string("g4SimHits"), std::string("MuonRPCHits")));
0073 
0074   ecalADCtoGevToken_ = esConsumes();
0075   rpcGeomToken_ = esConsumes();
0076   tTopoToken_ = esConsumes();
0077   hcaldbToken_ = esConsumes();
0078 
0079   // use value of first digit to determine default output level (inclusive)
0080   // 0 is none, 1 is basic, 2 is fill output, 3 is gather output
0081   verbosity %= 10;
0082 
0083   // print out Parameter Set information being used
0084   if (verbosity >= 0) {
0085     edm::LogInfo(MsgLoggerCat) << "\n===============================\n"
0086                                << "Initialized as EDAnalyzer with parameter values:\n"
0087                                << "    Name          = " << fName << "\n"
0088                                << "    Verbosity     = " << verbosity << "\n"
0089                                << "    Frequency     = " << frequency << "\n"
0090                                << "    GetProv       = " << getAllProvenances << "\n"
0091                                << "    PrintProv     = " << printProvenanceInfo << "\n"
0092                                << "    ECalEBSrc     = " << ECalEBSrc_.label() << ":" << ECalEBSrc_.instance() << "\n"
0093                                << "    ECalEESrc     = " << ECalEESrc_.label() << ":" << ECalEESrc_.instance() << "\n"
0094                                << "    ECalESSrc     = " << ECalESSrc_.label() << ":" << ECalESSrc_.instance() << "\n"
0095                                << "    HCalSrc       = " << HCalSrc_.label() << ":" << HCalSrc_.instance() << "\n"
0096                                << "    HCalDigi       = " << HCalDigi_.label() << ":" << HCalDigi_.instance() << "\n"
0097                                << "    SiStripSrc    = " << SiStripSrc_.label() << ":" << SiStripSrc_.instance() << "\n"
0098                                << "    SiPixelSrc    = " << SiPxlSrc_.label() << ":" << SiPxlSrc_.instance() << "\n"
0099                                << "    MuDTSrc       = " << MuDTSrc_.label() << ":" << MuDTSrc_.instance() << "\n"
0100                                << "    MuCSCStripSrc = " << MuCSCStripSrc_.label() << ":" << MuCSCStripSrc_.instance()
0101                                << "\n"
0102                                << "    MuCSCWireSrc  = " << MuCSCWireSrc_.label() << ":" << MuCSCWireSrc_.instance()
0103                                << "\n"
0104                                << "    MuRPCSrc      = " << MuRPCSrc_.label() << ":" << MuRPCSrc_.instance() << "\n"
0105                                << "===============================\n";
0106   }
0107 
0108   // set default constants
0109   // ECal
0110 
0111   ECalbarrelADCtoGeV_ = 0.035;
0112   ECalendcapADCtoGeV_ = 0.06;
0113 
0114   EcalMGPAGainRatio defaultRatios;
0115   ECalgainConv_[0] = 0.;
0116   ECalgainConv_[1] = 1.;
0117   ECalgainConv_[2] = defaultRatios.gain12Over6();
0118   ECalgainConv_[3] = ECalgainConv_[2] * (defaultRatios.gain6Over1());
0119 
0120   if (verbosity >= 0) {
0121     edm::LogInfo(MsgLoggerCat) << "Modified Calorimeter gain constants: g0 = " << ECalgainConv_[0]
0122                                << ", g1 = " << ECalgainConv_[1] << ", g2 = " << ECalgainConv_[2]
0123                                << ", g3 = " << ECalgainConv_[3];
0124   }
0125 }
0126 
0127 GlobalDigisAnalyzer::~GlobalDigisAnalyzer() {}
0128 
0129 void GlobalDigisAnalyzer::bookHistograms(DQMStore::IBooker &iBooker, edm::Run const &, edm::EventSetup const &) {
0130   // Si Strip
0131   std::string SiStripString[19] = {"TECW1",
0132                                    "TECW2",
0133                                    "TECW3",
0134                                    "TECW4",
0135                                    "TECW5",
0136                                    "TECW6",
0137                                    "TECW7",
0138                                    "TECW8",
0139                                    "TIBL1",
0140                                    "TIBL2",
0141                                    "TIBL3",
0142                                    "TIBL4",
0143                                    "TIDW1",
0144                                    "TIDW2",
0145                                    "TIDW3",
0146                                    "TOBL1",
0147                                    "TOBL2",
0148                                    "TOBL3",
0149                                    "TOBL4"};
0150   for (int i = 0; i < 19; ++i) {
0151     mehSiStripn[i] = nullptr;
0152     mehSiStripADC[i] = nullptr;
0153     mehSiStripStrip[i] = nullptr;
0154   }
0155   std::string hcharname, hchartitle;
0156   iBooker.setCurrentFolder("GlobalDigisV/SiStrips");
0157   for (int amend = 0; amend < 19; ++amend) {
0158     hcharname = "hSiStripn_" + SiStripString[amend];
0159     hchartitle = SiStripString[amend] + "  Digis";
0160     mehSiStripn[amend] = iBooker.book1D(hcharname, hchartitle, 5000, 0., 10000.);
0161     mehSiStripn[amend]->setAxisTitle("Number of Digis", 1);
0162     mehSiStripn[amend]->setAxisTitle("Count", 2);
0163 
0164     hcharname = "hSiStripADC_" + SiStripString[amend];
0165     hchartitle = SiStripString[amend] + " ADC";
0166     mehSiStripADC[amend] = iBooker.book1D(hcharname, hchartitle, 150, 0.0, 300.);
0167     mehSiStripADC[amend]->setAxisTitle("ADC", 1);
0168     mehSiStripADC[amend]->setAxisTitle("Count", 2);
0169 
0170     hcharname = "hSiStripStripADC_" + SiStripString[amend];
0171     hchartitle = SiStripString[amend] + " Strip";
0172     mehSiStripStrip[amend] = iBooker.book1D(hcharname, hchartitle, 200, 0.0, 800.);
0173     mehSiStripStrip[amend]->setAxisTitle("Strip Number", 1);
0174     mehSiStripStrip[amend]->setAxisTitle("Count", 2);
0175   }
0176 
0177   // HCal
0178   std::string HCalString[4] = {"HB", "HE", "HO", "HF"};
0179   float calnUpper[4] = {30000., 30000., 30000., 20000.};
0180   float calnLower[4] = {0., 0., 0., 0.};
0181   float SHEUpper[4] = {1., 1., 1., 1.};
0182   float SHEvAEEUpper[4] = {5000, 5000, 5000, 5000};
0183   float SHEvAEELower[4] = {-5000, -5000, -5000, -5000};
0184   int SHEvAEEnBins[4] = {200, 200, 200, 200};
0185   double ProfileUpper[4] = {1., 1., 1., 1.};
0186 
0187   for (int i = 0; i < 4; ++i) {
0188     mehHcaln[i] = nullptr;
0189     mehHcalAEE[i] = nullptr;
0190     mehHcalSHE[i] = nullptr;
0191     mehHcalAEESHE[i] = nullptr;
0192     mehHcalSHEvAEE[i] = nullptr;
0193   }
0194   iBooker.setCurrentFolder("GlobalDigisV/HCals");
0195 
0196   for (int amend = 0; amend < 4; ++amend) {
0197     hcharname = "hHcaln_" + HCalString[amend];
0198     hchartitle = HCalString[amend] + "  digis";
0199     mehHcaln[amend] = iBooker.book1D(hcharname, hchartitle, 10000, calnLower[amend], calnUpper[amend]);
0200     mehHcaln[amend]->setAxisTitle("Number of Digis", 1);
0201     mehHcaln[amend]->setAxisTitle("Count", 2);
0202 
0203     hcharname = "hHcalAEE_" + HCalString[amend];
0204     hchartitle = HCalString[amend] + "Cal AEE";
0205     mehHcalAEE[amend] = iBooker.book1D(hcharname, hchartitle, 60, -10., 50.);
0206     mehHcalAEE[amend]->setAxisTitle("Analog Equivalent Energy", 1);
0207     mehHcalAEE[amend]->setAxisTitle("Count", 2);
0208 
0209     hcharname = "hHcalSHE_" + HCalString[amend];
0210     hchartitle = HCalString[amend] + "Cal SHE";
0211     mehHcalSHE[amend] = iBooker.book1D(hcharname, hchartitle, 1000, 0.0, SHEUpper[amend]);
0212     mehHcalSHE[amend]->setAxisTitle("Simulated Hit Energy", 1);
0213     mehHcalSHE[amend]->setAxisTitle("Count", 2);
0214 
0215     hcharname = "hHcalAEESHE_" + HCalString[amend];
0216     hchartitle = HCalString[amend] + "Cal AEE/SHE";
0217     mehHcalAEESHE[amend] =
0218         iBooker.book1D(hcharname, hchartitle, SHEvAEEnBins[amend], SHEvAEELower[amend], SHEvAEEUpper[amend]);
0219     mehHcalAEESHE[amend]->setAxisTitle("ADC / SHE", 1);
0220     mehHcalAEESHE[amend]->setAxisTitle("Count", 2);
0221 
0222     hcharname = "hHcalSHEvAEE_" + HCalString[amend];
0223     hchartitle = HCalString[amend] + "Cal SHE vs. AEE";
0224     mehHcalSHEvAEE[amend] =
0225         iBooker.bookProfile(hcharname, hchartitle, 60, -10., 50., 100, 0., (float)ProfileUpper[amend], "");
0226     mehHcalSHEvAEE[amend]->setAxisTitle("AEE / SHE", 1);
0227     mehHcalSHEvAEE[amend]->setAxisTitle("SHE", 2);
0228   }
0229 
0230   // Ecal
0231   std::string ECalString[2] = {"EB", "EE"};
0232 
0233   for (int i = 0; i < 2; ++i) {
0234     mehEcaln[i] = nullptr;
0235     mehEcalAEE[i] = nullptr;
0236     mehEcalSHE[i] = nullptr;
0237     mehEcalMaxPos[i] = nullptr;
0238     mehEcalMultvAEE[i] = nullptr;
0239     mehEcalSHEvAEESHE[i] = nullptr;
0240   }
0241   iBooker.setCurrentFolder("GlobalDigisV/ECals");
0242 
0243   for (int amend = 0; amend < 2; ++amend) {
0244     hcharname = "hEcaln_" + ECalString[amend];
0245     hchartitle = ECalString[amend] + "  digis";
0246     mehEcaln[amend] = iBooker.book1D(hcharname, hchartitle, 3000, 0., 40000.);
0247     mehEcaln[amend]->setAxisTitle("Number of Digis", 1);
0248     mehEcaln[amend]->setAxisTitle("Count", 2);
0249 
0250     hcharname = "hEcalAEE_" + ECalString[amend];
0251     hchartitle = ECalString[amend] + "Cal AEE";
0252     mehEcalAEE[amend] = iBooker.book1D(hcharname, hchartitle, 1000, 0., 100.);
0253     mehEcalAEE[amend]->setAxisTitle("Analog Equivalent Energy", 1);
0254     mehEcalAEE[amend]->setAxisTitle("Count", 2);
0255 
0256     hcharname = "hEcalSHE_" + ECalString[amend];
0257     hchartitle = ECalString[amend] + "Cal SHE";
0258     mehEcalSHE[amend] = iBooker.book1D(hcharname, hchartitle, 500, 0., 50.);
0259     mehEcalSHE[amend]->setAxisTitle("Simulated Hit Energy", 1);
0260     mehEcalSHE[amend]->setAxisTitle("Count", 2);
0261 
0262     hcharname = "hEcalMaxPos_" + ECalString[amend];
0263     hchartitle = ECalString[amend] + "Cal MaxPos";
0264     mehEcalMaxPos[amend] = iBooker.book1D(hcharname, hchartitle, 10, 0., 10.);
0265     mehEcalMaxPos[amend]->setAxisTitle("Maximum Position", 1);
0266     mehEcalMaxPos[amend]->setAxisTitle("Count", 2);
0267 
0268     hcharname = "hEcalSHEvAEESHE_" + ECalString[amend];
0269     hchartitle = ECalString[amend] + "Cal SHE vs. AEE/SHE";
0270     mehEcalSHEvAEESHE[amend] = iBooker.bookProfile(hcharname, hchartitle, 1000, 0., 100., 500, 0., 50., "");
0271     mehEcalSHEvAEESHE[amend]->setAxisTitle("AEE / SHE", 1);
0272     mehEcalSHEvAEESHE[amend]->setAxisTitle("SHE", 2);
0273 
0274     hcharname = "hEcalMultvAEE_" + ECalString[amend];
0275     hchartitle = ECalString[amend] + "Cal Multi vs. AEE";
0276     mehEcalMultvAEE[amend] = iBooker.bookProfile(hcharname, hchartitle, 1000, 0., 100., 4000, 0., 40000., "");
0277     mehEcalMultvAEE[amend]->setAxisTitle("Analog Equivalent Energy", 1);
0278     mehEcalMultvAEE[amend]->setAxisTitle("Number of Digis", 2);
0279   }
0280   mehEScaln = nullptr;
0281 
0282   hcharname = "hEcaln_ES";
0283   hchartitle = "ESCAL  digis";
0284   mehEScaln = iBooker.book1D(hcharname, hchartitle, 1000, 0., 5000.);
0285   mehEScaln->setAxisTitle("Number of Digis", 1);
0286   mehEScaln->setAxisTitle("Count", 2);
0287 
0288   std::string ADCNumber[3] = {"0", "1", "2"};
0289   for (int i = 0; i < 3; ++i) {
0290     mehEScalADC[i] = nullptr;
0291     hcharname = "hEcalADC" + ADCNumber[i] + "_ES";
0292     hchartitle = "ESCAL  ADC" + ADCNumber[i];
0293     mehEScalADC[i] = iBooker.book1D(hcharname, hchartitle, 1500, 0., 1500.);
0294     mehEScalADC[i]->setAxisTitle("ADC" + ADCNumber[i], 1);
0295     mehEScalADC[i]->setAxisTitle("Count", 2);
0296   }
0297 
0298   // Si Pixels ***DONE***
0299   std::string SiPixelString[7] = {"BRL1", "BRL2", "BRL3", "FWD1n", "FWD1p", "FWD2n", "FWD2p"};
0300   for (int j = 0; j < 7; ++j) {
0301     mehSiPixeln[j] = nullptr;
0302     mehSiPixelADC[j] = nullptr;
0303     mehSiPixelRow[j] = nullptr;
0304     mehSiPixelCol[j] = nullptr;
0305   }
0306 
0307   iBooker.setCurrentFolder("GlobalDigisV/SiPixels");
0308   for (int amend = 0; amend < 7; ++amend) {
0309     hcharname = "hSiPixeln_" + SiPixelString[amend];
0310     hchartitle = SiPixelString[amend] + " Digis";
0311     if (amend < 3)
0312       mehSiPixeln[amend] = iBooker.book1D(hcharname, hchartitle, 500, 0., 1000.);
0313     else
0314       mehSiPixeln[amend] = iBooker.book1D(hcharname, hchartitle, 500, 0., 1000.);
0315     mehSiPixeln[amend]->setAxisTitle("Number of Digis", 1);
0316     mehSiPixeln[amend]->setAxisTitle("Count", 2);
0317 
0318     hcharname = "hSiPixelADC_" + SiPixelString[amend];
0319     hchartitle = SiPixelString[amend] + " ADC";
0320     mehSiPixelADC[amend] = iBooker.book1D(hcharname, hchartitle, 150, 0.0, 300.);
0321     mehSiPixelADC[amend]->setAxisTitle("ADC", 1);
0322     mehSiPixelADC[amend]->setAxisTitle("Count", 2);
0323 
0324     hcharname = "hSiPixelRow_" + SiPixelString[amend];
0325     hchartitle = SiPixelString[amend] + " Row";
0326     mehSiPixelRow[amend] = iBooker.book1D(hcharname, hchartitle, 100, 0.0, 100.);
0327     mehSiPixelRow[amend]->setAxisTitle("Row Number", 1);
0328     mehSiPixelRow[amend]->setAxisTitle("Count", 2);
0329 
0330     hcharname = "hSiPixelColumn_" + SiPixelString[amend];
0331     hchartitle = SiPixelString[amend] + " Column";
0332     mehSiPixelCol[amend] = iBooker.book1D(hcharname, hchartitle, 200, 0.0, 500.);
0333     mehSiPixelCol[amend]->setAxisTitle("Column Number", 1);
0334     mehSiPixelCol[amend]->setAxisTitle("Count", 2);
0335   }
0336 
0337   // Muons
0338   iBooker.setCurrentFolder("GlobalDigisV/Muons");
0339 
0340   // DT
0341   std::string MuonString[4] = {"MB1", "MB2", "MB3", "MB4"};
0342 
0343   for (int i = 0; i < 4; ++i) {
0344     mehDtMuonn[i] = nullptr;
0345     mehDtMuonLayer[i] = nullptr;
0346     mehDtMuonTime[i] = nullptr;
0347     mehDtMuonTimevLayer[i] = nullptr;
0348   }
0349 
0350   for (int j = 0; j < 4; ++j) {
0351     hcharname = "hDtMuonn_" + MuonString[j];
0352     hchartitle = MuonString[j] + "  digis";
0353     mehDtMuonn[j] = iBooker.book1D(hcharname, hchartitle, 250, 0., 500.);
0354     mehDtMuonn[j]->setAxisTitle("Number of Digis", 1);
0355     mehDtMuonn[j]->setAxisTitle("Count", 2);
0356 
0357     hcharname = "hDtLayer_" + MuonString[j];
0358     hchartitle = MuonString[j] + "  Layer";
0359     mehDtMuonLayer[j] = iBooker.book1D(hcharname, hchartitle, 12, 1., 13.);
0360     mehDtMuonLayer[j]->setAxisTitle("4 * (SuperLayer - 1) + Layer", 1);
0361     mehDtMuonLayer[j]->setAxisTitle("Count", 2);
0362 
0363     hcharname = "hDtMuonTime_" + MuonString[j];
0364     hchartitle = MuonString[j] + "  Time";
0365     mehDtMuonTime[j] = iBooker.book1D(hcharname, hchartitle, 300, 400., 1000.);
0366     mehDtMuonTime[j]->setAxisTitle("Time", 1);
0367     mehDtMuonTime[j]->setAxisTitle("Count", 2);
0368 
0369     hcharname = "hDtMuonTimevLayer_" + MuonString[j];
0370     hchartitle = MuonString[j] + "  Time vs. Layer";
0371     mehDtMuonTimevLayer[j] = iBooker.bookProfile(hcharname, hchartitle, 12, 1., 13., 300, 400., 1000., "");
0372     mehDtMuonTimevLayer[j]->setAxisTitle("4 * (SuperLayer - 1) + Layer", 1);
0373     mehDtMuonTimevLayer[j]->setAxisTitle("Time", 2);
0374   }
0375 
0376   // CSC
0377   mehCSCStripn = nullptr;
0378   hcharname = "hCSCStripn";
0379   hchartitle = "CSC Strip digis";
0380   mehCSCStripn = iBooker.book1D(hcharname, hchartitle, 250, 0., 500.);
0381   mehCSCStripn->setAxisTitle("Number of Digis", 1);
0382   mehCSCStripn->setAxisTitle("Count", 2);
0383 
0384   mehCSCStripADC = nullptr;
0385   hcharname = "hCSCStripADC";
0386   hchartitle = "CSC Strip ADC";
0387   mehCSCStripADC = iBooker.book1D(hcharname, hchartitle, 110, 0., 1100.);
0388   mehCSCStripADC->setAxisTitle("ADC", 1);
0389   mehCSCStripADC->setAxisTitle("Count", 2);
0390 
0391   mehCSCWiren = nullptr;
0392   hcharname = "hCSCWiren";
0393   hchartitle = "CSC Wire digis";
0394   mehCSCWiren = iBooker.book1D(hcharname, hchartitle, 250, 0., 500.);
0395   mehCSCWiren->setAxisTitle("Number of Digis", 1);
0396   mehCSCWiren->setAxisTitle("Count", 2);
0397 
0398   mehCSCWireTime = nullptr;
0399   hcharname = "hCSCWireTime";
0400   hchartitle = "CSC Wire Time";
0401   mehCSCWireTime = iBooker.book1D(hcharname, hchartitle, 10, 0., 10.);
0402   mehCSCWireTime->setAxisTitle("Time", 1);
0403   mehCSCWireTime->setAxisTitle("Count", 2);
0404 
0405   // RPC
0406   mehRPCMuonn = nullptr;
0407   hcharname = "hRPCMuonn";
0408   hchartitle = "RPC digis";
0409   mehCSCStripn = iBooker.book1D(hcharname, hchartitle, 250, 0., 500.);
0410   mehCSCStripn->setAxisTitle("Number of Digis", 1);
0411   mehCSCStripn->setAxisTitle("Count", 2);
0412 
0413   std::string MuonRPCString[5] = {"Wmin2", "Wmin1", "W0", "Wpu1", "Wpu2"};
0414   for (int i = 0; i < 5; ++i) {
0415     mehRPCRes[i] = nullptr;
0416   }
0417 
0418   for (int j = 0; j < 5; ++j) {
0419     hcharname = "hRPCRes_" + MuonRPCString[j];
0420     hchartitle = MuonRPCString[j] + "  Digi - Sim";
0421     mehRPCRes[j] = iBooker.book1D(hcharname, hchartitle, 200, -8., 8.);
0422     mehRPCRes[j]->setAxisTitle("Digi - Sim center of strip x", 1);
0423     mehRPCRes[j]->setAxisTitle("Count", 2);
0424   }
0425 }
0426 
0427 void GlobalDigisAnalyzer::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0428   std::string MsgLoggerCat = "GlobalDigisAnalyzer_analyze";
0429 
0430   // keep track of number of events processed
0431   ++count;
0432 
0433   // THIS BLOCK MIGRATED HERE FROM beginJob:
0434   // setup calorimeter constants from service
0435   const EcalADCToGeVConstant *agc = &iSetup.getData(ecalADCtoGevToken_);
0436   ECalbarrelADCtoGeV_ = agc->getEBValue();
0437   ECalendcapADCtoGeV_ = agc->getEEValue();
0438   if (verbosity >= 0) {
0439     edm::LogInfo(MsgLoggerCat) << "Modified Calorimeter ADCtoGeV constants: barrel = " << ECalbarrelADCtoGeV_
0440                                << ", endcap = " << ECalendcapADCtoGeV_;
0441   }
0442 
0443   // get event id information
0444   edm::RunNumber_t nrun = iEvent.id().run();
0445   edm::EventNumber_t nevt = iEvent.id().event();
0446 
0447   if (verbosity > 0) {
0448     edm::LogInfo(MsgLoggerCat) << "Processing run " << nrun << ", event " << nevt << " (" << count << " events total)";
0449   } else if (verbosity == 0) {
0450     if (nevt % frequency == 0 || nevt == 1) {
0451       edm::LogInfo(MsgLoggerCat) << "Processing run " << nrun << ", event " << nevt << " (" << count
0452                                  << " events total)";
0453     }
0454   }
0455 
0456   // look at information available in the event
0457   if (getAllProvenances) {
0458     std::vector<const edm::StableProvenance *> AllProv;
0459     iEvent.getAllStableProvenance(AllProv);
0460 
0461     if (verbosity >= 0)
0462       edm::LogInfo(MsgLoggerCat) << "Number of Provenances = " << AllProv.size();
0463 
0464     if (printProvenanceInfo && (verbosity >= 0)) {
0465       TString eventout("\nProvenance info:\n");
0466 
0467       for (unsigned int i = 0; i < AllProv.size(); ++i) {
0468         eventout += "\n       ******************************";
0469         eventout += "\n       Module       : ";
0470         eventout += AllProv[i]->moduleLabel();
0471         eventout += "\n       ProductID    : ";
0472         eventout += AllProv[i]->productID().id();
0473         eventout += "\n       ClassName    : ";
0474         eventout += AllProv[i]->className();
0475         eventout += "\n       InstanceName : ";
0476         eventout += AllProv[i]->productInstanceName();
0477         eventout += "\n       BranchName   : ";
0478         eventout += AllProv[i]->branchName();
0479       }
0480       eventout += "\n       ******************************\n";
0481       edm::LogInfo(MsgLoggerCat) << eventout << "\n";
0482       printProvenanceInfo = false;
0483     }
0484     getAllProvenances = false;
0485   }
0486 
0487   // call fill functions
0488   // gather Ecal information from event
0489   fillECal(iEvent, iSetup);
0490   // gather Hcal information from event
0491   fillHCal(iEvent, iSetup);
0492   // gather Track information from event
0493   fillTrk(iEvent, iSetup);
0494   // gather Muon information from event
0495   fillMuon(iEvent, iSetup);
0496 
0497   if (verbosity > 0)
0498     edm::LogInfo(MsgLoggerCat) << "Done gathering data from event.";
0499 
0500   if (verbosity > 2)
0501     edm::LogInfo(MsgLoggerCat) << "Saving event contents:";
0502 
0503   return;
0504 }
0505 
0506 void GlobalDigisAnalyzer::fillECal(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0507   std::string MsgLoggerCat = "GlobalDigisAnalyzer_fillECal";
0508 
0509   TString eventout;
0510   if (verbosity > 0)
0511     eventout = "\nGathering info:";
0512 
0513   // extract crossing frame from event
0514   edm::Handle<CrossingFrame<PCaloHit>> crossingFrame;
0515 
0516   ////////////////////////
0517   // extract EB information
0518   ////////////////////////
0519   bool isBarrel = true;
0520   edm::Handle<EBDigiCollection> EcalDigiEB;
0521   iEvent.getByToken(ECalEBSrc_Token_, EcalDigiEB);
0522   bool validDigiEB = true;
0523   if (!EcalDigiEB.isValid()) {
0524     LogDebug(MsgLoggerCat) << "Unable to find EcalDigiEB in event!";
0525     validDigiEB = false;
0526   }
0527   if (validDigiEB) {
0528     if (EcalDigiEB->empty())
0529       isBarrel = false;
0530 
0531     if (isBarrel) {
0532       // loop over simhits
0533       MapType ebSimMap;
0534       iEvent.getByToken(EBHits_Token_, crossingFrame);
0535       bool validXFrame = true;
0536       if (!crossingFrame.isValid()) {
0537         LogDebug(MsgLoggerCat) << "Unable to find cal barrel crossingFrame in event!";
0538         validXFrame = false;
0539       }
0540       if (validXFrame) {
0541         const MixCollection<PCaloHit> barrelHits(crossingFrame.product());
0542 
0543         // keep track of sum of simhit energy in each crystal
0544         for (auto const &iHit : barrelHits) {
0545           EBDetId ebid = EBDetId(iHit.id());
0546 
0547           uint32_t crystid = ebid.rawId();
0548           ebSimMap[crystid] += iHit.energy();
0549         }
0550       }
0551 
0552       // loop over digis
0553       const EBDigiCollection *barrelDigi = EcalDigiEB.product();
0554 
0555       std::vector<double> ebAnalogSignal;
0556       std::vector<double> ebADCCounts;
0557       std::vector<double> ebADCGains;
0558       ebAnalogSignal.reserve(EBDataFrame::MAXSAMPLES);
0559       ebADCCounts.reserve(EBDataFrame::MAXSAMPLES);
0560       ebADCGains.reserve(EBDataFrame::MAXSAMPLES);
0561 
0562       int i = 0;
0563       for (unsigned int digis = 0; digis < EcalDigiEB->size(); ++digis) {
0564         ++i;
0565 
0566         EBDataFrame ebdf = (*barrelDigi)[digis];
0567         int nrSamples = ebdf.size();
0568 
0569         EBDetId ebid = ebdf.id();
0570 
0571         double Emax = 0;
0572         int Pmax = 0;
0573         double pedestalPreSampleAnalog = 0.;
0574 
0575         for (int sample = 0; sample < nrSamples; ++sample) {
0576           ebAnalogSignal[sample] = 0.;
0577           ebADCCounts[sample] = 0.;
0578           ebADCGains[sample] = -1.;
0579         }
0580 
0581         // calculate maximum energy and pedestal
0582         for (int sample = 0; sample < nrSamples; ++sample) {
0583           EcalMGPASample thisSample = ebdf[sample];
0584           ebADCCounts[sample] = (thisSample.adc());
0585           ebADCGains[sample] = (thisSample.gainId());
0586           ebAnalogSignal[sample] = (ebADCCounts[sample] * ECalgainConv_[(int)ebADCGains[sample]] * ECalbarrelADCtoGeV_);
0587           if (Emax < ebAnalogSignal[sample]) {
0588             Emax = ebAnalogSignal[sample];
0589             Pmax = sample;
0590           }
0591           if (sample < 3) {
0592             pedestalPreSampleAnalog +=
0593                 ebADCCounts[sample] * ECalgainConv_[(int)ebADCGains[sample]] * ECalbarrelADCtoGeV_;
0594           }
0595         }
0596         pedestalPreSampleAnalog /= 3.;
0597 
0598         // calculate pedestal subtracted digi energy in the crystal
0599         double Erec = Emax - pedestalPreSampleAnalog * ECalgainConv_[(int)ebADCGains[Pmax]];
0600 
0601         // gather necessary information
0602         mehEcalMaxPos[0]->Fill(Pmax);
0603         mehEcalSHE[0]->Fill(ebSimMap[ebid.rawId()]);
0604         mehEcalAEE[0]->Fill(Erec);
0605         // Adding protection against FPE
0606         if (ebSimMap[ebid.rawId()] != 0) {
0607           mehEcalSHEvAEESHE[0]->Fill(Erec / ebSimMap[ebid.rawId()], ebSimMap[ebid.rawId()]);
0608         }
0609         // else {
0610         //  std::cout<<"Would have been an FPE! with
0611         //  ebSimMap[ebid.rawId()]==0\n";
0612         //}
0613 
0614         mehEcalMultvAEE[0]->Fill(Pmax, (float)i);
0615       }
0616 
0617       if (verbosity > 1) {
0618         eventout += "\n          Number of EBDigis collected:.............. ";
0619         eventout += i;
0620       }
0621       mehEcaln[0]->Fill((float)i);
0622     }
0623   }
0624 
0625   /////////////////////////
0626   // extract EE information
0627   ////////////////////////
0628   bool isEndCap = true;
0629   edm::Handle<EEDigiCollection> EcalDigiEE;
0630   iEvent.getByToken(ECalEESrc_Token_, EcalDigiEE);
0631   bool validDigiEE = true;
0632   if (!EcalDigiEE.isValid()) {
0633     LogDebug(MsgLoggerCat) << "Unable to find EcalDigiEE in event!";
0634     validDigiEE = false;
0635   }
0636   if (validDigiEE) {
0637     if (EcalDigiEE->empty())
0638       isEndCap = false;
0639 
0640     if (isEndCap) {
0641       // loop over simhits
0642       MapType eeSimMap;
0643       iEvent.getByToken(EEHits_Token_, crossingFrame);
0644       bool validXFrame = true;
0645       if (!crossingFrame.isValid()) {
0646         LogDebug(MsgLoggerCat) << "Unable to find cal endcap crossingFrame in event!";
0647         validXFrame = false;
0648       }
0649       if (validXFrame) {
0650         const MixCollection<PCaloHit> endcapHits(crossingFrame.product());
0651 
0652         // keep track of sum of simhit energy in each crystal
0653         for (auto const &iHit : endcapHits) {
0654           EEDetId eeid = EEDetId(iHit.id());
0655 
0656           uint32_t crystid = eeid.rawId();
0657           eeSimMap[crystid] += iHit.energy();
0658         }
0659       }
0660 
0661       // loop over digis
0662       const EEDigiCollection *endcapDigi = EcalDigiEE.product();
0663 
0664       std::vector<double> eeAnalogSignal;
0665       std::vector<double> eeADCCounts;
0666       std::vector<double> eeADCGains;
0667       eeAnalogSignal.reserve(EEDataFrame::MAXSAMPLES);
0668       eeADCCounts.reserve(EEDataFrame::MAXSAMPLES);
0669       eeADCGains.reserve(EEDataFrame::MAXSAMPLES);
0670 
0671       int inc = 0;
0672       for (unsigned int digis = 0; digis < EcalDigiEE->size(); ++digis) {
0673         ++inc;
0674 
0675         EEDataFrame eedf = (*endcapDigi)[digis];
0676         int nrSamples = eedf.size();
0677 
0678         EEDetId eeid = eedf.id();
0679 
0680         double Emax = 0;
0681         int Pmax = 0;
0682         double pedestalPreSampleAnalog = 0.;
0683 
0684         for (int sample = 0; sample < nrSamples; ++sample) {
0685           eeAnalogSignal[sample] = 0.;
0686           eeADCCounts[sample] = 0.;
0687           eeADCGains[sample] = -1.;
0688         }
0689 
0690         // calculate maximum enery and pedestal
0691         for (int sample = 0; sample < nrSamples; ++sample) {
0692           EcalMGPASample thisSample = eedf[sample];
0693 
0694           eeADCCounts[sample] = (thisSample.adc());
0695           eeADCGains[sample] = (thisSample.gainId());
0696           eeAnalogSignal[sample] = (eeADCCounts[sample] * ECalgainConv_[(int)eeADCGains[sample]] * ECalbarrelADCtoGeV_);
0697           if (Emax < eeAnalogSignal[sample]) {
0698             Emax = eeAnalogSignal[sample];
0699             Pmax = sample;
0700           }
0701           if (sample < 3) {
0702             pedestalPreSampleAnalog +=
0703                 eeADCCounts[sample] * ECalgainConv_[(int)eeADCGains[sample]] * ECalbarrelADCtoGeV_;
0704           }
0705         }
0706         pedestalPreSampleAnalog /= 3.;
0707 
0708         // calculate pedestal subtracted digi energy in the crystal
0709         double Erec = Emax - pedestalPreSampleAnalog * ECalgainConv_[(int)eeADCGains[Pmax]];
0710 
0711         // gather necessary information
0712         mehEcalMaxPos[1]->Fill(Pmax);
0713         mehEcalSHE[1]->Fill(eeSimMap[eeid.rawId()]);
0714         mehEcalAEE[1]->Fill(Erec);
0715         // Adding protection against FPE
0716         if (eeSimMap[eeid.rawId()] != 0) {
0717           mehEcalSHEvAEESHE[1]->Fill(Erec / eeSimMap[eeid.rawId()], eeSimMap[eeid.rawId()]);
0718         }
0719         // else{
0720         //  std::cout<<"Would have been an FPE! with
0721         //  eeSimMap[eeid.rawId()]==0\n";
0722         //}
0723         mehEcalMultvAEE[1]->Fill(Pmax, (float)inc);
0724       }
0725 
0726       if (verbosity > 1) {
0727         eventout += "\n          Number of EEDigis collected:.............. ";
0728         eventout += inc;
0729       }
0730 
0731       mehEcaln[1]->Fill((float)inc);
0732     }
0733   }
0734 
0735   /////////////////////////
0736   // extract ES information
0737   ////////////////////////
0738   bool isPreshower = true;
0739   edm::Handle<ESDigiCollection> EcalDigiES;
0740   iEvent.getByToken(ECalESSrc_Token_, EcalDigiES);
0741   bool validDigiES = true;
0742   if (!EcalDigiES.isValid()) {
0743     LogDebug(MsgLoggerCat) << "Unable to find EcalDigiES in event!";
0744     [[clang::suppress]] validDigiES = false;
0745   }
0746 
0747   // ONLY WHILE GEOMETRY IS REMOVED
0748   validDigiES = false;
0749 
0750   if (validDigiES) {
0751     if (EcalDigiES->empty())
0752       isPreshower = false;
0753 
0754     if (isPreshower) {
0755       // loop over simhits
0756       iEvent.getByToken(ESHits_Token_, crossingFrame);
0757       bool validXFrame = true;
0758       if (!crossingFrame.isValid()) {
0759         LogDebug(MsgLoggerCat) << "Unable to find cal preshower crossingFrame in event!";
0760         validXFrame = false;
0761       }
0762       if (validXFrame) {
0763         const MixCollection<PCaloHit> preshowerHits(crossingFrame.product());
0764 
0765         // keep track of sum of simhit energy in each crystal
0766         MapType esSimMap;
0767         for (auto const &iHit : preshowerHits) {
0768           ESDetId esid = ESDetId(iHit.id());
0769 
0770           uint32_t crystid = esid.rawId();
0771           esSimMap[crystid] += iHit.energy();
0772         }
0773       }
0774 
0775       // loop over digis
0776       const ESDigiCollection *preshowerDigi = EcalDigiES.product();
0777 
0778       std::vector<double> esADCCounts;
0779       esADCCounts.reserve(ESDataFrame::MAXSAMPLES);
0780 
0781       int i = 0;
0782       for (unsigned int digis = 0; digis < EcalDigiES->size(); ++digis) {
0783         ++i;
0784 
0785         ESDataFrame esdf = (*preshowerDigi)[digis];
0786         int nrSamples = esdf.size();
0787 
0788         for (int sample = 0; sample < nrSamples; ++sample) {
0789           esADCCounts[sample] = 0.;
0790         }
0791 
0792         // gether ADC counts
0793         for (int sample = 0; sample < nrSamples; ++sample) {
0794           ESSample thisSample = esdf[sample];
0795           esADCCounts[sample] = (thisSample.adc());
0796         }
0797 
0798         mehEScalADC[0]->Fill(esADCCounts[0]);
0799         mehEScalADC[1]->Fill(esADCCounts[1]);
0800         mehEScalADC[2]->Fill(esADCCounts[2]);
0801       }
0802 
0803       if (verbosity > 1) {
0804         eventout += "\n          Number of ESDigis collected:.............. ";
0805         eventout += i;
0806       }
0807 
0808       mehEScaln->Fill((float)i);
0809     }
0810   }
0811   if (verbosity > 0)
0812     edm::LogInfo(MsgLoggerCat) << eventout << "\n";
0813 
0814   return;
0815 }
0816 
0817 void GlobalDigisAnalyzer::fillHCal(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0818   std::string MsgLoggerCat = "GlobalDigisAnalyzer_fillHCal";
0819 
0820   TString eventout;
0821   if (verbosity > 0)
0822     eventout = "\nGathering info:";
0823 
0824   // get calibration info
0825   const auto &HCalconditions = iSetup.getHandle(hcaldbToken_);
0826   if (!HCalconditions.isValid()) {
0827     edm::LogWarning(MsgLoggerCat) << "Unable to find HCalconditions in event!";
0828     return;
0829   }
0830   // HcalCalibrations calibrations;
0831   CaloSamples tool;
0832 
0833   ///////////////////////
0834   // extract simhit info
0835   //////////////////////
0836   edm::Handle<edm::PCaloHitContainer> hcalHits;
0837   iEvent.getByToken(HCalSrc_Token_, hcalHits);
0838   bool validhcalHits = true;
0839   if (!hcalHits.isValid()) {
0840     LogDebug(MsgLoggerCat) << "Unable to find hcalHits in event!";
0841     validhcalHits = false;
0842   }
0843   MapType fHBEnergySimHits;
0844   MapType fHEEnergySimHits;
0845   MapType fHOEnergySimHits;
0846   MapType fHFEnergySimHits;
0847   if (validhcalHits) {
0848     const edm::PCaloHitContainer *simhitResult = hcalHits.product();
0849 
0850     for (std::vector<PCaloHit>::const_iterator simhits = simhitResult->begin(); simhits != simhitResult->end();
0851          ++simhits) {
0852       HcalDetId detId(simhits->id());
0853       uint32_t cellid = detId.rawId();
0854 
0855       if (detId.subdet() == sdHcalBrl) {
0856         fHBEnergySimHits[cellid] += simhits->energy();
0857       }
0858       if (detId.subdet() == sdHcalEC) {
0859         fHEEnergySimHits[cellid] += simhits->energy();
0860       }
0861       if (detId.subdet() == sdHcalOut) {
0862         fHOEnergySimHits[cellid] += simhits->energy();
0863       }
0864       if (detId.subdet() == sdHcalFwd) {
0865         fHFEnergySimHits[cellid] += simhits->energy();
0866       }
0867     }
0868   }
0869 
0870   ////////////////////////
0871   // get HBHE information
0872   ///////////////////////
0873   edm::Handle<edm::SortedCollection<HBHEDataFrame>> hbhe;
0874   iEvent.getByToken(HBHEDigi_Token_, hbhe);
0875   bool validHBHE = true;
0876   if (!hbhe.isValid()) {
0877     LogDebug(MsgLoggerCat) << "Unable to find HBHEDataFrame in event!";
0878     validHBHE = false;
0879   }
0880 
0881   if (validHBHE) {
0882     edm::SortedCollection<HBHEDataFrame>::const_iterator ihbhe;
0883 
0884     int iHB = 0;
0885     int iHE = 0;
0886     for (ihbhe = hbhe->begin(); ihbhe != hbhe->end(); ++ihbhe) {
0887       HcalDetId cell(ihbhe->id());
0888 
0889       if ((cell.subdet() == sdHcalBrl) || (cell.subdet() == sdHcalEC)) {
0890         // HCalconditions->makeHcalCalibration(cell, &calibrations);
0891         const HcalCalibrations &calibrations = HCalconditions->getHcalCalibrations(cell);
0892         const HcalQIECoder *channelCoder = HCalconditions->getHcalCoder(cell);
0893         const HcalQIEShape *shape = HCalconditions->getHcalShape(channelCoder);
0894         HcalCoderDb coder(*channelCoder, *shape);
0895         coder.adc2fC(*ihbhe, tool);
0896 
0897         // get HB info
0898         if (cell.subdet() == sdHcalBrl) {
0899           ++iHB;
0900           float fDigiSum = 0.0;
0901           for (int ii = 0; ii < tool.size(); ++ii) {
0902             // default ped is 4.5
0903             int capid = (*ihbhe)[ii].capid();
0904             fDigiSum += (tool[ii] - calibrations.pedestal(capid));
0905           }
0906 
0907           mehHcalSHE[0]->Fill(fHFEnergySimHits[cell.rawId()]);
0908           mehHcalAEE[0]->Fill(fDigiSum);
0909           // Adding protection against FPE
0910           if (fHFEnergySimHits[cell.rawId()] != 0) {
0911             mehHcalAEESHE[0]->Fill(fDigiSum / fHFEnergySimHits[cell.rawId()]);
0912           }
0913           // else {
0914           //  std::cout<<"It would have been an FPE! with
0915           //  fHFEnergySimHits[cell.rawId()]==0!\n";
0916           //}
0917           mehHcalSHEvAEE[0]->Fill(fDigiSum, fHFEnergySimHits[cell.rawId()]);
0918         }
0919 
0920         // get HE info
0921         if (cell.subdet() == sdHcalEC) {
0922           ++iHE;
0923           float fDigiSum = 0.0;
0924           for (int ii = 0; ii < tool.size(); ++ii) {
0925             int capid = (*ihbhe)[ii].capid();
0926             fDigiSum += (tool[ii] - calibrations.pedestal(capid));
0927           }
0928 
0929           mehHcalSHE[1]->Fill(fHFEnergySimHits[cell.rawId()]);
0930           mehHcalAEE[1]->Fill(fDigiSum);
0931           // Adding protection against FPE
0932           if (fHFEnergySimHits[cell.rawId()] != 0) {
0933             mehHcalAEESHE[1]->Fill(fDigiSum / fHFEnergySimHits[cell.rawId()]);
0934           }
0935           // else{
0936           //  std::cout<<"It would have been an FPE! with
0937           //  fHFEnergySimHits[cell.rawId()]==0!\n";
0938           //}
0939           mehHcalSHEvAEE[1]->Fill(fDigiSum, fHFEnergySimHits[cell.rawId()]);
0940         }
0941       }
0942     }
0943 
0944     if (verbosity > 1) {
0945       eventout += "\n          Number of HBDigis collected:.............. ";
0946       eventout += iHB;
0947     }
0948     mehHcaln[0]->Fill((float)iHB);
0949 
0950     if (verbosity > 1) {
0951       eventout += "\n          Number of HEDigis collected:.............. ";
0952       eventout += iHE;
0953     }
0954     mehHcaln[1]->Fill((float)iHE);
0955   }
0956 
0957   ////////////////////////
0958   // get HO information
0959   ///////////////////////
0960   edm::Handle<edm::SortedCollection<HODataFrame>> ho;
0961   iEvent.getByToken(HODigi_Token_, ho);
0962   bool validHO = true;
0963   if (!ho.isValid()) {
0964     LogDebug(MsgLoggerCat) << "Unable to find HODataFrame in event!";
0965     validHO = false;
0966   }
0967   if (validHO) {
0968     edm::SortedCollection<HODataFrame>::const_iterator iho;
0969 
0970     int iHO = 0;
0971     for (iho = ho->begin(); iho != ho->end(); ++iho) {
0972       HcalDetId cell(iho->id());
0973 
0974       if (cell.subdet() == sdHcalOut) {
0975         // HCalconditions->makeHcalCalibration(cell, &calibrations);
0976         const HcalCalibrations &calibrations = HCalconditions->getHcalCalibrations(cell);
0977         const HcalQIECoder *channelCoder = HCalconditions->getHcalCoder(cell);
0978         const HcalQIEShape *shape = HCalconditions->getHcalShape(channelCoder);
0979         HcalCoderDb coder(*channelCoder, *shape);
0980         coder.adc2fC(*iho, tool);
0981 
0982         ++iHO;
0983         float fDigiSum = 0.0;
0984         for (int ii = 0; ii < tool.size(); ++ii) {
0985           // default ped is 4.5
0986           int capid = (*iho)[ii].capid();
0987           fDigiSum += (tool[ii] - calibrations.pedestal(capid));
0988         }
0989 
0990         mehHcalSHE[2]->Fill(fHFEnergySimHits[cell.rawId()]);
0991         mehHcalAEE[2]->Fill(fDigiSum);
0992         // Adding protection against FPE
0993         if (fHFEnergySimHits[cell.rawId()] != 0) {
0994           mehHcalAEESHE[2]->Fill(fDigiSum / fHFEnergySimHits[cell.rawId()]);
0995         }
0996         // else{
0997         //  std::cout<<"It would have been an FPE! with
0998         //  fHFEnergySimHits[cell.rawId()]==0!\n";
0999         //}
1000         mehHcalSHEvAEE[2]->Fill(fDigiSum, fHFEnergySimHits[cell.rawId()]);
1001       }
1002     }
1003 
1004     if (verbosity > 1) {
1005       eventout += "\n          Number of HODigis collected:.............. ";
1006       eventout += iHO;
1007     }
1008     mehHcaln[2]->Fill((float)iHO);
1009   }
1010 
1011   ////////////////////////
1012   // get HF information
1013   ///////////////////////
1014   edm::Handle<edm::SortedCollection<HFDataFrame>> hf;
1015   iEvent.getByToken(HFDigi_Token_, hf);
1016   bool validHF = true;
1017   if (!hf.isValid()) {
1018     LogDebug(MsgLoggerCat) << "Unable to find HFDataFrame in event!";
1019     validHF = false;
1020   }
1021   if (validHF) {
1022     edm::SortedCollection<HFDataFrame>::const_iterator ihf;
1023 
1024     int iHF = 0;
1025     for (ihf = hf->begin(); ihf != hf->end(); ++ihf) {
1026       HcalDetId cell(ihf->id());
1027 
1028       if (cell.subdet() == sdHcalFwd) {
1029         // HCalconditions->makeHcalCalibration(cell, &calibrations);
1030         const HcalCalibrations &calibrations = HCalconditions->getHcalCalibrations(cell);
1031         const HcalQIECoder *channelCoder = HCalconditions->getHcalCoder(cell);
1032         const HcalQIEShape *shape = HCalconditions->getHcalShape(channelCoder);
1033         HcalCoderDb coder(*channelCoder, *shape);
1034         coder.adc2fC(*ihf, tool);
1035 
1036         ++iHF;
1037         float fDigiSum = 0.0;
1038         for (int ii = 0; ii < tool.size(); ++ii) {
1039           // default ped is 1.73077
1040           int capid = (*ihf)[ii].capid();
1041           fDigiSum += (tool[ii] - calibrations.pedestal(capid));
1042         }
1043 
1044         mehHcalSHE[3]->Fill(fHFEnergySimHits[cell.rawId()]);
1045         mehHcalAEE[3]->Fill(fDigiSum);
1046         // Adding protection against FPE
1047         if (fHFEnergySimHits[cell.rawId()] != 0) {
1048           mehHcalAEESHE[3]->Fill(fDigiSum / fHFEnergySimHits[cell.rawId()]);
1049         }
1050         // else{
1051         //  std::cout<<"It would have been an FPE! with
1052         //  fHFEnergySimHits[cell.rawId()]==0!\n";
1053         //}
1054         mehHcalSHEvAEE[3]->Fill(fDigiSum, fHFEnergySimHits[cell.rawId()]);
1055       }
1056     }
1057 
1058     if (verbosity > 1) {
1059       eventout += "\n          Number of HFDigis collected:.............. ";
1060       eventout += iHF;
1061     }
1062     mehHcaln[3]->Fill((float)iHF);
1063   }
1064 
1065   if (verbosity > 0)
1066     edm::LogInfo(MsgLoggerCat) << eventout << "\n";
1067 
1068   return;
1069 }
1070 
1071 void GlobalDigisAnalyzer::fillTrk(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
1072   // Retrieve tracker topology from geometry
1073   const TrackerTopology *const tTopo = &iSetup.getData(tTopoToken_);
1074   std::string MsgLoggerCat = "GlobalDigisAnalyzer_fillTrk";
1075   TString eventout;
1076   if (verbosity > 0)
1077     eventout = "\nGathering info:";
1078 
1079   // get strip information
1080   edm::Handle<edm::DetSetVector<SiStripDigi>> stripDigis;
1081   iEvent.getByToken(SiStripSrc_Token_, stripDigis);
1082   bool validstripDigis = true;
1083   if (!stripDigis.isValid()) {
1084     LogDebug(MsgLoggerCat) << "Unable to find stripDigis in event!";
1085     validstripDigis = false;
1086   }
1087 
1088   if (validstripDigis) {
1089     int nStripBrl = 0, nStripFwd = 0;
1090     edm::DetSetVector<SiStripDigi>::const_iterator DSViter;
1091     for (DSViter = stripDigis->begin(); DSViter != stripDigis->end(); ++DSViter) {
1092       unsigned int id = DSViter->id;
1093       DetId detId(id);
1094       edm::DetSet<SiStripDigi>::const_iterator begin = DSViter->data.begin();
1095       edm::DetSet<SiStripDigi>::const_iterator end = DSViter->data.end();
1096       edm::DetSet<SiStripDigi>::const_iterator iter;
1097 
1098       // get TIB
1099       if (detId.subdetId() == sdSiTIB) {
1100         for (iter = begin; iter != end; ++iter) {
1101           ++nStripBrl;
1102           if (tTopo->tibLayer(id) == 1) {
1103             mehSiStripADC[0]->Fill((*iter).adc());
1104             mehSiStripStrip[0]->Fill((*iter).strip());
1105           }
1106           if (tTopo->tibLayer(id) == 2) {
1107             mehSiStripADC[1]->Fill((*iter).adc());
1108             mehSiStripStrip[1]->Fill((*iter).strip());
1109           }
1110           if (tTopo->tibLayer(id) == 3) {
1111             mehSiStripADC[2]->Fill((*iter).adc());
1112             mehSiStripStrip[2]->Fill((*iter).strip());
1113           }
1114           if (tTopo->tibLayer(id) == 4) {
1115             mehSiStripADC[3]->Fill((*iter).adc());
1116             mehSiStripStrip[3]->Fill((*iter).strip());
1117           }
1118         }
1119       }
1120 
1121       // get TOB
1122       if (detId.subdetId() == sdSiTOB) {
1123         for (iter = begin; iter != end; ++iter) {
1124           ++nStripBrl;
1125           if (tTopo->tobLayer(id) == 1) {
1126             mehSiStripADC[4]->Fill((*iter).adc());
1127             mehSiStripStrip[4]->Fill((*iter).strip());
1128           }
1129           if (tTopo->tobLayer(id) == 2) {
1130             mehSiStripADC[5]->Fill((*iter).adc());
1131             mehSiStripStrip[5]->Fill((*iter).strip());
1132           }
1133           if (tTopo->tobLayer(id) == 3) {
1134             mehSiStripADC[6]->Fill((*iter).adc());
1135             mehSiStripStrip[6]->Fill((*iter).strip());
1136           }
1137           if (tTopo->tobLayer(id) == 4) {
1138             mehSiStripADC[7]->Fill((*iter).adc());
1139             mehSiStripStrip[7]->Fill((*iter).strip());
1140           }
1141         }
1142       }
1143 
1144       // get TID
1145       if (detId.subdetId() == sdSiTID) {
1146         for (iter = begin; iter != end; ++iter) {
1147           ++nStripFwd;
1148           if (tTopo->tidWheel(id) == 1) {
1149             mehSiStripADC[8]->Fill((*iter).adc());
1150             mehSiStripStrip[8]->Fill((*iter).strip());
1151           }
1152           if (tTopo->tidWheel(id) == 2) {
1153             mehSiStripADC[9]->Fill((*iter).adc());
1154             mehSiStripStrip[9]->Fill((*iter).strip());
1155           }
1156           if (tTopo->tidWheel(id) == 3) {
1157             mehSiStripADC[10]->Fill((*iter).adc());
1158             mehSiStripStrip[10]->Fill((*iter).strip());
1159           }
1160         }
1161       }
1162 
1163       // get TEC
1164       if (detId.subdetId() == sdSiTEC) {
1165         for (iter = begin; iter != end; ++iter) {
1166           ++nStripFwd;
1167           if (tTopo->tecWheel(id) == 1) {
1168             mehSiStripADC[11]->Fill((*iter).adc());
1169             mehSiStripStrip[11]->Fill((*iter).strip());
1170           }
1171           if (tTopo->tecWheel(id) == 2) {
1172             mehSiStripADC[12]->Fill((*iter).adc());
1173             mehSiStripStrip[12]->Fill((*iter).strip());
1174           }
1175           if (tTopo->tecWheel(id) == 3) {
1176             mehSiStripADC[13]->Fill((*iter).adc());
1177             mehSiStripStrip[13]->Fill((*iter).strip());
1178           }
1179           if (tTopo->tecWheel(id) == 4) {
1180             mehSiStripADC[14]->Fill((*iter).adc());
1181             mehSiStripStrip[14]->Fill((*iter).strip());
1182           }
1183           if (tTopo->tecWheel(id) == 5) {
1184             mehSiStripADC[15]->Fill((*iter).adc());
1185             mehSiStripStrip[15]->Fill((*iter).strip());
1186           }
1187           if (tTopo->tecWheel(id) == 6) {
1188             mehSiStripADC[16]->Fill((*iter).adc());
1189             mehSiStripStrip[16]->Fill((*iter).strip());
1190           }
1191           if (tTopo->tecWheel(id) == 7) {
1192             mehSiStripADC[17]->Fill((*iter).adc());
1193             mehSiStripStrip[17]->Fill((*iter).strip());
1194           }
1195           if (tTopo->tecWheel(id) == 8) {
1196             mehSiStripADC[18]->Fill((*iter).adc());
1197             mehSiStripStrip[18]->Fill((*iter).strip());
1198           }
1199         }
1200       }
1201     }  // end loop over DataSetVector
1202 
1203     if (verbosity > 1) {
1204       eventout += "\n          Number of BrlStripDigis collected:........ ";
1205       eventout += nStripBrl;
1206     }
1207     for (int i = 0; i < 8; ++i) {
1208       mehSiStripn[i]->Fill((float)nStripBrl);
1209     }
1210 
1211     if (verbosity > 1) {
1212       eventout += "\n          Number of FrwdStripDigis collected:....... ";
1213       eventout += nStripFwd;
1214     }
1215     for (int i = 8; i < 19; ++i) {
1216       mehSiStripn[i]->Fill((float)nStripFwd);
1217     }
1218   }
1219 
1220   // get pixel information
1221   edm::Handle<edm::DetSetVector<PixelDigi>> pixelDigis;
1222   iEvent.getByToken(SiPxlSrc_Token_, pixelDigis);
1223   bool validpixelDigis = true;
1224   if (!pixelDigis.isValid()) {
1225     LogDebug(MsgLoggerCat) << "Unable to find pixelDigis in event!";
1226     validpixelDigis = false;
1227   }
1228   if (validpixelDigis) {
1229     int nPxlBrl = 0, nPxlFwd = 0;
1230     edm::DetSetVector<PixelDigi>::const_iterator DPViter;
1231     for (DPViter = pixelDigis->begin(); DPViter != pixelDigis->end(); ++DPViter) {
1232       unsigned int id = DPViter->id;
1233       DetId detId(id);
1234       edm::DetSet<PixelDigi>::const_iterator begin = DPViter->data.begin();
1235       edm::DetSet<PixelDigi>::const_iterator end = DPViter->data.end();
1236       edm::DetSet<PixelDigi>::const_iterator iter;
1237 
1238       // get Barrel pixels
1239       if (detId.subdetId() == sdPxlBrl) {
1240         for (iter = begin; iter != end; ++iter) {
1241           ++nPxlBrl;
1242           if (tTopo->pxbLayer(id) == 1) {
1243             mehSiPixelADC[0]->Fill((*iter).adc());
1244             mehSiPixelRow[0]->Fill((*iter).row());
1245             mehSiPixelCol[0]->Fill((*iter).column());
1246           }
1247           if (tTopo->pxbLayer(id) == 2) {
1248             mehSiPixelADC[1]->Fill((*iter).adc());
1249             mehSiPixelRow[1]->Fill((*iter).row());
1250             mehSiPixelCol[1]->Fill((*iter).column());
1251           }
1252           if (tTopo->pxbLayer(id) == 3) {
1253             mehSiPixelADC[2]->Fill((*iter).adc());
1254             mehSiPixelRow[2]->Fill((*iter).row());
1255             mehSiPixelCol[2]->Fill((*iter).column());
1256           }
1257         }
1258       }
1259 
1260       // get Forward pixels
1261       if (detId.subdetId() == sdPxlFwd) {
1262         for (iter = begin; iter != end; ++iter) {
1263           ++nPxlFwd;
1264           if (tTopo->pxfDisk(id) == 1) {
1265             if (tTopo->pxfSide(id) == 1) {
1266               mehSiPixelADC[4]->Fill((*iter).adc());
1267               mehSiPixelRow[4]->Fill((*iter).row());
1268               mehSiPixelCol[4]->Fill((*iter).column());
1269             }
1270             if (tTopo->pxfSide(id) == 2) {
1271               mehSiPixelADC[3]->Fill((*iter).adc());
1272               mehSiPixelRow[3]->Fill((*iter).row());
1273               mehSiPixelCol[3]->Fill((*iter).column());
1274             }
1275           }
1276           if (tTopo->pxfDisk(id) == 2) {
1277             if (tTopo->pxfSide(id) == 1) {
1278               mehSiPixelADC[6]->Fill((*iter).adc());
1279               mehSiPixelRow[6]->Fill((*iter).row());
1280               mehSiPixelCol[6]->Fill((*iter).column());
1281             }
1282             if (tTopo->pxfSide(id) == 2) {
1283               mehSiPixelADC[5]->Fill((*iter).adc());
1284               mehSiPixelRow[5]->Fill((*iter).row());
1285               mehSiPixelCol[5]->Fill((*iter).column());
1286             }
1287           }
1288         }
1289       }
1290     }
1291 
1292     if (verbosity > 1) {
1293       eventout += "\n          Number of BrlPixelDigis collected:........ ";
1294       eventout += nPxlBrl;
1295     }
1296     for (int i = 0; i < 3; ++i) {
1297       mehSiPixeln[i]->Fill((float)nPxlBrl);
1298     }
1299 
1300     if (verbosity > 1) {
1301       eventout += "\n          Number of FrwdPixelDigis collected:....... ";
1302       eventout += nPxlFwd;
1303     }
1304 
1305     for (int i = 3; i < 7; ++i) {
1306       mehSiPixeln[i]->Fill((float)nPxlFwd);
1307     }
1308   }
1309 
1310   if (verbosity > 0)
1311     edm::LogInfo(MsgLoggerCat) << eventout << "\n";
1312 
1313   return;
1314 }
1315 
1316 void GlobalDigisAnalyzer::fillMuon(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
1317   std::string MsgLoggerCat = "GlobalDigisAnalyzer_fillMuon";
1318 
1319   TString eventout;
1320   if (verbosity > 0)
1321     eventout = "\nGathering info:";
1322 
1323   // get DT information
1324   edm::Handle<DTDigiCollection> dtDigis;
1325   iEvent.getByToken(MuDTSrc_Token_, dtDigis);
1326   bool validdtDigis = true;
1327   if (!dtDigis.isValid()) {
1328     LogDebug(MsgLoggerCat) << "Unable to find dtDigis in event!";
1329     validdtDigis = false;
1330   }
1331   if (validdtDigis) {
1332     int nDt0 = 0;
1333     int nDt1 = 0;
1334     int nDt2 = 0;
1335     int nDt3 = 0;
1336     int nDt = 0;
1337     DTDigiCollection::DigiRangeIterator detUnitIt;
1338     for (detUnitIt = dtDigis->begin(); detUnitIt != dtDigis->end(); ++detUnitIt) {
1339       const DTLayerId &id = (*detUnitIt).first;
1340       const DTDigiCollection::Range &range = (*detUnitIt).second;
1341       for (DTDigiCollection::const_iterator digiIt = range.first; digiIt != range.second; ++digiIt) {
1342         ++nDt;
1343 
1344         DTWireId wireId(id, (*digiIt).wire());
1345         if (wireId.station() == 1) {
1346           mehDtMuonLayer[0]->Fill(id.layer());
1347           mehDtMuonTime[0]->Fill((*digiIt).time());
1348           mehDtMuonTimevLayer[0]->Fill(id.layer(), (*digiIt).time());
1349           ++nDt0;
1350         }
1351         if (wireId.station() == 2) {
1352           mehDtMuonLayer[1]->Fill(id.layer());
1353           mehDtMuonTime[1]->Fill((*digiIt).time());
1354           mehDtMuonTimevLayer[1]->Fill(id.layer(), (*digiIt).time());
1355           ++nDt1;
1356         }
1357         if (wireId.station() == 3) {
1358           mehDtMuonLayer[2]->Fill(id.layer());
1359           mehDtMuonTime[2]->Fill((*digiIt).time());
1360           mehDtMuonTimevLayer[2]->Fill(id.layer(), (*digiIt).time());
1361           ++nDt2;
1362         }
1363         if (wireId.station() == 4) {
1364           mehDtMuonLayer[3]->Fill(id.layer());
1365           mehDtMuonTime[3]->Fill((*digiIt).time());
1366           mehDtMuonTimevLayer[3]->Fill(id.layer(), (*digiIt).time());
1367           ++nDt3;
1368         }
1369       }
1370     }
1371     mehDtMuonn[0]->Fill((float)nDt0);
1372     mehDtMuonn[1]->Fill((float)nDt1);
1373     mehDtMuonn[2]->Fill((float)nDt2);
1374     mehDtMuonn[3]->Fill((float)nDt3);
1375 
1376     if (verbosity > 1) {
1377       eventout += "\n          Number of DtMuonDigis collected:.......... ";
1378       eventout += nDt;
1379     }
1380   }
1381 
1382   // get CSC Strip information
1383   edm::Handle<CSCStripDigiCollection> strips;
1384   iEvent.getByToken(MuCSCStripSrc_Token_, strips);
1385   bool validstrips = true;
1386   if (!strips.isValid()) {
1387     LogDebug(MsgLoggerCat) << "Unable to find muon strips in event!";
1388     validstrips = false;
1389   }
1390 
1391   if (validstrips) {
1392     int nStrips = 0;
1393     for (CSCStripDigiCollection::DigiRangeIterator j = strips->begin(); j != strips->end(); ++j) {
1394       std::vector<CSCStripDigi>::const_iterator digiItr = (*j).second.first;
1395       std::vector<CSCStripDigi>::const_iterator last = (*j).second.second;
1396 
1397       for (; digiItr != last; ++digiItr) {
1398         ++nStrips;
1399 
1400         // average pedestals
1401         std::vector<int> adcCounts = digiItr->getADCCounts();
1402         theCSCStripPedestalSum += adcCounts[0];
1403         theCSCStripPedestalSum += adcCounts[1];
1404         theCSCStripPedestalCount += 2;
1405 
1406         // if there are enough pedestal statistics
1407         if (theCSCStripPedestalCount > 100) {
1408           float pedestal = theCSCStripPedestalSum / theCSCStripPedestalCount;
1409           if (adcCounts[5] > (pedestal + 100))
1410             mehCSCStripADC->Fill(adcCounts[4] - pedestal);
1411         }
1412       }
1413     }
1414 
1415     if (verbosity > 1) {
1416       eventout += "\n          Number of CSCStripDigis collected:........ ";
1417       eventout += nStrips;
1418     }
1419     mehCSCStripn->Fill((float)nStrips);
1420   }
1421 
1422   // get CSC Wire information
1423   edm::Handle<CSCWireDigiCollection> wires;
1424   iEvent.getByToken(MuCSCWireSrc_Token_, wires);
1425   bool validwires = true;
1426   if (!wires.isValid()) {
1427     LogDebug(MsgLoggerCat) << "Unable to find muon wires in event!";
1428     validwires = false;
1429   }
1430 
1431   if (validwires) {
1432     int nWires = 0;
1433     for (CSCWireDigiCollection::DigiRangeIterator j = wires->begin(); j != wires->end(); ++j) {
1434       std::vector<CSCWireDigi>::const_iterator digiItr = (*j).second.first;
1435       std::vector<CSCWireDigi>::const_iterator endDigi = (*j).second.second;
1436 
1437       for (; digiItr != endDigi; ++digiItr) {
1438         ++nWires;
1439         mehCSCWireTime->Fill(digiItr->getTimeBin());
1440       }
1441     }
1442 
1443     if (verbosity > 1) {
1444       eventout += "\n          Number of CSCWireDigis collected:......... ";
1445       eventout += nWires;
1446     }
1447     mehCSCWiren->Fill((float)nWires);
1448   }
1449 
1450   // get RPC information
1451   // Get the RPC Geometry
1452   const auto &rpcGeom = iSetup.getHandle(rpcGeomToken_);
1453   if (!rpcGeom.isValid()) {
1454     edm::LogWarning(MsgLoggerCat) << "Unable to find RPCGeometryRecord in event!";
1455     return;
1456   }
1457 
1458   edm::Handle<edm::PSimHitContainer> rpcsimHit;
1459   iEvent.getByToken(RPCSimHit_Token_, rpcsimHit);
1460   bool validrpcsim = true;
1461   if (!rpcsimHit.isValid()) {
1462     LogDebug(MsgLoggerCat) << "Unable to find rpcsimHit in event!";
1463     validrpcsim = false;
1464   }
1465 
1466   edm::Handle<RPCDigiCollection> rpcDigis;
1467   iEvent.getByToken(MuRPCSrc_Token_, rpcDigis);
1468   bool validrpcdigi = true;
1469   if (!rpcDigis.isValid()) {
1470     LogDebug(MsgLoggerCat) << "Unable to find rpcDigis in event!";
1471     [[clang::suppress]] validrpcdigi = false;
1472   }
1473 
1474   // ONLY UNTIL PROBLEM WITH RPC DIGIS IS FIGURED OUT
1475   validrpcdigi = false;
1476 
1477   // Loop on simhits
1478   edm::PSimHitContainer::const_iterator rpcsimIt;
1479   std::map<RPCDetId, std::vector<double>> allsims;
1480 
1481   if (validrpcsim) {
1482     for (rpcsimIt = rpcsimHit->begin(); rpcsimIt != rpcsimHit->end(); rpcsimIt++) {
1483       RPCDetId Rsid = (RPCDetId)(*rpcsimIt).detUnitId();
1484       int ptype = rpcsimIt->particleType();
1485 
1486       if (ptype == 13 || ptype == -13) {
1487         std::vector<double> buff;
1488         if (allsims.find(Rsid) != allsims.end()) {
1489           buff = allsims[Rsid];
1490         }
1491         buff.push_back(rpcsimIt->localPosition().x());
1492         allsims[Rsid] = buff;
1493       }
1494     }
1495   }
1496 
1497   // CRASH HAPPENS SOMEWHERE HERE IN FOR DECLARATION
1498   // WHAT IS WRONG WITH rpcDigis?????
1499   if (validrpcdigi) {
1500     int nRPC = 0;
1501     RPCDigiCollection::DigiRangeIterator rpcdetUnitIt;
1502     for (rpcdetUnitIt = rpcDigis->begin(); rpcdetUnitIt != rpcDigis->end(); ++rpcdetUnitIt) {
1503       const RPCDetId Rsid = (*rpcdetUnitIt).first;
1504       const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(Rsid));
1505       const RPCDigiCollection::Range &range = (*rpcdetUnitIt).second;
1506 
1507       std::vector<double> sims;
1508       if (allsims.find(Rsid) != allsims.end()) {
1509         sims = allsims[Rsid];
1510       }
1511 
1512       int ndigi = 0;
1513       for (RPCDigiCollection::const_iterator rpcdigiIt = range.first; rpcdigiIt != range.second; ++rpcdigiIt) {
1514         ++ndigi;
1515         ++nRPC;
1516       }
1517 
1518       if (sims.size() == 1 && ndigi == 1) {
1519         double dis = roll->centreOfStrip(range.first->strip()).x() - sims[0];
1520 
1521         if (Rsid.region() == 0) {
1522           if (Rsid.ring() == -2)
1523             mehRPCRes[0]->Fill((float)dis);
1524           else if (Rsid.ring() == -1)
1525             mehRPCRes[1]->Fill((float)dis);
1526           else if (Rsid.ring() == 0)
1527             mehRPCRes[2]->Fill((float)dis);
1528           else if (Rsid.ring() == 1)
1529             mehRPCRes[3]->Fill((float)dis);
1530           else if (Rsid.ring() == 2)
1531             mehRPCRes[4]->Fill((float)dis);
1532         }
1533       }
1534     }
1535 
1536     if (verbosity > 1) {
1537       eventout += "\n          Number of RPCDigis collected:.............. ";
1538       eventout += nRPC;
1539     }
1540     mehRPCMuonn->Fill(float(nRPC));
1541   }
1542 
1543   if (verbosity > 0)
1544     edm::LogInfo(MsgLoggerCat) << eventout << "\n";
1545 
1546   return;
1547 }