Back to home page

Project CMSSW displayed by LXR

 
 

    


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

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 pedestalPreSample = 0.;
0574         double pedestalPreSampleAnalog = 0.;
0575 
0576         for (int sample = 0; sample < nrSamples; ++sample) {
0577           ebAnalogSignal[sample] = 0.;
0578           ebADCCounts[sample] = 0.;
0579           ebADCGains[sample] = -1.;
0580         }
0581 
0582         // calculate maximum energy and pedestal
0583         for (int sample = 0; sample < nrSamples; ++sample) {
0584           EcalMGPASample thisSample = ebdf[sample];
0585           ebADCCounts[sample] = (thisSample.adc());
0586           ebADCGains[sample] = (thisSample.gainId());
0587           ebAnalogSignal[sample] = (ebADCCounts[sample] * ECalgainConv_[(int)ebADCGains[sample]] * ECalbarrelADCtoGeV_);
0588           if (Emax < ebAnalogSignal[sample]) {
0589             Emax = ebAnalogSignal[sample];
0590             Pmax = sample;
0591           }
0592           if (sample < 3) {
0593             pedestalPreSample += ebADCCounts[sample];
0594             pedestalPreSampleAnalog +=
0595                 ebADCCounts[sample] * ECalgainConv_[(int)ebADCGains[sample]] * ECalbarrelADCtoGeV_;
0596           }
0597         }
0598         pedestalPreSample /= 3.;
0599         pedestalPreSampleAnalog /= 3.;
0600 
0601         // calculate pedestal subtracted digi energy in the crystal
0602         double Erec = Emax - pedestalPreSampleAnalog * ECalgainConv_[(int)ebADCGains[Pmax]];
0603 
0604         // gather necessary information
0605         mehEcalMaxPos[0]->Fill(Pmax);
0606         mehEcalSHE[0]->Fill(ebSimMap[ebid.rawId()]);
0607         mehEcalAEE[0]->Fill(Erec);
0608         // Adding protection against FPE
0609         if (ebSimMap[ebid.rawId()] != 0) {
0610           mehEcalSHEvAEESHE[0]->Fill(Erec / ebSimMap[ebid.rawId()], ebSimMap[ebid.rawId()]);
0611         }
0612         // else {
0613         //  std::cout<<"Would have been an FPE! with
0614         //  ebSimMap[ebid.rawId()]==0\n";
0615         //}
0616 
0617         mehEcalMultvAEE[0]->Fill(Pmax, (float)i);
0618       }
0619 
0620       if (verbosity > 1) {
0621         eventout += "\n          Number of EBDigis collected:.............. ";
0622         eventout += i;
0623       }
0624       mehEcaln[0]->Fill((float)i);
0625     }
0626   }
0627 
0628   /////////////////////////
0629   // extract EE information
0630   ////////////////////////
0631   bool isEndCap = true;
0632   edm::Handle<EEDigiCollection> EcalDigiEE;
0633   iEvent.getByToken(ECalEESrc_Token_, EcalDigiEE);
0634   bool validDigiEE = true;
0635   if (!EcalDigiEE.isValid()) {
0636     LogDebug(MsgLoggerCat) << "Unable to find EcalDigiEE in event!";
0637     validDigiEE = false;
0638   }
0639   if (validDigiEE) {
0640     if (EcalDigiEE->empty())
0641       isEndCap = false;
0642 
0643     if (isEndCap) {
0644       // loop over simhits
0645       MapType eeSimMap;
0646       iEvent.getByToken(EEHits_Token_, crossingFrame);
0647       bool validXFrame = true;
0648       if (!crossingFrame.isValid()) {
0649         LogDebug(MsgLoggerCat) << "Unable to find cal endcap crossingFrame in event!";
0650         validXFrame = false;
0651       }
0652       if (validXFrame) {
0653         const MixCollection<PCaloHit> endcapHits(crossingFrame.product());
0654 
0655         // keep track of sum of simhit energy in each crystal
0656         for (auto const &iHit : endcapHits) {
0657           EEDetId eeid = EEDetId(iHit.id());
0658 
0659           uint32_t crystid = eeid.rawId();
0660           eeSimMap[crystid] += iHit.energy();
0661         }
0662       }
0663 
0664       // loop over digis
0665       const EEDigiCollection *endcapDigi = EcalDigiEE.product();
0666 
0667       std::vector<double> eeAnalogSignal;
0668       std::vector<double> eeADCCounts;
0669       std::vector<double> eeADCGains;
0670       eeAnalogSignal.reserve(EEDataFrame::MAXSAMPLES);
0671       eeADCCounts.reserve(EEDataFrame::MAXSAMPLES);
0672       eeADCGains.reserve(EEDataFrame::MAXSAMPLES);
0673 
0674       int inc = 0;
0675       for (unsigned int digis = 0; digis < EcalDigiEE->size(); ++digis) {
0676         ++inc;
0677 
0678         EEDataFrame eedf = (*endcapDigi)[digis];
0679         int nrSamples = eedf.size();
0680 
0681         EEDetId eeid = eedf.id();
0682 
0683         double Emax = 0;
0684         int Pmax = 0;
0685         double pedestalPreSample = 0.;
0686         double pedestalPreSampleAnalog = 0.;
0687 
0688         for (int sample = 0; sample < nrSamples; ++sample) {
0689           eeAnalogSignal[sample] = 0.;
0690           eeADCCounts[sample] = 0.;
0691           eeADCGains[sample] = -1.;
0692         }
0693 
0694         // calculate maximum enery and pedestal
0695         for (int sample = 0; sample < nrSamples; ++sample) {
0696           EcalMGPASample thisSample = eedf[sample];
0697 
0698           eeADCCounts[sample] = (thisSample.adc());
0699           eeADCGains[sample] = (thisSample.gainId());
0700           eeAnalogSignal[sample] = (eeADCCounts[sample] * ECalgainConv_[(int)eeADCGains[sample]] * ECalbarrelADCtoGeV_);
0701           if (Emax < eeAnalogSignal[sample]) {
0702             Emax = eeAnalogSignal[sample];
0703             Pmax = sample;
0704           }
0705           if (sample < 3) {
0706             pedestalPreSample += eeADCCounts[sample];
0707             pedestalPreSampleAnalog +=
0708                 eeADCCounts[sample] * ECalgainConv_[(int)eeADCGains[sample]] * ECalbarrelADCtoGeV_;
0709           }
0710         }
0711         pedestalPreSample /= 3.;
0712         pedestalPreSampleAnalog /= 3.;
0713 
0714         // calculate pedestal subtracted digi energy in the crystal
0715         double Erec = Emax - pedestalPreSampleAnalog * ECalgainConv_[(int)eeADCGains[Pmax]];
0716 
0717         // gather necessary information
0718         mehEcalMaxPos[1]->Fill(Pmax);
0719         mehEcalSHE[1]->Fill(eeSimMap[eeid.rawId()]);
0720         mehEcalAEE[1]->Fill(Erec);
0721         // Adding protection against FPE
0722         if (eeSimMap[eeid.rawId()] != 0) {
0723           mehEcalSHEvAEESHE[1]->Fill(Erec / eeSimMap[eeid.rawId()], eeSimMap[eeid.rawId()]);
0724         }
0725         // else{
0726         //  std::cout<<"Would have been an FPE! with
0727         //  eeSimMap[eeid.rawId()]==0\n";
0728         //}
0729         mehEcalMultvAEE[1]->Fill(Pmax, (float)inc);
0730       }
0731 
0732       if (verbosity > 1) {
0733         eventout += "\n          Number of EEDigis collected:.............. ";
0734         eventout += inc;
0735       }
0736 
0737       mehEcaln[1]->Fill((float)inc);
0738     }
0739   }
0740 
0741   /////////////////////////
0742   // extract ES information
0743   ////////////////////////
0744   bool isPreshower = true;
0745   edm::Handle<ESDigiCollection> EcalDigiES;
0746   iEvent.getByToken(ECalESSrc_Token_, EcalDigiES);
0747   bool validDigiES = true;
0748   if (!EcalDigiES.isValid()) {
0749     LogDebug(MsgLoggerCat) << "Unable to find EcalDigiES in event!";
0750     validDigiES = false;
0751   }
0752 
0753   // ONLY WHILE GEOMETRY IS REMOVED
0754   validDigiES = false;
0755 
0756   if (validDigiES) {
0757     if (EcalDigiES->empty())
0758       isPreshower = false;
0759 
0760     if (isPreshower) {
0761       // loop over simhits
0762       iEvent.getByToken(ESHits_Token_, crossingFrame);
0763       bool validXFrame = true;
0764       if (!crossingFrame.isValid()) {
0765         LogDebug(MsgLoggerCat) << "Unable to find cal preshower crossingFrame in event!";
0766         validXFrame = false;
0767       }
0768       if (validXFrame) {
0769         const MixCollection<PCaloHit> preshowerHits(crossingFrame.product());
0770 
0771         // keep track of sum of simhit energy in each crystal
0772         MapType esSimMap;
0773         for (auto const &iHit : preshowerHits) {
0774           ESDetId esid = ESDetId(iHit.id());
0775 
0776           uint32_t crystid = esid.rawId();
0777           esSimMap[crystid] += iHit.energy();
0778         }
0779       }
0780 
0781       // loop over digis
0782       const ESDigiCollection *preshowerDigi = EcalDigiES.product();
0783 
0784       std::vector<double> esADCCounts;
0785       esADCCounts.reserve(ESDataFrame::MAXSAMPLES);
0786 
0787       int i = 0;
0788       for (unsigned int digis = 0; digis < EcalDigiES->size(); ++digis) {
0789         ++i;
0790 
0791         ESDataFrame esdf = (*preshowerDigi)[digis];
0792         int nrSamples = esdf.size();
0793 
0794         for (int sample = 0; sample < nrSamples; ++sample) {
0795           esADCCounts[sample] = 0.;
0796         }
0797 
0798         // gether ADC counts
0799         for (int sample = 0; sample < nrSamples; ++sample) {
0800           ESSample thisSample = esdf[sample];
0801           esADCCounts[sample] = (thisSample.adc());
0802         }
0803 
0804         mehEScalADC[0]->Fill(esADCCounts[0]);
0805         mehEScalADC[1]->Fill(esADCCounts[1]);
0806         mehEScalADC[2]->Fill(esADCCounts[2]);
0807       }
0808 
0809       if (verbosity > 1) {
0810         eventout += "\n          Number of ESDigis collected:.............. ";
0811         eventout += i;
0812       }
0813 
0814       mehEScaln->Fill((float)i);
0815     }
0816   }
0817   if (verbosity > 0)
0818     edm::LogInfo(MsgLoggerCat) << eventout << "\n";
0819 
0820   return;
0821 }
0822 
0823 void GlobalDigisAnalyzer::fillHCal(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0824   std::string MsgLoggerCat = "GlobalDigisAnalyzer_fillHCal";
0825 
0826   TString eventout;
0827   if (verbosity > 0)
0828     eventout = "\nGathering info:";
0829 
0830   // get calibration info
0831   const auto &HCalconditions = iSetup.getHandle(hcaldbToken_);
0832   if (!HCalconditions.isValid()) {
0833     edm::LogWarning(MsgLoggerCat) << "Unable to find HCalconditions in event!";
0834     return;
0835   }
0836   // HcalCalibrations calibrations;
0837   CaloSamples tool;
0838 
0839   ///////////////////////
0840   // extract simhit info
0841   //////////////////////
0842   edm::Handle<edm::PCaloHitContainer> hcalHits;
0843   iEvent.getByToken(HCalSrc_Token_, hcalHits);
0844   bool validhcalHits = true;
0845   if (!hcalHits.isValid()) {
0846     LogDebug(MsgLoggerCat) << "Unable to find hcalHits in event!";
0847     validhcalHits = false;
0848   }
0849   MapType fHBEnergySimHits;
0850   MapType fHEEnergySimHits;
0851   MapType fHOEnergySimHits;
0852   MapType fHFEnergySimHits;
0853   if (validhcalHits) {
0854     const edm::PCaloHitContainer *simhitResult = hcalHits.product();
0855 
0856     for (std::vector<PCaloHit>::const_iterator simhits = simhitResult->begin(); simhits != simhitResult->end();
0857          ++simhits) {
0858       HcalDetId detId(simhits->id());
0859       uint32_t cellid = detId.rawId();
0860 
0861       if (detId.subdet() == sdHcalBrl) {
0862         fHBEnergySimHits[cellid] += simhits->energy();
0863       }
0864       if (detId.subdet() == sdHcalEC) {
0865         fHEEnergySimHits[cellid] += simhits->energy();
0866       }
0867       if (detId.subdet() == sdHcalOut) {
0868         fHOEnergySimHits[cellid] += simhits->energy();
0869       }
0870       if (detId.subdet() == sdHcalFwd) {
0871         fHFEnergySimHits[cellid] += simhits->energy();
0872       }
0873     }
0874   }
0875 
0876   ////////////////////////
0877   // get HBHE information
0878   ///////////////////////
0879   edm::Handle<edm::SortedCollection<HBHEDataFrame>> hbhe;
0880   iEvent.getByToken(HBHEDigi_Token_, hbhe);
0881   bool validHBHE = true;
0882   if (!hbhe.isValid()) {
0883     LogDebug(MsgLoggerCat) << "Unable to find HBHEDataFrame in event!";
0884     validHBHE = false;
0885   }
0886 
0887   if (validHBHE) {
0888     edm::SortedCollection<HBHEDataFrame>::const_iterator ihbhe;
0889 
0890     int iHB = 0;
0891     int iHE = 0;
0892     for (ihbhe = hbhe->begin(); ihbhe != hbhe->end(); ++ihbhe) {
0893       HcalDetId cell(ihbhe->id());
0894 
0895       if ((cell.subdet() == sdHcalBrl) || (cell.subdet() == sdHcalEC)) {
0896         // HCalconditions->makeHcalCalibration(cell, &calibrations);
0897         const HcalCalibrations &calibrations = HCalconditions->getHcalCalibrations(cell);
0898         const HcalQIECoder *channelCoder = HCalconditions->getHcalCoder(cell);
0899         const HcalQIEShape *shape = HCalconditions->getHcalShape(channelCoder);
0900         HcalCoderDb coder(*channelCoder, *shape);
0901         coder.adc2fC(*ihbhe, tool);
0902 
0903         // get HB info
0904         if (cell.subdet() == sdHcalBrl) {
0905           ++iHB;
0906           float fDigiSum = 0.0;
0907           for (int ii = 0; ii < tool.size(); ++ii) {
0908             // default ped is 4.5
0909             int capid = (*ihbhe)[ii].capid();
0910             fDigiSum += (tool[ii] - calibrations.pedestal(capid));
0911           }
0912 
0913           mehHcalSHE[0]->Fill(fHFEnergySimHits[cell.rawId()]);
0914           mehHcalAEE[0]->Fill(fDigiSum);
0915           // Adding protection against FPE
0916           if (fHFEnergySimHits[cell.rawId()] != 0) {
0917             mehHcalAEESHE[0]->Fill(fDigiSum / fHFEnergySimHits[cell.rawId()]);
0918           }
0919           // else {
0920           //  std::cout<<"It would have been an FPE! with
0921           //  fHFEnergySimHits[cell.rawId()]==0!\n";
0922           //}
0923           mehHcalSHEvAEE[0]->Fill(fDigiSum, fHFEnergySimHits[cell.rawId()]);
0924         }
0925 
0926         // get HE info
0927         if (cell.subdet() == sdHcalEC) {
0928           ++iHE;
0929           float fDigiSum = 0.0;
0930           for (int ii = 0; ii < tool.size(); ++ii) {
0931             int capid = (*ihbhe)[ii].capid();
0932             fDigiSum += (tool[ii] - calibrations.pedestal(capid));
0933           }
0934 
0935           mehHcalSHE[1]->Fill(fHFEnergySimHits[cell.rawId()]);
0936           mehHcalAEE[1]->Fill(fDigiSum);
0937           // Adding protection against FPE
0938           if (fHFEnergySimHits[cell.rawId()] != 0) {
0939             mehHcalAEESHE[1]->Fill(fDigiSum / fHFEnergySimHits[cell.rawId()]);
0940           }
0941           // else{
0942           //  std::cout<<"It would have been an FPE! with
0943           //  fHFEnergySimHits[cell.rawId()]==0!\n";
0944           //}
0945           mehHcalSHEvAEE[1]->Fill(fDigiSum, fHFEnergySimHits[cell.rawId()]);
0946         }
0947       }
0948     }
0949 
0950     if (verbosity > 1) {
0951       eventout += "\n          Number of HBDigis collected:.............. ";
0952       eventout += iHB;
0953     }
0954     mehHcaln[0]->Fill((float)iHB);
0955 
0956     if (verbosity > 1) {
0957       eventout += "\n          Number of HEDigis collected:.............. ";
0958       eventout += iHE;
0959     }
0960     mehHcaln[1]->Fill((float)iHE);
0961   }
0962 
0963   ////////////////////////
0964   // get HO information
0965   ///////////////////////
0966   edm::Handle<edm::SortedCollection<HODataFrame>> ho;
0967   iEvent.getByToken(HODigi_Token_, ho);
0968   bool validHO = true;
0969   if (!ho.isValid()) {
0970     LogDebug(MsgLoggerCat) << "Unable to find HODataFrame in event!";
0971     validHO = false;
0972   }
0973   if (validHO) {
0974     edm::SortedCollection<HODataFrame>::const_iterator iho;
0975 
0976     int iHO = 0;
0977     for (iho = ho->begin(); iho != ho->end(); ++iho) {
0978       HcalDetId cell(iho->id());
0979 
0980       if (cell.subdet() == sdHcalOut) {
0981         // HCalconditions->makeHcalCalibration(cell, &calibrations);
0982         const HcalCalibrations &calibrations = HCalconditions->getHcalCalibrations(cell);
0983         const HcalQIECoder *channelCoder = HCalconditions->getHcalCoder(cell);
0984         const HcalQIEShape *shape = HCalconditions->getHcalShape(channelCoder);
0985         HcalCoderDb coder(*channelCoder, *shape);
0986         coder.adc2fC(*iho, tool);
0987 
0988         ++iHO;
0989         float fDigiSum = 0.0;
0990         for (int ii = 0; ii < tool.size(); ++ii) {
0991           // default ped is 4.5
0992           int capid = (*iho)[ii].capid();
0993           fDigiSum += (tool[ii] - calibrations.pedestal(capid));
0994         }
0995 
0996         mehHcalSHE[2]->Fill(fHFEnergySimHits[cell.rawId()]);
0997         mehHcalAEE[2]->Fill(fDigiSum);
0998         // Adding protection against FPE
0999         if (fHFEnergySimHits[cell.rawId()] != 0) {
1000           mehHcalAEESHE[2]->Fill(fDigiSum / fHFEnergySimHits[cell.rawId()]);
1001         }
1002         // else{
1003         //  std::cout<<"It would have been an FPE! with
1004         //  fHFEnergySimHits[cell.rawId()]==0!\n";
1005         //}
1006         mehHcalSHEvAEE[2]->Fill(fDigiSum, fHFEnergySimHits[cell.rawId()]);
1007       }
1008     }
1009 
1010     if (verbosity > 1) {
1011       eventout += "\n          Number of HODigis collected:.............. ";
1012       eventout += iHO;
1013     }
1014     mehHcaln[2]->Fill((float)iHO);
1015   }
1016 
1017   ////////////////////////
1018   // get HF information
1019   ///////////////////////
1020   edm::Handle<edm::SortedCollection<HFDataFrame>> hf;
1021   iEvent.getByToken(HFDigi_Token_, hf);
1022   bool validHF = true;
1023   if (!hf.isValid()) {
1024     LogDebug(MsgLoggerCat) << "Unable to find HFDataFrame in event!";
1025     validHF = false;
1026   }
1027   if (validHF) {
1028     edm::SortedCollection<HFDataFrame>::const_iterator ihf;
1029 
1030     int iHF = 0;
1031     for (ihf = hf->begin(); ihf != hf->end(); ++ihf) {
1032       HcalDetId cell(ihf->id());
1033 
1034       if (cell.subdet() == sdHcalFwd) {
1035         // HCalconditions->makeHcalCalibration(cell, &calibrations);
1036         const HcalCalibrations &calibrations = HCalconditions->getHcalCalibrations(cell);
1037         const HcalQIECoder *channelCoder = HCalconditions->getHcalCoder(cell);
1038         const HcalQIEShape *shape = HCalconditions->getHcalShape(channelCoder);
1039         HcalCoderDb coder(*channelCoder, *shape);
1040         coder.adc2fC(*ihf, tool);
1041 
1042         ++iHF;
1043         float fDigiSum = 0.0;
1044         for (int ii = 0; ii < tool.size(); ++ii) {
1045           // default ped is 1.73077
1046           int capid = (*ihf)[ii].capid();
1047           fDigiSum += (tool[ii] - calibrations.pedestal(capid));
1048         }
1049 
1050         mehHcalSHE[3]->Fill(fHFEnergySimHits[cell.rawId()]);
1051         mehHcalAEE[3]->Fill(fDigiSum);
1052         // Adding protection against FPE
1053         if (fHFEnergySimHits[cell.rawId()] != 0) {
1054           mehHcalAEESHE[3]->Fill(fDigiSum / fHFEnergySimHits[cell.rawId()]);
1055         }
1056         // else{
1057         //  std::cout<<"It would have been an FPE! with
1058         //  fHFEnergySimHits[cell.rawId()]==0!\n";
1059         //}
1060         mehHcalSHEvAEE[3]->Fill(fDigiSum, fHFEnergySimHits[cell.rawId()]);
1061       }
1062     }
1063 
1064     if (verbosity > 1) {
1065       eventout += "\n          Number of HFDigis collected:.............. ";
1066       eventout += iHF;
1067     }
1068     mehHcaln[3]->Fill((float)iHF);
1069   }
1070 
1071   if (verbosity > 0)
1072     edm::LogInfo(MsgLoggerCat) << eventout << "\n";
1073 
1074   return;
1075 }
1076 
1077 void GlobalDigisAnalyzer::fillTrk(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
1078   // Retrieve tracker topology from geometry
1079   const TrackerTopology *const tTopo = &iSetup.getData(tTopoToken_);
1080   std::string MsgLoggerCat = "GlobalDigisAnalyzer_fillTrk";
1081   TString eventout;
1082   if (verbosity > 0)
1083     eventout = "\nGathering info:";
1084 
1085   // get strip information
1086   edm::Handle<edm::DetSetVector<SiStripDigi>> stripDigis;
1087   iEvent.getByToken(SiStripSrc_Token_, stripDigis);
1088   bool validstripDigis = true;
1089   if (!stripDigis.isValid()) {
1090     LogDebug(MsgLoggerCat) << "Unable to find stripDigis in event!";
1091     validstripDigis = false;
1092   }
1093 
1094   if (validstripDigis) {
1095     int nStripBrl = 0, nStripFwd = 0;
1096     edm::DetSetVector<SiStripDigi>::const_iterator DSViter;
1097     for (DSViter = stripDigis->begin(); DSViter != stripDigis->end(); ++DSViter) {
1098       unsigned int id = DSViter->id;
1099       DetId detId(id);
1100       edm::DetSet<SiStripDigi>::const_iterator begin = DSViter->data.begin();
1101       edm::DetSet<SiStripDigi>::const_iterator end = DSViter->data.end();
1102       edm::DetSet<SiStripDigi>::const_iterator iter;
1103 
1104       // get TIB
1105       if (detId.subdetId() == sdSiTIB) {
1106         for (iter = begin; iter != end; ++iter) {
1107           ++nStripBrl;
1108           if (tTopo->tibLayer(id) == 1) {
1109             mehSiStripADC[0]->Fill((*iter).adc());
1110             mehSiStripStrip[0]->Fill((*iter).strip());
1111           }
1112           if (tTopo->tibLayer(id) == 2) {
1113             mehSiStripADC[1]->Fill((*iter).adc());
1114             mehSiStripStrip[1]->Fill((*iter).strip());
1115           }
1116           if (tTopo->tibLayer(id) == 3) {
1117             mehSiStripADC[2]->Fill((*iter).adc());
1118             mehSiStripStrip[2]->Fill((*iter).strip());
1119           }
1120           if (tTopo->tibLayer(id) == 4) {
1121             mehSiStripADC[3]->Fill((*iter).adc());
1122             mehSiStripStrip[3]->Fill((*iter).strip());
1123           }
1124         }
1125       }
1126 
1127       // get TOB
1128       if (detId.subdetId() == sdSiTOB) {
1129         for (iter = begin; iter != end; ++iter) {
1130           ++nStripBrl;
1131           if (tTopo->tobLayer(id) == 1) {
1132             mehSiStripADC[4]->Fill((*iter).adc());
1133             mehSiStripStrip[4]->Fill((*iter).strip());
1134           }
1135           if (tTopo->tobLayer(id) == 2) {
1136             mehSiStripADC[5]->Fill((*iter).adc());
1137             mehSiStripStrip[5]->Fill((*iter).strip());
1138           }
1139           if (tTopo->tobLayer(id) == 3) {
1140             mehSiStripADC[6]->Fill((*iter).adc());
1141             mehSiStripStrip[6]->Fill((*iter).strip());
1142           }
1143           if (tTopo->tobLayer(id) == 4) {
1144             mehSiStripADC[7]->Fill((*iter).adc());
1145             mehSiStripStrip[7]->Fill((*iter).strip());
1146           }
1147         }
1148       }
1149 
1150       // get TID
1151       if (detId.subdetId() == sdSiTID) {
1152         for (iter = begin; iter != end; ++iter) {
1153           ++nStripFwd;
1154           if (tTopo->tidWheel(id) == 1) {
1155             mehSiStripADC[8]->Fill((*iter).adc());
1156             mehSiStripStrip[8]->Fill((*iter).strip());
1157           }
1158           if (tTopo->tidWheel(id) == 2) {
1159             mehSiStripADC[9]->Fill((*iter).adc());
1160             mehSiStripStrip[9]->Fill((*iter).strip());
1161           }
1162           if (tTopo->tidWheel(id) == 3) {
1163             mehSiStripADC[10]->Fill((*iter).adc());
1164             mehSiStripStrip[10]->Fill((*iter).strip());
1165           }
1166         }
1167       }
1168 
1169       // get TEC
1170       if (detId.subdetId() == sdSiTEC) {
1171         for (iter = begin; iter != end; ++iter) {
1172           ++nStripFwd;
1173           if (tTopo->tecWheel(id) == 1) {
1174             mehSiStripADC[11]->Fill((*iter).adc());
1175             mehSiStripStrip[11]->Fill((*iter).strip());
1176           }
1177           if (tTopo->tecWheel(id) == 2) {
1178             mehSiStripADC[12]->Fill((*iter).adc());
1179             mehSiStripStrip[12]->Fill((*iter).strip());
1180           }
1181           if (tTopo->tecWheel(id) == 3) {
1182             mehSiStripADC[13]->Fill((*iter).adc());
1183             mehSiStripStrip[13]->Fill((*iter).strip());
1184           }
1185           if (tTopo->tecWheel(id) == 4) {
1186             mehSiStripADC[14]->Fill((*iter).adc());
1187             mehSiStripStrip[14]->Fill((*iter).strip());
1188           }
1189           if (tTopo->tecWheel(id) == 5) {
1190             mehSiStripADC[15]->Fill((*iter).adc());
1191             mehSiStripStrip[15]->Fill((*iter).strip());
1192           }
1193           if (tTopo->tecWheel(id) == 6) {
1194             mehSiStripADC[16]->Fill((*iter).adc());
1195             mehSiStripStrip[16]->Fill((*iter).strip());
1196           }
1197           if (tTopo->tecWheel(id) == 7) {
1198             mehSiStripADC[17]->Fill((*iter).adc());
1199             mehSiStripStrip[17]->Fill((*iter).strip());
1200           }
1201           if (tTopo->tecWheel(id) == 8) {
1202             mehSiStripADC[18]->Fill((*iter).adc());
1203             mehSiStripStrip[18]->Fill((*iter).strip());
1204           }
1205         }
1206       }
1207     }  // end loop over DataSetVector
1208 
1209     if (verbosity > 1) {
1210       eventout += "\n          Number of BrlStripDigis collected:........ ";
1211       eventout += nStripBrl;
1212     }
1213     for (int i = 0; i < 8; ++i) {
1214       mehSiStripn[i]->Fill((float)nStripBrl);
1215     }
1216 
1217     if (verbosity > 1) {
1218       eventout += "\n          Number of FrwdStripDigis collected:....... ";
1219       eventout += nStripFwd;
1220     }
1221     for (int i = 8; i < 19; ++i) {
1222       mehSiStripn[i]->Fill((float)nStripFwd);
1223     }
1224   }
1225 
1226   // get pixel information
1227   edm::Handle<edm::DetSetVector<PixelDigi>> pixelDigis;
1228   iEvent.getByToken(SiPxlSrc_Token_, pixelDigis);
1229   bool validpixelDigis = true;
1230   if (!pixelDigis.isValid()) {
1231     LogDebug(MsgLoggerCat) << "Unable to find pixelDigis in event!";
1232     validpixelDigis = false;
1233   }
1234   if (validpixelDigis) {
1235     int nPxlBrl = 0, nPxlFwd = 0;
1236     edm::DetSetVector<PixelDigi>::const_iterator DPViter;
1237     for (DPViter = pixelDigis->begin(); DPViter != pixelDigis->end(); ++DPViter) {
1238       unsigned int id = DPViter->id;
1239       DetId detId(id);
1240       edm::DetSet<PixelDigi>::const_iterator begin = DPViter->data.begin();
1241       edm::DetSet<PixelDigi>::const_iterator end = DPViter->data.end();
1242       edm::DetSet<PixelDigi>::const_iterator iter;
1243 
1244       // get Barrel pixels
1245       if (detId.subdetId() == sdPxlBrl) {
1246         for (iter = begin; iter != end; ++iter) {
1247           ++nPxlBrl;
1248           if (tTopo->pxbLayer(id) == 1) {
1249             mehSiPixelADC[0]->Fill((*iter).adc());
1250             mehSiPixelRow[0]->Fill((*iter).row());
1251             mehSiPixelCol[0]->Fill((*iter).column());
1252           }
1253           if (tTopo->pxbLayer(id) == 2) {
1254             mehSiPixelADC[1]->Fill((*iter).adc());
1255             mehSiPixelRow[1]->Fill((*iter).row());
1256             mehSiPixelCol[1]->Fill((*iter).column());
1257           }
1258           if (tTopo->pxbLayer(id) == 3) {
1259             mehSiPixelADC[2]->Fill((*iter).adc());
1260             mehSiPixelRow[2]->Fill((*iter).row());
1261             mehSiPixelCol[2]->Fill((*iter).column());
1262           }
1263         }
1264       }
1265 
1266       // get Forward pixels
1267       if (detId.subdetId() == sdPxlFwd) {
1268         for (iter = begin; iter != end; ++iter) {
1269           ++nPxlFwd;
1270           if (tTopo->pxfDisk(id) == 1) {
1271             if (tTopo->pxfSide(id) == 1) {
1272               mehSiPixelADC[4]->Fill((*iter).adc());
1273               mehSiPixelRow[4]->Fill((*iter).row());
1274               mehSiPixelCol[4]->Fill((*iter).column());
1275             }
1276             if (tTopo->pxfSide(id) == 2) {
1277               mehSiPixelADC[3]->Fill((*iter).adc());
1278               mehSiPixelRow[3]->Fill((*iter).row());
1279               mehSiPixelCol[3]->Fill((*iter).column());
1280             }
1281           }
1282           if (tTopo->pxfDisk(id) == 2) {
1283             if (tTopo->pxfSide(id) == 1) {
1284               mehSiPixelADC[6]->Fill((*iter).adc());
1285               mehSiPixelRow[6]->Fill((*iter).row());
1286               mehSiPixelCol[6]->Fill((*iter).column());
1287             }
1288             if (tTopo->pxfSide(id) == 2) {
1289               mehSiPixelADC[5]->Fill((*iter).adc());
1290               mehSiPixelRow[5]->Fill((*iter).row());
1291               mehSiPixelCol[5]->Fill((*iter).column());
1292             }
1293           }
1294         }
1295       }
1296     }
1297 
1298     if (verbosity > 1) {
1299       eventout += "\n          Number of BrlPixelDigis collected:........ ";
1300       eventout += nPxlBrl;
1301     }
1302     for (int i = 0; i < 3; ++i) {
1303       mehSiPixeln[i]->Fill((float)nPxlBrl);
1304     }
1305 
1306     if (verbosity > 1) {
1307       eventout += "\n          Number of FrwdPixelDigis collected:....... ";
1308       eventout += nPxlFwd;
1309     }
1310 
1311     for (int i = 3; i < 7; ++i) {
1312       mehSiPixeln[i]->Fill((float)nPxlFwd);
1313     }
1314   }
1315 
1316   if (verbosity > 0)
1317     edm::LogInfo(MsgLoggerCat) << eventout << "\n";
1318 
1319   return;
1320 }
1321 
1322 void GlobalDigisAnalyzer::fillMuon(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
1323   std::string MsgLoggerCat = "GlobalDigisAnalyzer_fillMuon";
1324 
1325   TString eventout;
1326   if (verbosity > 0)
1327     eventout = "\nGathering info:";
1328 
1329   // get DT information
1330   edm::Handle<DTDigiCollection> dtDigis;
1331   iEvent.getByToken(MuDTSrc_Token_, dtDigis);
1332   bool validdtDigis = true;
1333   if (!dtDigis.isValid()) {
1334     LogDebug(MsgLoggerCat) << "Unable to find dtDigis in event!";
1335     validdtDigis = false;
1336   }
1337   if (validdtDigis) {
1338     int nDt0 = 0;
1339     int nDt1 = 0;
1340     int nDt2 = 0;
1341     int nDt3 = 0;
1342     int nDt = 0;
1343     DTDigiCollection::DigiRangeIterator detUnitIt;
1344     for (detUnitIt = dtDigis->begin(); detUnitIt != dtDigis->end(); ++detUnitIt) {
1345       const DTLayerId &id = (*detUnitIt).first;
1346       const DTDigiCollection::Range &range = (*detUnitIt).second;
1347       for (DTDigiCollection::const_iterator digiIt = range.first; digiIt != range.second; ++digiIt) {
1348         ++nDt;
1349 
1350         DTWireId wireId(id, (*digiIt).wire());
1351         if (wireId.station() == 1) {
1352           mehDtMuonLayer[0]->Fill(id.layer());
1353           mehDtMuonTime[0]->Fill((*digiIt).time());
1354           mehDtMuonTimevLayer[0]->Fill(id.layer(), (*digiIt).time());
1355           ++nDt0;
1356         }
1357         if (wireId.station() == 2) {
1358           mehDtMuonLayer[1]->Fill(id.layer());
1359           mehDtMuonTime[1]->Fill((*digiIt).time());
1360           mehDtMuonTimevLayer[1]->Fill(id.layer(), (*digiIt).time());
1361           ++nDt1;
1362         }
1363         if (wireId.station() == 3) {
1364           mehDtMuonLayer[2]->Fill(id.layer());
1365           mehDtMuonTime[2]->Fill((*digiIt).time());
1366           mehDtMuonTimevLayer[2]->Fill(id.layer(), (*digiIt).time());
1367           ++nDt2;
1368         }
1369         if (wireId.station() == 4) {
1370           mehDtMuonLayer[3]->Fill(id.layer());
1371           mehDtMuonTime[3]->Fill((*digiIt).time());
1372           mehDtMuonTimevLayer[3]->Fill(id.layer(), (*digiIt).time());
1373           ++nDt3;
1374         }
1375       }
1376     }
1377     mehDtMuonn[0]->Fill((float)nDt0);
1378     mehDtMuonn[1]->Fill((float)nDt1);
1379     mehDtMuonn[2]->Fill((float)nDt2);
1380     mehDtMuonn[3]->Fill((float)nDt3);
1381 
1382     if (verbosity > 1) {
1383       eventout += "\n          Number of DtMuonDigis collected:.......... ";
1384       eventout += nDt;
1385     }
1386   }
1387 
1388   // get CSC Strip information
1389   edm::Handle<CSCStripDigiCollection> strips;
1390   iEvent.getByToken(MuCSCStripSrc_Token_, strips);
1391   bool validstrips = true;
1392   if (!strips.isValid()) {
1393     LogDebug(MsgLoggerCat) << "Unable to find muon strips in event!";
1394     validstrips = false;
1395   }
1396 
1397   if (validstrips) {
1398     int nStrips = 0;
1399     for (CSCStripDigiCollection::DigiRangeIterator j = strips->begin(); j != strips->end(); ++j) {
1400       std::vector<CSCStripDigi>::const_iterator digiItr = (*j).second.first;
1401       std::vector<CSCStripDigi>::const_iterator last = (*j).second.second;
1402 
1403       for (; digiItr != last; ++digiItr) {
1404         ++nStrips;
1405 
1406         // average pedestals
1407         std::vector<int> adcCounts = digiItr->getADCCounts();
1408         theCSCStripPedestalSum += adcCounts[0];
1409         theCSCStripPedestalSum += adcCounts[1];
1410         theCSCStripPedestalCount += 2;
1411 
1412         // if there are enough pedestal statistics
1413         if (theCSCStripPedestalCount > 100) {
1414           float pedestal = theCSCStripPedestalSum / theCSCStripPedestalCount;
1415           if (adcCounts[5] > (pedestal + 100))
1416             mehCSCStripADC->Fill(adcCounts[4] - pedestal);
1417         }
1418       }
1419     }
1420 
1421     if (verbosity > 1) {
1422       eventout += "\n          Number of CSCStripDigis collected:........ ";
1423       eventout += nStrips;
1424     }
1425     mehCSCStripn->Fill((float)nStrips);
1426   }
1427 
1428   // get CSC Wire information
1429   edm::Handle<CSCWireDigiCollection> wires;
1430   iEvent.getByToken(MuCSCWireSrc_Token_, wires);
1431   bool validwires = true;
1432   if (!wires.isValid()) {
1433     LogDebug(MsgLoggerCat) << "Unable to find muon wires in event!";
1434     validwires = false;
1435   }
1436 
1437   if (validwires) {
1438     int nWires = 0;
1439     for (CSCWireDigiCollection::DigiRangeIterator j = wires->begin(); j != wires->end(); ++j) {
1440       std::vector<CSCWireDigi>::const_iterator digiItr = (*j).second.first;
1441       std::vector<CSCWireDigi>::const_iterator endDigi = (*j).second.second;
1442 
1443       for (; digiItr != endDigi; ++digiItr) {
1444         ++nWires;
1445         mehCSCWireTime->Fill(digiItr->getTimeBin());
1446       }
1447     }
1448 
1449     if (verbosity > 1) {
1450       eventout += "\n          Number of CSCWireDigis collected:......... ";
1451       eventout += nWires;
1452     }
1453     mehCSCWiren->Fill((float)nWires);
1454   }
1455 
1456   // get RPC information
1457   // Get the RPC Geometry
1458   const auto &rpcGeom = iSetup.getHandle(rpcGeomToken_);
1459   if (!rpcGeom.isValid()) {
1460     edm::LogWarning(MsgLoggerCat) << "Unable to find RPCGeometryRecord in event!";
1461     return;
1462   }
1463 
1464   edm::Handle<edm::PSimHitContainer> rpcsimHit;
1465   iEvent.getByToken(RPCSimHit_Token_, rpcsimHit);
1466   bool validrpcsim = true;
1467   if (!rpcsimHit.isValid()) {
1468     LogDebug(MsgLoggerCat) << "Unable to find rpcsimHit in event!";
1469     validrpcsim = false;
1470   }
1471 
1472   edm::Handle<RPCDigiCollection> rpcDigis;
1473   iEvent.getByToken(MuRPCSrc_Token_, rpcDigis);
1474   bool validrpcdigi = true;
1475   if (!rpcDigis.isValid()) {
1476     LogDebug(MsgLoggerCat) << "Unable to find rpcDigis in event!";
1477     validrpcdigi = false;
1478   }
1479 
1480   // ONLY UNTIL PROBLEM WITH RPC DIGIS IS FIGURED OUT
1481   validrpcdigi = false;
1482 
1483   // Loop on simhits
1484   edm::PSimHitContainer::const_iterator rpcsimIt;
1485   std::map<RPCDetId, std::vector<double>> allsims;
1486 
1487   if (validrpcsim) {
1488     for (rpcsimIt = rpcsimHit->begin(); rpcsimIt != rpcsimHit->end(); rpcsimIt++) {
1489       RPCDetId Rsid = (RPCDetId)(*rpcsimIt).detUnitId();
1490       int ptype = rpcsimIt->particleType();
1491 
1492       if (ptype == 13 || ptype == -13) {
1493         std::vector<double> buff;
1494         if (allsims.find(Rsid) != allsims.end()) {
1495           buff = allsims[Rsid];
1496         }
1497         buff.push_back(rpcsimIt->localPosition().x());
1498         allsims[Rsid] = buff;
1499       }
1500     }
1501   }
1502 
1503   // CRASH HAPPENS SOMEWHERE HERE IN FOR DECLARATION
1504   // WHAT IS WRONG WITH rpcDigis?????
1505   if (validrpcdigi) {
1506     int nRPC = 0;
1507     RPCDigiCollection::DigiRangeIterator rpcdetUnitIt;
1508     for (rpcdetUnitIt = rpcDigis->begin(); rpcdetUnitIt != rpcDigis->end(); ++rpcdetUnitIt) {
1509       const RPCDetId Rsid = (*rpcdetUnitIt).first;
1510       const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(Rsid));
1511       const RPCDigiCollection::Range &range = (*rpcdetUnitIt).second;
1512 
1513       std::vector<double> sims;
1514       if (allsims.find(Rsid) != allsims.end()) {
1515         sims = allsims[Rsid];
1516       }
1517 
1518       int ndigi = 0;
1519       for (RPCDigiCollection::const_iterator rpcdigiIt = range.first; rpcdigiIt != range.second; ++rpcdigiIt) {
1520         ++ndigi;
1521         ++nRPC;
1522       }
1523 
1524       if (sims.size() == 1 && ndigi == 1) {
1525         double dis = roll->centreOfStrip(range.first->strip()).x() - sims[0];
1526 
1527         if (Rsid.region() == 0) {
1528           if (Rsid.ring() == -2)
1529             mehRPCRes[0]->Fill((float)dis);
1530           else if (Rsid.ring() == -1)
1531             mehRPCRes[1]->Fill((float)dis);
1532           else if (Rsid.ring() == 0)
1533             mehRPCRes[2]->Fill((float)dis);
1534           else if (Rsid.ring() == 1)
1535             mehRPCRes[3]->Fill((float)dis);
1536           else if (Rsid.ring() == 2)
1537             mehRPCRes[4]->Fill((float)dis);
1538         }
1539       }
1540     }
1541 
1542     if (verbosity > 1) {
1543       eventout += "\n          Number of RPCDigis collected:.............. ";
1544       eventout += nRPC;
1545     }
1546     mehRPCMuonn->Fill(float(nRPC));
1547   }
1548 
1549   if (verbosity > 0)
1550     edm::LogInfo(MsgLoggerCat) << eventout << "\n";
1551 
1552   return;
1553 }