GlobalDigisProducer.cc

Go to the documentation of this file.

#include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
#include "Geometry/Records/interface/TrackerTopologyRcd.h"
#include "Validation/GlobalDigis/interface/GlobalDigisProducer.h"

GlobalDigisProducer::GlobalDigisProducer(const edm::ParameterSet &iPSet)
    : fName(""),
      verbosity(0),
      frequency(0),
      label(""),
      getAllProvenances(false),
      printProvenanceInfo(false),
      theCSCStripPedestalSum(0),
      theCSCStripPedestalCount(0),
      count(0) {
  std::string MsgLoggerCat = "GlobalDigisProducer_GlobalDigisProducer";

  // get information from parameter set
  fName = iPSet.getUntrackedParameter<std::string>("Name");
  verbosity = iPSet.getUntrackedParameter<int>("Verbosity");
  frequency = iPSet.getUntrackedParameter<int>("Frequency");
  label = iPSet.getParameter<std::string>("Label");
  edm::ParameterSet m_Prov = iPSet.getParameter<edm::ParameterSet>("ProvenanceLookup");
  getAllProvenances = m_Prov.getUntrackedParameter<bool>("GetAllProvenances");
  printProvenanceInfo = m_Prov.getUntrackedParameter<bool>("PrintProvenanceInfo");

  // get Labels to use to extract information
  ECalEBSrc_ = iPSet.getParameter<edm::InputTag>("ECalEBSrc");
  ECalEESrc_ = iPSet.getParameter<edm::InputTag>("ECalEESrc");
  ECalESSrc_ = iPSet.getParameter<edm::InputTag>("ECalESSrc");
  HCalSrc_ = iPSet.getParameter<edm::InputTag>("HCalSrc");
  HCalDigi_ = iPSet.getParameter<edm::InputTag>("HCalDigi");
  SiStripSrc_ = iPSet.getParameter<edm::InputTag>("SiStripSrc");
  SiPxlSrc_ = iPSet.getParameter<edm::InputTag>("SiPxlSrc");
  MuDTSrc_ = iPSet.getParameter<edm::InputTag>("MuDTSrc");
  MuCSCStripSrc_ = iPSet.getParameter<edm::InputTag>("MuCSCStripSrc");
  MuCSCWireSrc_ = iPSet.getParameter<edm::InputTag>("MuCSCWireSrc");
  // fix for consumes
  ECalEBSrc_Token_ = consumes<EBDigiCollection>(iPSet.getParameter<edm::InputTag>("ECalEBSrc"));
  ECalEESrc_Token_ = consumes<EEDigiCollection>(iPSet.getParameter<edm::InputTag>("ECalEESrc"));
  ECalESSrc_Token_ = consumes<ESDigiCollection>(iPSet.getParameter<edm::InputTag>("ECalESSrc"));
  HCalSrc_Token_ = consumes<edm::PCaloHitContainer>(iPSet.getParameter<edm::InputTag>("HCalSrc"));
  HBHEDigi_Token_ = consumes<edm::SortedCollection<HBHEDataFrame>>(iPSet.getParameter<edm::InputTag>("HCalDigi"));
  HODigi_Token_ = consumes<edm::SortedCollection<HODataFrame>>(iPSet.getParameter<edm::InputTag>("HCalDigi"));
  HFDigi_Token_ = consumes<edm::SortedCollection<HFDataFrame>>(iPSet.getParameter<edm::InputTag>("HCalDigi"));
  SiStripSrc_Token_ = consumes<edm::DetSetVector<SiStripDigi>>(iPSet.getParameter<edm::InputTag>("SiStripSrc"));
  SiPxlSrc_Token_ = consumes<edm::DetSetVector<PixelDigi>>(iPSet.getParameter<edm::InputTag>("SiPxlSrc"));
  MuDTSrc_Token_ = consumes<DTDigiCollection>(iPSet.getParameter<edm::InputTag>("MuDTSrc"));
  MuCSCStripSrc_Token_ = consumes<CSCStripDigiCollection>(iPSet.getParameter<edm::InputTag>("MuCSCStripSrc"));
  MuCSCWireSrc_Token_ = consumes<CSCWireDigiCollection>(iPSet.getParameter<edm::InputTag>("MuCSCWireSrc"));
  //
  const std::string barrelHitsName("EcalHitsEB");
  const std::string endcapHitsName("EcalHitsEE");
  const std::string preshowerHitsName("EcalHitsES");
  EBHits_Token_ = consumes<CrossingFrame<PCaloHit>>(edm::InputTag(std::string("mix"), std::string("barrelHitsName")));
  EEHits_Token_ = consumes<CrossingFrame<PCaloHit>>(edm::InputTag(std::string("mix"), std::string("endcapHitsName")));
  ESHits_Token_ =
      consumes<CrossingFrame<PCaloHit>>(edm::InputTag(std::string("mix"), std::string("preshowerHitsName")));

  // use value of first digit to determine default output level (inclusive)
  // 0 is none, 1 is basic, 2 is fill output, 3 is gather output
  verbosity %= 10;

  // create persistent object
  produces<PGlobalDigi>(label);

  // print out Parameter Set information being used
  if (verbosity >= 0) {
    edm::LogInfo(MsgLoggerCat) << "\n===============================\n"
                               << "Initialized as EDProducer with parameter values:\n"
                               << "    Name          = " << fName << "\n"
                               << "    Verbosity     = " << verbosity << "\n"
                               << "    Frequency     = " << frequency << "\n"
                               << "    Label         = " << label << "\n"
                               << "    GetProv       = " << getAllProvenances << "\n"
                               << "    PrintProv     = " << printProvenanceInfo << "\n"
                               << "    ECalEBSrc     = " << ECalEBSrc_.label() << ":" << ECalEBSrc_.instance() << "\n"
                               << "    ECalEESrc     = " << ECalEESrc_.label() << ":" << ECalEESrc_.instance() << "\n"
                               << "    ECalESSrc     = " << ECalESSrc_.label() << ":" << ECalESSrc_.instance() << "\n"
                               << "    HCalSrc       = " << HCalSrc_.label() << ":" << HCalSrc_.instance() << "\n"
                               << "    HCalDigi       = " << HCalDigi_.label() << ":" << HCalDigi_.instance() << "\n"
                               << "    SiStripSrc    = " << SiStripSrc_.label() << ":" << SiStripSrc_.instance() << "\n"
                               << "    SiPixelSrc    = " << SiPxlSrc_.label() << ":" << SiPxlSrc_.instance() << "\n"
                               << "    MuDTSrc       = " << MuDTSrc_.label() << ":" << MuDTSrc_.instance() << "\n"
                               << "    MuCSCStripSrc = " << MuCSCStripSrc_.label() << ":" << MuCSCStripSrc_.instance()
                               << "\n"
                               << "    MuCSCWireSrc  = " << MuCSCWireSrc_.label() << ":" << MuCSCWireSrc_.instance()
                               << "\n"
                               << "===============================\n";
  }

  // set default constants
  // ECal
  ECalgainConv_[0] = 0.;
  ECalgainConv_[1] = 1.;
  ECalgainConv_[2] = 2.;
  ECalgainConv_[3] = 12.;
  ECalbarrelADCtoGeV_ = 0.035;
  ECalendcapADCtoGeV_ = 0.06;
}

GlobalDigisProducer::~GlobalDigisProducer() {}

void GlobalDigisProducer::beginJob(void) {
  std::string MsgLoggerCat = "GlobalDigisProducer_beginJob";

  //   // setup calorimeter constants from service
  //   edm::ESHandle<EcalADCToGeVConstant> pAgc;
  //   iSetup.get<EcalADCToGeVConstantRcd>().get(pAgc);
  //   const EcalADCToGeVConstant* agc = pAgc.product();

  EcalMGPAGainRatio *defaultRatios = new EcalMGPAGainRatio();

  ECalgainConv_[0] = 0.;
  ECalgainConv_[1] = 1.;
  ECalgainConv_[2] = defaultRatios->gain12Over6();
  ECalgainConv_[3] = ECalgainConv_[2] * (defaultRatios->gain6Over1());

  delete defaultRatios;

  //   ECalbarrelADCtoGeV_ = agc->getEBValue();
  //   ECalendcapADCtoGeV_ = agc->getEEValue();

  if (verbosity >= 0) {
    edm::LogInfo(MsgLoggerCat) << "Modified Calorimeter gain constants: g0 = " << ECalgainConv_[0]
                               << ", g1 = " << ECalgainConv_[1] << ", g2 = " << ECalgainConv_[2]
                               << ", g3 = " << ECalgainConv_[3];
    //     edm::LogInfo(MsgLoggerCat)
    //       << "Modified Calorimeter ADCtoGeV constants: barrel = "
    //       << ECalbarrelADCtoGeV_ << ", endcap = " << ECalendcapADCtoGeV_;
  }

  // clear storage vectors
  clear();
  return;
}

void GlobalDigisProducer::endJob() {
  std::string MsgLoggerCat = "GlobalDigisProducer_endJob";
  if (verbosity >= 0)
    edm::LogInfo(MsgLoggerCat) << "Terminating having processed " << count << " events.";
  return;
}

void GlobalDigisProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
  std::string MsgLoggerCat = "GlobalDigisProducer_produce";

  // keep track of number of events processed
  ++count;

  // THIS BLOCK MIGRATED HERE FROM beginJob:
  // setup calorimeter constants from service
  edm::ESHandle<EcalADCToGeVConstant> pAgc;
  iSetup.get<EcalADCToGeVConstantRcd>().get(pAgc);
  const EcalADCToGeVConstant *agc = pAgc.product();
  ECalbarrelADCtoGeV_ = agc->getEBValue();
  ECalendcapADCtoGeV_ = agc->getEEValue();
  if (verbosity >= 0) {
    edm::LogInfo(MsgLoggerCat) << "Modified Calorimeter ADCtoGeV constants: barrel = " << ECalbarrelADCtoGeV_
                               << ", endcap = " << ECalendcapADCtoGeV_;
  }

  // get event id information
  edm::RunNumber_t nrun = iEvent.id().run();
  edm::EventNumber_t nevt = iEvent.id().event();

  if (verbosity > 0) {
    edm::LogInfo(MsgLoggerCat) << "Processing run " << nrun << ", event " << nevt << " (" << count << " events total)";
  } else if (verbosity == 0) {
    if (nevt % frequency == 0 || nevt == 1) {
      edm::LogInfo(MsgLoggerCat) << "Processing run " << nrun << ", event " << nevt << " (" << count
                                 << " events total)";
    }
  }

  // clear event holders
  clear();

  // look at information available in the event
  if (getAllProvenances) {
    std::vector<const edm::StableProvenance *> AllProv;
    iEvent.getAllStableProvenance(AllProv);

    if (verbosity >= 0)
      edm::LogInfo(MsgLoggerCat) << "Number of Provenances = " << AllProv.size();

    if (printProvenanceInfo && (verbosity >= 0)) {
      TString eventout("\nProvenance info:\n");

      for (unsigned int i = 0; i < AllProv.size(); ++i) {
        eventout += "\n       ******************************";
        eventout += "\n       Module       : ";
        // eventout += (AllProv[i]->product).moduleLabel();
        eventout += AllProv[i]->moduleLabel();
        eventout += "\n       ProductID    : ";
        // eventout += (AllProv[i]->product).productID_.id_;
        eventout += AllProv[i]->productID().id();
        eventout += "\n       ClassName    : ";
        // eventout += (AllProv[i]->product).fullClassName_;
        eventout += AllProv[i]->className();
        eventout += "\n       InstanceName : ";
        // eventout += (AllProv[i]->product).productInstanceName_;
        eventout += AllProv[i]->productInstanceName();
        eventout += "\n       BranchName   : ";
        // eventout += (AllProv[i]->product).branchName_;
        eventout += AllProv[i]->branchName();
      }
      eventout += "\n       ******************************\n";
      edm::LogInfo(MsgLoggerCat) << eventout << "\n";
      printProvenanceInfo = false;
    }
    getAllProvenances = false;
  }

  // call fill functions
  // gather Ecal information from event
  fillECal(iEvent, iSetup);
  // gather Hcal information from event
  fillHCal(iEvent, iSetup);
  // gather Track information from event
  fillTrk(iEvent, iSetup);
  // gather Muon information from event
  fillMuon(iEvent, iSetup);

  if (verbosity > 0)
    edm::LogInfo(MsgLoggerCat) << "Done gathering data from event.";

  // produce object to put into event
  std::unique_ptr<PGlobalDigi> pOut(new PGlobalDigi);

  if (verbosity > 2)
    edm::LogInfo(MsgLoggerCat) << "Saving event contents:";

  // call store functions
  // store ECal information in produce
  storeECal(*pOut);
  // store HCal information in produce
  storeHCal(*pOut);
  // store Track information in produce
  storeTrk(*pOut);
  // store Muon information in produce
  storeMuon(*pOut);

  // store information in event
  iEvent.put(std::move(pOut), label);

  return;
}

void GlobalDigisProducer::fillECal(edm::Event &iEvent, const edm::EventSetup &iSetup) {
  std::string MsgLoggerCat = "GlobalDigisProducer_fillECal";

  TString eventout;
  if (verbosity > 0)
    eventout = "\nGathering info:";

  // extract crossing frame from event
  // edm::Handle<CrossingFrame> crossingFrame;
  edm::Handle<CrossingFrame<PCaloHit>> crossingFrame;
  // iEvent.getByType(crossingFrame);
  // if (!crossingFrame.isValid()) {
  //  edm::LogWarning(MsgLoggerCat)
  //    << "Unable to crossingFrame in event!";
  //  return;
  //}

  // extract EB information
  bool isBarrel = true;
  edm::Handle<EBDigiCollection> EcalDigiEB;
  iEvent.getByToken(ECalEBSrc_Token_, EcalDigiEB);
  if (!EcalDigiEB.isValid()) {
    edm::LogWarning(MsgLoggerCat) << "Unable to find EcalDigiEB in event!";
    return;
  }
  if (EcalDigiEB->empty())
    isBarrel = false;

  if (isBarrel) {
    // loop over simhits
    iEvent.getByToken(EBHits_Token_, crossingFrame);
    if (!crossingFrame.isValid()) {
      edm::LogWarning(MsgLoggerCat) << "Unable to find cal barrel crossingFrame in event!";
      return;
    }
    // std::unique_ptr<MixCollection<PCaloHit> >
    // barrelHits(new MixCollection<PCaloHit>
    //       (crossingFrame.product(), barrelHitsName));
    std::unique_ptr<MixCollection<PCaloHit>> barrelHits(new MixCollection<PCaloHit>(crossingFrame.product()));

    // keep track of sum of simhit energy in each crystal
    MapType ebSimMap;
    for (MixCollection<PCaloHit>::MixItr hitItr = barrelHits->begin(); hitItr != barrelHits->end(); ++hitItr) {
      EBDetId ebid = EBDetId(hitItr->id());

      uint32_t crystid = ebid.rawId();
      ebSimMap[crystid] += hitItr->energy();
    }

    // loop over digis
    const EBDigiCollection *barrelDigi = EcalDigiEB.product();

    std::vector<double> ebAnalogSignal;
    std::vector<double> ebADCCounts;
    std::vector<double> ebADCGains;
    ebAnalogSignal.reserve(EBDataFrame::MAXSAMPLES);
    ebADCCounts.reserve(EBDataFrame::MAXSAMPLES);
    ebADCGains.reserve(EBDataFrame::MAXSAMPLES);

    int i = 0;
    for (unsigned int digis = 0; digis < EcalDigiEB->size(); ++digis) {
      // for (std::vector<EBDataFrame>::const_iterator digis =
      //   barrelDigi->begin();
      // digis != barrelDigi->end();
      // ++digis) {

      ++i;

      EBDataFrame ebdf = (*barrelDigi)[digis];
      int nrSamples = ebdf.size();

      EBDetId ebid = ebdf.id();
      // EBDetId ebid = digis->id();

      double Emax = 0;
      int Pmax = 0;
      double pedestalPreSample = 0.;
      double pedestalPreSampleAnalog = 0.;

      for (int sample = 0; sample < nrSamples; ++sample) {
        // for (int sample = 0; sample < digis->size(); ++sample) {
        ebAnalogSignal[sample] = 0.;
        ebADCCounts[sample] = 0.;
        ebADCGains[sample] = -1.;
      }

      // calculate maximum energy and pedestal
      for (int sample = 0; sample < nrSamples; ++sample) {
        // for (int sample = 0; sample < digis->size(); ++sample) {

        EcalMGPASample thisSample = ebdf[sample];
        ebADCCounts[sample] = (thisSample.adc());
        ebADCGains[sample] = (thisSample.gainId());
        ebAnalogSignal[sample] = (ebADCCounts[sample] * ECalgainConv_[(int)ebADCGains[sample]] * ECalbarrelADCtoGeV_);
        if (Emax < ebAnalogSignal[sample]) {
          Emax = ebAnalogSignal[sample];
          Pmax = sample;
        }
        if (sample < 3) {
          pedestalPreSample += ebADCCounts[sample];
          pedestalPreSampleAnalog += ebADCCounts[sample] * ECalgainConv_[(int)ebADCGains[sample]] * ECalbarrelADCtoGeV_;
        }
      }
      pedestalPreSample /= 3.;
      pedestalPreSampleAnalog /= 3.;

      // calculate pedestal subtracted digi energy in the crystal
      double Erec = Emax - pedestalPreSampleAnalog * ECalgainConv_[(int)ebADCGains[Pmax]];

      // gather necessary information
      EBCalAEE.push_back(Erec);
      EBCalSHE.push_back(ebSimMap[ebid.rawId()]);
      EBCalmaxPos.push_back(Pmax);
    }

    if (verbosity > 1) {
      eventout += "\n          Number of EBDigis collected:.............. ";
      eventout += i;
    }
  }

  // extract EE information
  bool isEndCap = true;
  edm::Handle<EEDigiCollection> EcalDigiEE;
  iEvent.getByToken(ECalEESrc_Token_, EcalDigiEE);
  if (!EcalDigiEE.isValid()) {
    edm::LogWarning(MsgLoggerCat) << "Unable to find EcalDigiEE in event!";
    return;
  }
  if (EcalDigiEE->empty())
    isEndCap = false;

  if (isEndCap) {
    // loop over simhits
    iEvent.getByToken(EEHits_Token_, crossingFrame);
    if (!crossingFrame.isValid()) {
      edm::LogWarning(MsgLoggerCat) << "Unable to find cal endcap crossingFrame in event!";
      return;
    }
    // std::unique_ptr<MixCollection<PCaloHit> >
    //  endcapHits(new MixCollection<PCaloHit>
    //   (crossingFrame.product(), endcapHitsName));
    std::unique_ptr<MixCollection<PCaloHit>> endcapHits(new MixCollection<PCaloHit>(crossingFrame.product()));

    // keep track of sum of simhit energy in each crystal
    MapType eeSimMap;
    for (MixCollection<PCaloHit>::MixItr hitItr = endcapHits->begin(); hitItr != endcapHits->end(); ++hitItr) {
      EEDetId eeid = EEDetId(hitItr->id());

      uint32_t crystid = eeid.rawId();
      eeSimMap[crystid] += hitItr->energy();
    }

    // loop over digis
    const EEDigiCollection *endcapDigi = EcalDigiEE.product();

    std::vector<double> eeAnalogSignal;
    std::vector<double> eeADCCounts;
    std::vector<double> eeADCGains;
    eeAnalogSignal.reserve(EEDataFrame::MAXSAMPLES);
    eeADCCounts.reserve(EEDataFrame::MAXSAMPLES);
    eeADCGains.reserve(EEDataFrame::MAXSAMPLES);

    int i = 0;
    // for (std::vector<EEDataFrame>::const_iterator digis =
    //   endcapDigi->begin();
    // digis != endcapDigi->end();
    // ++digis) {
    for (unsigned int digis = 0; digis < EcalDigiEE->size(); ++digis) {
      ++i;

      EEDataFrame eedf = (*endcapDigi)[digis];
      int nrSamples = eedf.size();

      EEDetId eeid = eedf.id();
      // EEDetId eeid = digis->id();

      double Emax = 0;
      int Pmax = 0;
      double pedestalPreSample = 0.;
      double pedestalPreSampleAnalog = 0.;

      for (int sample = 0; sample < nrSamples; ++sample) {
        // for (int sample = 0; sample < digis->size(); ++sample) {
        eeAnalogSignal[sample] = 0.;
        eeADCCounts[sample] = 0.;
        eeADCGains[sample] = -1.;
      }

      // calculate maximum enery and pedestal
      for (int sample = 0; sample < nrSamples; ++sample) {
        // for (int sample = 0; sample < digis->size(); ++sample) {

        EcalMGPASample thisSample = eedf[sample];

        eeADCCounts[sample] = (thisSample.adc());
        eeADCGains[sample] = (thisSample.gainId());
        eeAnalogSignal[sample] = (eeADCCounts[sample] * ECalgainConv_[(int)eeADCGains[sample]] * ECalbarrelADCtoGeV_);
        if (Emax < eeAnalogSignal[sample]) {
          Emax = eeAnalogSignal[sample];
          Pmax = sample;
        }
        if (sample < 3) {
          pedestalPreSample += eeADCCounts[sample];
          pedestalPreSampleAnalog += eeADCCounts[sample] * ECalgainConv_[(int)eeADCGains[sample]] * ECalbarrelADCtoGeV_;
        }
      }
      pedestalPreSample /= 3.;
      pedestalPreSampleAnalog /= 3.;

      // calculate pedestal subtracted digi energy in the crystal
      double Erec = Emax - pedestalPreSampleAnalog * ECalgainConv_[(int)eeADCGains[Pmax]];

      // gather necessary information
      EECalAEE.push_back(Erec);
      EECalSHE.push_back(eeSimMap[eeid.rawId()]);
      EECalmaxPos.push_back(Pmax);
    }

    if (verbosity > 1) {
      eventout += "\n          Number of EEDigis collected:.............. ";
      eventout += i;
    }
  }

  // extract ES information
  bool isPreshower = true;
  edm::Handle<ESDigiCollection> EcalDigiES;
  iEvent.getByToken(ECalESSrc_Token_, EcalDigiES);
  if (!EcalDigiES.isValid()) {
    edm::LogWarning(MsgLoggerCat) << "Unable to find EcalDigiES in event!";
    return;
  }
  if (EcalDigiES->empty())
    isPreshower = false;

  if (isPreshower) {
    // loop over simhits
    iEvent.getByToken(ESHits_Token_, crossingFrame);
    if (!crossingFrame.isValid()) {
      edm::LogWarning(MsgLoggerCat) << "Unable to find cal preshower crossingFrame in event!";
      return;
    }
    // std::unique_ptr<MixCollection<PCaloHit> >
    //  preshowerHits(new MixCollection<PCaloHit>
    //       (crossingFrame.product(), preshowerHitsName));
    std::unique_ptr<MixCollection<PCaloHit>> preshowerHits(new MixCollection<PCaloHit>(crossingFrame.product()));

    // keep track of sum of simhit energy in each crystal
    MapType esSimMap;
    for (MixCollection<PCaloHit>::MixItr hitItr = preshowerHits->begin(); hitItr != preshowerHits->end(); ++hitItr) {
      ESDetId esid = ESDetId(hitItr->id());

      uint32_t crystid = esid.rawId();
      esSimMap[crystid] += hitItr->energy();
    }

    // loop over digis
    const ESDigiCollection *preshowerDigi = EcalDigiES.product();

    std::vector<double> esADCCounts;
    esADCCounts.reserve(ESDataFrame::MAXSAMPLES);

    int i = 0;
    for (unsigned int digis = 0; digis < EcalDigiES->size(); ++digis) {
      // for (std::vector<ESDataFrame>::const_iterator digis =
      //   preshowerDigi->begin();
      // digis != preshowerDigi->end();
      // ++digis) {

      ++i;

      ESDataFrame esdf = (*preshowerDigi)[digis];
      int nrSamples = esdf.size();

      ESDetId esid = esdf.id();
      // ESDetId esid = digis->id();

      for (int sample = 0; sample < nrSamples; ++sample) {
        // for (int sample = 0; sample < digis->size(); ++sample) {
        esADCCounts[sample] = 0.;
      }

      // gether ADC counts
      for (int sample = 0; sample < nrSamples; ++sample) {
        ESSample thisSample = esdf[sample];
        // for (int sample = 0; sample < digis->size(); ++sample) {
        esADCCounts[sample] = (thisSample.adc());
      }

      ESCalADC0.push_back(esADCCounts[0]);
      ESCalADC1.push_back(esADCCounts[1]);
      ESCalADC2.push_back(esADCCounts[2]);
      ESCalSHE.push_back(esSimMap[esid.rawId()]);
    }

    if (verbosity > 1) {
      eventout += "\n          Number of ESDigis collected:.............. ";
      eventout += i;
    }
  }

  if (verbosity > 0)
    edm::LogInfo(MsgLoggerCat) << eventout << "\n";

  return;
}

void GlobalDigisProducer::storeECal(PGlobalDigi &product) {
  std::string MsgLoggerCat = "GlobalDigisProducer_storeECal";

  if (verbosity > 2) {
    TString eventout("\n         nEBDigis     = ");
    eventout += EBCalmaxPos.size();
    for (unsigned int i = 0; i < EBCalmaxPos.size(); ++i) {
      eventout += "\n      (maxPos, AEE, SHE) = (";
      eventout += EBCalmaxPos[i];
      eventout += ", ";
      eventout += EBCalAEE[i];
      eventout += ", ";
      eventout += EBCalSHE[i];
      eventout += ")";
    }
    eventout += "\n         nEEDigis     = ";
    eventout += EECalmaxPos.size();
    for (unsigned int i = 0; i < EECalmaxPos.size(); ++i) {
      eventout += "\n      (maxPos, AEE, SHE) = (";
      eventout += EECalmaxPos[i];
      eventout += ", ";
      eventout += EECalAEE[i];
      eventout += ", ";
      eventout += EECalSHE[i];
      eventout += ")";
    }
    eventout += "\n         nESDigis          = ";
    eventout += ESCalADC0.size();
    for (unsigned int i = 0; i < ESCalADC0.size(); ++i) {
      eventout += "\n      (ADC0, ADC1, ADC2, SHE) = (";
      eventout += ESCalADC0[i];
      eventout += ", ";
      eventout += ESCalADC1[i];
      eventout += ", ";
      eventout += ESCalADC2[i];
      eventout += ", ";
      eventout += ESCalSHE[i];
      eventout += ")";
    }
    edm::LogInfo(MsgLoggerCat) << eventout << "\n";
  }

  product.putEBCalDigis(EBCalmaxPos, EBCalAEE, EBCalSHE);
  product.putEECalDigis(EECalmaxPos, EECalAEE, EECalSHE);
  product.putESCalDigis(ESCalADC0, ESCalADC1, ESCalADC2, ESCalSHE);

  return;
}

void GlobalDigisProducer::fillHCal(edm::Event &iEvent, const edm::EventSetup &iSetup) {
  std::string MsgLoggerCat = "GlobalDigisProducer_fillHCal";

  TString eventout;
  if (verbosity > 0)
    eventout = "\nGathering info:";

  // get calibration info
  edm::ESHandle<HcalDbService> HCalconditions;
  iSetup.get<HcalDbRecord>().get(HCalconditions);
  if (!HCalconditions.isValid()) {
    edm::LogWarning(MsgLoggerCat) << "Unable to find HCalconditions in event!";
    return;
  }
  // HcalCalibrations calibrations;
  CaloSamples tool;

  // extract simhit info
  edm::Handle<edm::PCaloHitContainer> hcalHits;
  iEvent.getByToken(HCalSrc_Token_, hcalHits);
  if (!hcalHits.isValid()) {
    edm::LogWarning(MsgLoggerCat) << "Unable to find hcalHits in event!";
    return;
  }
  const edm::PCaloHitContainer *simhitResult = hcalHits.product();

  MapType fHBEnergySimHits;
  MapType fHEEnergySimHits;
  MapType fHOEnergySimHits;
  MapType fHFEnergySimHits;
  for (std::vector<PCaloHit>::const_iterator simhits = simhitResult->begin(); simhits != simhitResult->end();
       ++simhits) {
    HcalDetId detId(simhits->id());
    uint32_t cellid = detId.rawId();

    if (detId.subdet() == sdHcalBrl) {
      fHBEnergySimHits[cellid] += simhits->energy();
    }
    if (detId.subdet() == sdHcalEC) {
      fHEEnergySimHits[cellid] += simhits->energy();
    }
    if (detId.subdet() == sdHcalOut) {
      fHOEnergySimHits[cellid] += simhits->energy();
    }
    if (detId.subdet() == sdHcalFwd) {
      fHFEnergySimHits[cellid] += simhits->energy();
    }
  }

  // get HBHE information
  edm::Handle<edm::SortedCollection<HBHEDataFrame>> hbhe;
  iEvent.getByToken(HBHEDigi_Token_, hbhe);
  if (!hbhe.isValid()) {
    edm::LogWarning(MsgLoggerCat) << "Unable to find HBHEDataFrame in event!";
    return;
  }
  edm::SortedCollection<HBHEDataFrame>::const_iterator ihbhe;

  int iHB = 0;
  int iHE = 0;
  for (ihbhe = hbhe->begin(); ihbhe != hbhe->end(); ++ihbhe) {
    HcalDetId cell(ihbhe->id());

    if ((cell.subdet() == sdHcalBrl) || (cell.subdet() == sdHcalEC)) {
      // HCalconditions->makeHcalCalibration(cell, &calibrations);
      const HcalCalibrations &calibrations = HCalconditions->getHcalCalibrations(cell);
      const HcalQIECoder *channelCoder = HCalconditions->getHcalCoder(cell);
      const HcalQIEShape *shape = HCalconditions->getHcalShape(channelCoder);

      HcalCoderDb coder(*channelCoder, *shape);
      coder.adc2fC(*ihbhe, tool);

      // get HB info
      if (cell.subdet() == sdHcalBrl) {
        ++iHB;
        float fDigiSum = 0.0;
        for (int ii = 0; ii < tool.size(); ++ii) {
          // default ped is 4.5
          int capid = (*ihbhe)[ii].capid();
          fDigiSum += (tool[ii] - calibrations.pedestal(capid));
        }

        HBCalAEE.push_back(fDigiSum);
        HBCalSHE.push_back(fHBEnergySimHits[cell.rawId()]);
      }

      // get HE info
      if (cell.subdet() == sdHcalEC) {
        ++iHE;
        float fDigiSum = 0.0;
        for (int ii = 0; ii < tool.size(); ++ii) {
          int capid = (*ihbhe)[ii].capid();
          fDigiSum += (tool[ii] - calibrations.pedestal(capid));
        }

        HECalAEE.push_back(fDigiSum);
        HECalSHE.push_back(fHEEnergySimHits[cell.rawId()]);
      }
    }
  }

  if (verbosity > 1) {
    eventout += "\n          Number of HBDigis collected:.............. ";
    eventout += iHB;
  }

  if (verbosity > 1) {
    eventout += "\n          Number of HEDigis collected:.............. ";
    eventout += iHE;
  }

  // get HO information
  edm::Handle<edm::SortedCollection<HODataFrame>> ho;
  iEvent.getByToken(HODigi_Token_, ho);
  if (!ho.isValid()) {
    edm::LogWarning(MsgLoggerCat) << "Unable to find HODataFrame in event!";
    return;
  }
  edm::SortedCollection<HODataFrame>::const_iterator iho;

  int iHO = 0;
  for (iho = ho->begin(); iho != ho->end(); ++iho) {
    HcalDetId cell(iho->id());

    if (cell.subdet() == sdHcalOut) {
      // HCalconditions->makeHcalCalibration(cell, &calibrations);
      const HcalCalibrations &calibrations = HCalconditions->getHcalCalibrations(cell);
      const HcalQIECoder *channelCoder = HCalconditions->getHcalCoder(cell);
      const HcalQIEShape *shape = HCalconditions->getHcalShape(channelCoder);

      HcalCoderDb coder(*channelCoder, *shape);
      coder.adc2fC(*iho, tool);

      ++iHO;
      float fDigiSum = 0.0;
      for (int ii = 0; ii < tool.size(); ++ii) {
        // default ped is 4.5
        int capid = (*iho)[ii].capid();
        fDigiSum += (tool[ii] - calibrations.pedestal(capid));
      }

      HOCalAEE.push_back(fDigiSum);
      HOCalSHE.push_back(fHOEnergySimHits[cell.rawId()]);
    }
  }

  if (verbosity > 1) {
    eventout += "\n          Number of HODigis collected:.............. ";
    eventout += iHO;
  }

  // get HF information
  edm::Handle<edm::SortedCollection<HFDataFrame>> hf;
  iEvent.getByToken(HFDigi_Token_, hf);
  if (!hf.isValid()) {
    edm::LogWarning(MsgLoggerCat) << "Unable to find HFDataFrame in event!";
    return;
  }
  edm::SortedCollection<HFDataFrame>::const_iterator ihf;

  int iHF = 0;
  for (ihf = hf->begin(); ihf != hf->end(); ++ihf) {
    HcalDetId cell(ihf->id());

    if (cell.subdet() == sdHcalFwd) {
      // HCalconditions->makeHcalCalibration(cell, &calibrations);
      const HcalCalibrations &calibrations = HCalconditions->getHcalCalibrations(cell);
      const HcalQIECoder *channelCoder = HCalconditions->getHcalCoder(cell);
      const HcalQIEShape *shape = HCalconditions->getHcalShape(channelCoder);

      HcalCoderDb coder(*channelCoder, *shape);
      coder.adc2fC(*ihf, tool);

      ++iHF;
      float fDigiSum = 0.0;
      for (int ii = 0; ii < tool.size(); ++ii) {
        // default ped is 1.73077
        int capid = (*ihf)[ii].capid();
        fDigiSum += (tool[ii] - calibrations.pedestal(capid));
      }

      HFCalAEE.push_back(fDigiSum);
      HFCalSHE.push_back(fHFEnergySimHits[cell.rawId()]);
    }
  }

  if (verbosity > 1) {
    eventout += "\n          Number of HFDigis collected:.............. ";
    eventout += iHF;
  }

  if (verbosity > 0)
    edm::LogInfo(MsgLoggerCat) << eventout << "\n";

  return;
}

void GlobalDigisProducer::storeHCal(PGlobalDigi &product) {
  std::string MsgLoggerCat = "GlobalDigisProducer_storeHCal";

  if (verbosity > 2) {
    TString eventout("\n         nHBDigis     = ");
    eventout += HBCalAEE.size();
    for (unsigned int i = 0; i < HBCalAEE.size(); ++i) {
      eventout += "\n      (AEE, SHE) = (";
      eventout += HBCalAEE[i];
      eventout += ", ";
      eventout += HBCalSHE[i];
      eventout += ")";
    }
    eventout += "\n         nHEDigis     = ";
    eventout += HECalAEE.size();
    for (unsigned int i = 0; i < HECalAEE.size(); ++i) {
      eventout += "\n      (AEE, SHE) = (";
      eventout += HECalAEE[i];
      eventout += ", ";
      eventout += HECalSHE[i];
      eventout += ")";
    }
    eventout += "\n         nHODigis     = ";
    eventout += HOCalAEE.size();
    for (unsigned int i = 0; i < HOCalAEE.size(); ++i) {
      eventout += "\n      (AEE, SHE) = (";
      eventout += HOCalAEE[i];
      eventout += ", ";
      eventout += HOCalSHE[i];
      eventout += ")";
    }
    eventout += "\n         nHFDigis     = ";
    eventout += HFCalAEE.size();
    for (unsigned int i = 0; i < HFCalAEE.size(); ++i) {
      eventout += "\n      (AEE, SHE) = (";
      eventout += HFCalAEE[i];
      eventout += ", ";
      eventout += HFCalSHE[i];
      eventout += ")";
    }

    edm::LogInfo(MsgLoggerCat) << eventout << "\n";
  }

  product.putHBCalDigis(HBCalAEE, HBCalSHE);
  product.putHECalDigis(HECalAEE, HECalSHE);
  product.putHOCalDigis(HOCalAEE, HOCalSHE);
  product.putHFCalDigis(HFCalAEE, HFCalSHE);

  return;
}

void GlobalDigisProducer::fillTrk(edm::Event &iEvent, const edm::EventSetup &iSetup) {
  // Retrieve tracker topology from geometry
  edm::ESHandle<TrackerTopology> tTopoHandle;
  iSetup.get<TrackerTopologyRcd>().get(tTopoHandle);
  const TrackerTopology *const tTopo = tTopoHandle.product();

  std::string MsgLoggerCat = "GlobalDigisProducer_fillTrk";

  TString eventout;
  if (verbosity > 0)
    eventout = "\nGathering info:";

  // get strip information
  edm::Handle<edm::DetSetVector<SiStripDigi>> stripDigis;
  iEvent.getByToken(SiStripSrc_Token_, stripDigis);
  if (!stripDigis.isValid()) {
    edm::LogWarning(MsgLoggerCat) << "Unable to find stripDigis in event!";
    return;
  }

  int nStripBrl = 0, nStripFwd = 0;
  edm::DetSetVector<SiStripDigi>::const_iterator DSViter;
  for (DSViter = stripDigis->begin(); DSViter != stripDigis->end(); ++DSViter) {
    unsigned int id = DSViter->id;
    DetId detId(id);
    edm::DetSet<SiStripDigi>::const_iterator begin = DSViter->data.begin();
    edm::DetSet<SiStripDigi>::const_iterator end = DSViter->data.end();
    edm::DetSet<SiStripDigi>::const_iterator iter;

    // get TIB
    if (detId.subdetId() == sdSiTIB) {
      for (iter = begin; iter != end; ++iter) {
        ++nStripBrl;
        if (tTopo->tibLayer(id) == 1) {
          TIBL1ADC.push_back((*iter).adc());
          TIBL1Strip.push_back((*iter).strip());
        }
        if (tTopo->tibLayer(id) == 2) {
          TIBL2ADC.push_back((*iter).adc());
          TIBL2Strip.push_back((*iter).strip());
        }
        if (tTopo->tibLayer(id) == 3) {
          TIBL3ADC.push_back((*iter).adc());
          TIBL3Strip.push_back((*iter).strip());
        }
        if (tTopo->tibLayer(id) == 4) {
          TIBL4ADC.push_back((*iter).adc());
          TIBL4Strip.push_back((*iter).strip());
        }
      }
    }

    // get TOB
    if (detId.subdetId() == sdSiTOB) {
      for (iter = begin; iter != end; ++iter) {
        ++nStripBrl;
        if (tTopo->tobLayer(id) == 1) {
          TOBL1ADC.push_back((*iter).adc());
          TOBL1Strip.push_back((*iter).strip());
        }
        if (tTopo->tobLayer(id) == 2) {
          TOBL2ADC.push_back((*iter).adc());
          TOBL2Strip.push_back((*iter).strip());
        }
        if (tTopo->tobLayer(id) == 3) {
          TOBL3ADC.push_back((*iter).adc());
          TOBL3Strip.push_back((*iter).strip());
        }
        if (tTopo->tobLayer(id) == 4) {
          TOBL4ADC.push_back((*iter).adc());
          TOBL4Strip.push_back((*iter).strip());
        }
      }
    }

    // get TID
    if (detId.subdetId() == sdSiTID) {
      for (iter = begin; iter != end; ++iter) {
        ++nStripFwd;
        if (tTopo->tidWheel(id) == 1) {
          TIDW1ADC.push_back((*iter).adc());
          TIDW1Strip.push_back((*iter).strip());
        }
        if (tTopo->tidWheel(id) == 2) {
          TIDW2ADC.push_back((*iter).adc());
          TIDW2Strip.push_back((*iter).strip());
        }
        if (tTopo->tidWheel(id) == 3) {
          TIDW3ADC.push_back((*iter).adc());
          TIDW3Strip.push_back((*iter).strip());
        }
      }
    }

    // get TEC
    if (detId.subdetId() == sdSiTEC) {
      for (iter = begin; iter != end; ++iter) {
        ++nStripFwd;
        if (tTopo->tecWheel(id) == 1) {
          TECW1ADC.push_back((*iter).adc());
          TECW1Strip.push_back((*iter).strip());
        }
        if (tTopo->tecWheel(id) == 2) {
          TECW2ADC.push_back((*iter).adc());
          TECW2Strip.push_back((*iter).strip());
        }
        if (tTopo->tecWheel(id) == 3) {
          TECW3ADC.push_back((*iter).adc());
          TECW3Strip.push_back((*iter).strip());
        }
        if (tTopo->tecWheel(id) == 4) {
          TECW4ADC.push_back((*iter).adc());
          TECW4Strip.push_back((*iter).strip());
        }
        if (tTopo->tecWheel(id) == 5) {
          TECW5ADC.push_back((*iter).adc());
          TECW5Strip.push_back((*iter).strip());
        }
        if (tTopo->tecWheel(id) == 6) {
          TECW6ADC.push_back((*iter).adc());
          TECW6Strip.push_back((*iter).strip());
        }
        if (tTopo->tecWheel(id) == 7) {
          TECW7ADC.push_back((*iter).adc());
          TECW7Strip.push_back((*iter).strip());
        }
        if (tTopo->tecWheel(id) == 8) {
          TECW8ADC.push_back((*iter).adc());
          TECW8Strip.push_back((*iter).strip());
        }
      }
    }
  }  // end loop over DataSetVector

  if (verbosity > 1) {
    eventout += "\n          Number of BrlStripDigis collected:........ ";
    eventout += nStripBrl;
  }

  if (verbosity > 1) {
    eventout += "\n          Number of FrwdStripDigis collected:....... ";
    eventout += nStripFwd;
  }

  // get pixel information
  edm::Handle<edm::DetSetVector<PixelDigi>> pixelDigis;
  iEvent.getByToken(SiPxlSrc_Token_, pixelDigis);
  if (!pixelDigis.isValid()) {
    edm::LogWarning(MsgLoggerCat) << "Unable to find pixelDigis in event!";
    return;
  }

  int nPxlBrl = 0, nPxlFwd = 0;
  edm::DetSetVector<PixelDigi>::const_iterator DPViter;
  for (DPViter = pixelDigis->begin(); DPViter != pixelDigis->end(); ++DPViter) {
    unsigned int id = DPViter->id;
    DetId detId(id);
    edm::DetSet<PixelDigi>::const_iterator begin = DPViter->data.begin();
    edm::DetSet<PixelDigi>::const_iterator end = DPViter->data.end();
    edm::DetSet<PixelDigi>::const_iterator iter;

    // get Barrel pixels
    if (detId.subdetId() == sdPxlBrl) {
      for (iter = begin; iter != end; ++iter) {
        ++nPxlBrl;
        if (tTopo->pxbLayer(id) == 1) {
          BRL1ADC.push_back((*iter).adc());
          BRL1Row.push_back((*iter).row());
          BRL1Col.push_back((*iter).column());
        }
        if (tTopo->pxbLayer(id) == 2) {
          BRL2ADC.push_back((*iter).adc());
          BRL2Row.push_back((*iter).row());
          BRL2Col.push_back((*iter).column());
        }
        if (tTopo->pxbLayer(id) == 3) {
          BRL3ADC.push_back((*iter).adc());
          BRL3Row.push_back((*iter).row());
          BRL3Col.push_back((*iter).column());
        }
      }
    }

    // get Forward pixels
    if (detId.subdetId() == sdPxlFwd) {
      for (iter = begin; iter != end; ++iter) {
        ++nPxlFwd;
        if (tTopo->pxfDisk(id) == 1) {
          if (tTopo->pxfSide(id) == 1) {
            FWD1nADC.push_back((*iter).adc());
            FWD1nRow.push_back((*iter).row());
            FWD1nCol.push_back((*iter).column());
          }
          if (tTopo->pxfSide(id) == 2) {
            FWD1pADC.push_back((*iter).adc());
            FWD1pRow.push_back((*iter).row());
            FWD1pCol.push_back((*iter).column());
          }
        }
        if (tTopo->pxfDisk(id) == 2) {
          if (tTopo->pxfSide(id) == 1) {
            FWD2nADC.push_back((*iter).adc());
            FWD2nRow.push_back((*iter).row());
            FWD2nCol.push_back((*iter).column());
          }
          if (tTopo->pxfSide(id) == 2) {
            FWD2pADC.push_back((*iter).adc());
            FWD2pRow.push_back((*iter).row());
            FWD2pCol.push_back((*iter).column());
          }
        }
      }
    }
  }

  if (verbosity > 1) {
    eventout += "\n          Number of BrlPixelDigis collected:........ ";
    eventout += nPxlBrl;
  }

  if (verbosity > 1) {
    eventout += "\n          Number of FrwdPixelDigis collected:....... ";
    eventout += nPxlFwd;
  }

  if (verbosity > 0)
    edm::LogInfo(MsgLoggerCat) << eventout << "\n";

  return;
}

void GlobalDigisProducer::storeTrk(PGlobalDigi &product) {
  std::string MsgLoggerCat = "GlobalDigisProducer_storeTrk";

  if (verbosity > 2) {
    // strip output
    TString eventout("\n         nTIBL1     = ");
    eventout += TIBL1ADC.size();
    for (unsigned int i = 0; i < TIBL1ADC.size(); ++i) {
      eventout += "\n      (ADC, strip) = (";
      eventout += TIBL1ADC[i];
      eventout += ", ";
      eventout += TIBL1Strip[i];
      eventout += ")";
    }
    eventout += "\n         nTIBL2     = ";
    eventout += TIBL2ADC.size();
    for (unsigned int i = 0; i < TIBL2ADC.size(); ++i) {
      eventout += "\n      (ADC, strip) = (";
      eventout += TIBL2ADC[i];
      eventout += ", ";
      eventout += TIBL2Strip[i];
      eventout += ")";
    }
    eventout += "\n         nTIBL3     = ";
    eventout += TIBL3ADC.size();
    for (unsigned int i = 0; i < TIBL3ADC.size(); ++i) {
      eventout += "\n      (ADC, strip) = (";
      eventout += TIBL3ADC[i];
      eventout += ", ";
      eventout += TIBL3Strip[i];
      eventout += ")";
    }
    eventout += "\n         nTIBL4     = ";
    eventout += TIBL4ADC.size();
    for (unsigned int i = 0; i < TIBL4ADC.size(); ++i) {
      eventout += "\n      (ADC, strip) = (";
      eventout += TIBL4ADC[i];
      eventout += ", ";
      eventout += TIBL4Strip[i];
      eventout += ")";
    }
    eventout += "\n         nTOBL1     = ";
    eventout += TOBL1ADC.size();
    for (unsigned int i = 0; i < TOBL1ADC.size(); ++i) {
      eventout += "\n      (ADC, strip) = (";
      eventout += TOBL1ADC[i];
      eventout += ", ";
      eventout += TOBL1Strip[i];
      eventout += ")";
    }
    eventout += "\n         nTOBL2     = ";
    eventout += TOBL2ADC.size();
    for (unsigned int i = 0; i < TOBL2ADC.size(); ++i) {
      eventout += "\n      (ADC, strip) = (";
      eventout += TOBL2ADC[i];
      eventout += ", ";
      eventout += TOBL2Strip[i];
      eventout += ")";
    }
    eventout += "\n         nTOBL3     = ";
    eventout += TOBL3ADC.size();
    for (unsigned int i = 0; i < TOBL3ADC.size(); ++i) {
      eventout += "\n      (ADC, strip) = (";
      eventout += TOBL3ADC[i];
      eventout += ", ";
      eventout += TOBL3Strip[i];
      eventout += ")";
    }
    eventout += "\n         nTOBL4     = ";
    eventout += TOBL4ADC.size();
    for (unsigned int i = 0; i < TOBL4ADC.size(); ++i) {
      eventout += "\n      (ADC, strip) = (";
      eventout += TOBL4ADC[i];
      eventout += ", ";
      eventout += TOBL4Strip[i];
      eventout += ")";
    }
    eventout += "\n         nTIDW1     = ";
    eventout += TIDW1ADC.size();
    for (unsigned int i = 0; i < TIDW1ADC.size(); ++i) {
      eventout += "\n      (ADC, strip) = (";
      eventout += TIDW1ADC[i];
      eventout += ", ";
      eventout += TIDW1Strip[i];
      eventout += ")";
    }
    eventout += "\n         nTIDW2     = ";
    eventout += TIDW2ADC.size();
    for (unsigned int i = 0; i < TIDW2ADC.size(); ++i) {
      eventout += "\n      (ADC, strip) = (";
      eventout += TIDW2ADC[i];
      eventout += ", ";
      eventout += TIDW2Strip[i];
      eventout += ")";
    }
    eventout += "\n         nTIDW3     = ";
    eventout += TIDW3ADC.size();
    for (unsigned int i = 0; i < TIDW3ADC.size(); ++i) {
      eventout += "\n      (ADC, strip) = (";
      eventout += TIDW3ADC[i];
      eventout += ", ";
      eventout += TIDW3Strip[i];
      eventout += ")";
    }
    eventout += "\n         nTECW1     = ";
    eventout += TECW1ADC.size();
    for (unsigned int i = 0; i < TECW1ADC.size(); ++i) {
      eventout += "\n      (ADC, strip) = (";
      eventout += TECW1ADC[i];
      eventout += ", ";
      eventout += TECW1Strip[i];
      eventout += ")";
    }
    eventout += "\n         nTECW2     = ";
    eventout += TECW2ADC.size();
    for (unsigned int i = 0; i < TECW2ADC.size(); ++i) {
      eventout += "\n      (ADC, strip) = (";
      eventout += TECW2ADC[i];
      eventout += ", ";
      eventout += TECW2Strip[i];
      eventout += ")";
    }
    eventout += "\n         nTECW3     = ";
    eventout += TECW3ADC.size();
    for (unsigned int i = 0; i < TECW3ADC.size(); ++i) {
      eventout += "\n      (ADC, strip) = (";
      eventout += TECW3ADC[i];
      eventout += ", ";
      eventout += TECW3Strip[i];
      eventout += ")";
    }
    eventout += "\n         nTECW4     = ";
    eventout += TECW4ADC.size();
    for (unsigned int i = 0; i < TECW4ADC.size(); ++i) {
      eventout += "\n      (ADC, strip) = (";
      eventout += TECW4ADC[i];
      eventout += ", ";
      eventout += TECW4Strip[i];
      eventout += ")";
    }
    eventout += "\n         nTECW5     = ";
    eventout += TECW5ADC.size();
    for (unsigned int i = 0; i < TECW5ADC.size(); ++i) {
      eventout += "\n      (ADC, strip) = (";
      eventout += TECW5ADC[i];
      eventout += ", ";
      eventout += TECW5Strip[i];
      eventout += ")";
    }
    eventout += "\n         nTECW6     = ";
    eventout += TECW6ADC.size();
    for (unsigned int i = 0; i < TECW6ADC.size(); ++i) {
      eventout += "\n      (ADC, strip) = (";
      eventout += TECW6ADC[i];
      eventout += ", ";
      eventout += TECW6Strip[i];
      eventout += ")";
    }
    eventout += "\n         nTECW7     = ";
    eventout += TECW7ADC.size();
    for (unsigned int i = 0; i < TECW7ADC.size(); ++i) {
      eventout += "\n      (ADC, strip) = (";
      eventout += TECW7ADC[i];
      eventout += ", ";
      eventout += TECW7Strip[i];
      eventout += ")";
    }
    eventout += "\n         nTECW8     = ";
    eventout += TECW8ADC.size();
    for (unsigned int i = 0; i < TECW8ADC.size(); ++i) {
      eventout += "\n      (ADC, strip) = (";
      eventout += TECW8ADC[i];
      eventout += ", ";
      eventout += TECW8Strip[i];
      eventout += ")";
    }

    // pixel output
    eventout += "\n         nBRL1     = ";
    eventout += BRL1ADC.size();
    for (unsigned int i = 0; i < BRL1ADC.size(); ++i) {
      eventout += "\n      (ADC, row, column) = (";
      eventout += BRL1ADC[i];
      eventout += ", ";
      eventout += BRL1Row[i];
      eventout += ", ";
      eventout += BRL1Col[i];
      eventout += ")";
    }
    eventout += "\n         nBRL2     = ";
    eventout += BRL2ADC.size();
    for (unsigned int i = 0; i < BRL2ADC.size(); ++i) {
      eventout += "\n      (ADC, row, column) = (";
      eventout += BRL2ADC[i];
      eventout += ", ";
      eventout += BRL2Row[i];
      eventout += ", ";
      eventout += BRL2Col[i];
      eventout += ")";
    }
    eventout += "\n         nBRL3     = ";
    eventout += BRL3ADC.size();
    for (unsigned int i = 0; i < BRL3ADC.size(); ++i) {
      eventout += "\n      (ADC, row, column) = (";
      eventout += BRL3ADC[i];
      eventout += ", ";
      eventout += BRL3Row[i];
      eventout += ", ";
      eventout += BRL3Col[i];
      eventout += ")";
    }
    eventout += "\n         nFWD1p     = ";
    eventout += FWD1pADC.size();
    for (unsigned int i = 0; i < FWD1pADC.size(); ++i) {
      eventout += "\n      (ADC, row, column) = (";
      eventout += FWD1pADC[i];
      eventout += ", ";
      eventout += FWD1pRow[i];
      eventout += ", ";
      eventout += FWD1pCol[i];
      eventout += ")";
    }
    eventout += "\n         nFWD1p     = ";
    eventout += FWD1nADC.size();
    for (unsigned int i = 0; i < FWD1nADC.size(); ++i) {
      eventout += "\n      (ADC, row, column) = (";
      eventout += FWD1nADC[i];
      eventout += ", ";
      eventout += FWD1nRow[i];
      eventout += ", ";
      eventout += FWD1nCol[i];
      eventout += ")";
    }
    eventout += "\n         nFWD1p     = ";
    eventout += FWD2pADC.size();
    for (unsigned int i = 0; i < FWD2pADC.size(); ++i) {
      eventout += "\n      (ADC, row, column) = (";
      eventout += FWD2pADC[i];
      eventout += ", ";
      eventout += FWD2pRow[i];
      eventout += ", ";
      eventout += FWD2pCol[i];
      eventout += ")";
    }
    eventout += "\n         nFWD2p     = ";
    eventout += FWD2nADC.size();
    for (unsigned int i = 0; i < FWD2nADC.size(); ++i) {
      eventout += "\n      (ADC, row, column) = (";
      eventout += FWD2nADC[i];
      eventout += ", ";
      eventout += FWD2nRow[i];
      eventout += ", ";
      eventout += FWD2nCol[i];
      eventout += ")";
    }

    edm::LogInfo(MsgLoggerCat) << eventout << "\n";
  }

  // strip output
  product.putTIBL1Digis(TIBL1ADC, TIBL1Strip);
  product.putTIBL2Digis(TIBL2ADC, TIBL2Strip);
  product.putTIBL3Digis(TIBL3ADC, TIBL3Strip);
  product.putTIBL4Digis(TIBL4ADC, TIBL4Strip);
  product.putTOBL1Digis(TOBL1ADC, TOBL1Strip);
  product.putTOBL2Digis(TOBL2ADC, TOBL2Strip);
  product.putTOBL3Digis(TOBL3ADC, TOBL3Strip);
  product.putTOBL4Digis(TOBL4ADC, TOBL4Strip);
  product.putTIDW1Digis(TIDW1ADC, TIDW1Strip);
  product.putTIDW2Digis(TIDW2ADC, TIDW2Strip);
  product.putTIDW3Digis(TIDW3ADC, TIDW3Strip);
  product.putTECW1Digis(TECW1ADC, TECW1Strip);
  product.putTECW2Digis(TECW2ADC, TECW2Strip);
  product.putTECW3Digis(TECW3ADC, TECW3Strip);
  product.putTECW4Digis(TECW4ADC, TECW4Strip);
  product.putTECW5Digis(TECW5ADC, TECW5Strip);
  product.putTECW6Digis(TECW6ADC, TECW6Strip);
  product.putTECW7Digis(TECW7ADC, TECW7Strip);
  product.putTECW8Digis(TECW8ADC, TECW8Strip);

  // pixel output
  product.putBRL1Digis(BRL1ADC, BRL1Row, BRL1Col);
  product.putBRL2Digis(BRL2ADC, BRL2Row, BRL2Col);
  product.putBRL3Digis(BRL3ADC, BRL3Row, BRL3Col);
  product.putFWD1pDigis(FWD1pADC, FWD1pRow, FWD1pCol);
  product.putFWD1nDigis(FWD1nADC, FWD1nRow, FWD1nCol);
  product.putFWD2pDigis(FWD2pADC, FWD2pRow, FWD2pCol);
  product.putFWD2nDigis(FWD2nADC, FWD2nRow, FWD2nCol);

  return;
}

void GlobalDigisProducer::fillMuon(edm::Event &iEvent, const edm::EventSetup &iSetup) {
  std::string MsgLoggerCat = "GlobalDigisProducer_fillMuon";

  TString eventout;
  if (verbosity > 0)
    eventout = "\nGathering info:";

  // get DT information
  edm::Handle<DTDigiCollection> dtDigis;
  iEvent.getByToken(MuDTSrc_Token_, dtDigis);
  if (!dtDigis.isValid()) {
    edm::LogWarning(MsgLoggerCat) << "Unable to find dtDigis in event!";
    return;
  }

  int nDt = 0;
  DTDigiCollection::DigiRangeIterator detUnitIt;
  for (detUnitIt = dtDigis->begin(); detUnitIt != dtDigis->end(); ++detUnitIt) {
    const DTLayerId &id = (*detUnitIt).first;
    const DTDigiCollection::Range &range = (*detUnitIt).second;

    for (DTDigiCollection::const_iterator digiIt = range.first; digiIt != range.second; ++digiIt) {
      ++nDt;

      DTWireId wireId(id, (*digiIt).wire());
      if (wireId.station() == 1) {
        MB1SLayer.push_back(id.superlayer());
        MB1Time.push_back((*digiIt).time());
        MB1Layer.push_back(id.layer());
      }
      if (wireId.station() == 2) {
        MB2SLayer.push_back(id.superlayer());
        MB2Time.push_back((*digiIt).time());
        MB2Layer.push_back(id.layer());
      }
      if (wireId.station() == 3) {
        MB3SLayer.push_back(id.superlayer());
        MB3Time.push_back((*digiIt).time());
        MB3Layer.push_back(id.layer());
      }
      if (wireId.station() == 4) {
        MB4SLayer.push_back(id.superlayer());
        MB4Time.push_back((*digiIt).time());
        MB4Layer.push_back(id.layer());
      }
    }
  }

  if (verbosity > 1) {
    eventout += "\n          Number of DtMuonDigis collected:.......... ";
    eventout += nDt;
  }

  // get CSC Strip information
  edm::Handle<CSCStripDigiCollection> strips;
  iEvent.getByToken(MuCSCStripSrc_Token_, strips);
  if (!strips.isValid()) {
    edm::LogWarning(MsgLoggerCat) << "Unable to find muon strips in event!";
    return;
  }

  int nStrips = 0;
  for (CSCStripDigiCollection::DigiRangeIterator j = strips->begin(); j != strips->end(); ++j) {
    std::vector<CSCStripDigi>::const_iterator digiItr = (*j).second.first;
    std::vector<CSCStripDigi>::const_iterator last = (*j).second.second;

    for (; digiItr != last; ++digiItr) {
      ++nStrips;

      // average pedestals
      std::vector<int> adcCounts = digiItr->getADCCounts();
      theCSCStripPedestalSum += adcCounts[0];
      theCSCStripPedestalSum += adcCounts[1];
      theCSCStripPedestalCount += 2;

      // if there are enough pedestal statistics
      if (theCSCStripPedestalCount > 100) {
        float pedestal = theCSCStripPedestalSum / theCSCStripPedestalCount;
        if (adcCounts[5] > (pedestal + 100))
          CSCStripADC.push_back(adcCounts[4] - pedestal);
      }
    }
  }

  if (verbosity > 1) {
    eventout += "\n          Number of CSCStripDigis collected:........ ";
    eventout += nStrips;
  }

  // get CSC Wire information
  edm::Handle<CSCWireDigiCollection> wires;
  iEvent.getByToken(MuCSCWireSrc_Token_, wires);
  if (!wires.isValid()) {
    edm::LogWarning(MsgLoggerCat) << "Unable to find muon wires in event!";
    return;
  }

  int nWires = 0;
  for (CSCWireDigiCollection::DigiRangeIterator j = wires->begin(); j != wires->end(); ++j) {
    std::vector<CSCWireDigi>::const_iterator digiItr = (*j).second.first;
    std::vector<CSCWireDigi>::const_iterator endDigi = (*j).second.second;

    for (; digiItr != endDigi; ++digiItr) {
      ++nWires;

      CSCWireTime.push_back(digiItr->getTimeBin());
    }
  }

  if (verbosity > 1) {
    eventout += "\n          Number of CSCWireDigis collected:......... ";
    eventout += nWires;
  }

  if (verbosity > 0)
    edm::LogInfo(MsgLoggerCat) << eventout << "\n";

  return;
}

void GlobalDigisProducer::storeMuon(PGlobalDigi &product) {
  std::string MsgLoggerCat = "GlobalDigisProducer_storeMuon";

  if (verbosity > 2) {
    // dt output
    TString eventout("\n         nMB1     = ");
    eventout += MB1SLayer.size();
    for (unsigned int i = 0; i < MB1SLayer.size(); ++i) {
      eventout += "\n      (slayer, time, layer) = (";
      eventout += MB1SLayer[i];
      eventout += ", ";
      eventout += MB1Time[i];
      eventout += ", ";
      eventout += MB1Layer[i];
      eventout += ")";
    }
    eventout += "\n         nMB2     = ";
    eventout += MB2SLayer.size();
    for (unsigned int i = 0; i < MB2SLayer.size(); ++i) {
      eventout += "\n      (slayer, time, layer) = (";
      eventout += MB2SLayer[i];
      eventout += ", ";
      eventout += MB2Time[i];
      eventout += ", ";
      eventout += MB2Layer[i];
      eventout += ")";
    }
    eventout += "\n         nMB3     = ";
    eventout += MB3SLayer.size();
    for (unsigned int i = 0; i < MB3SLayer.size(); ++i) {
      eventout += "\n      (slayer, time, layer) = (";
      eventout += MB3SLayer[i];
      eventout += ", ";
      eventout += MB3Time[i];
      eventout += ", ";
      eventout += MB3Layer[i];
      eventout += ")";
    }
    eventout += "\n         nMB2     = ";
    eventout += MB4SLayer.size();
    for (unsigned int i = 0; i < MB4SLayer.size(); ++i) {
      eventout += "\n      (slayer, time, layer) = (";
      eventout += MB4SLayer[i];
      eventout += ", ";
      eventout += MB4Time[i];
      eventout += ", ";
      eventout += MB4Layer[i];
      eventout += ")";
    }

    // CSC Strip
    eventout += "\n         nCSCStrip     = ";
    eventout += CSCStripADC.size();
    for (unsigned int i = 0; i < CSCStripADC.size(); ++i) {
      eventout += "\n      (adc) = (";
      eventout += CSCStripADC[i];
      eventout += ")";
    }

    // CSC Wire
    eventout += "\n         nCSCWire     = ";
    eventout += CSCWireTime.size();
    for (unsigned int i = 0; i < CSCWireTime.size(); ++i) {
      eventout += "\n      (time) = (";
      eventout += CSCWireTime[i];
      eventout += ")";
    }

    edm::LogInfo(MsgLoggerCat) << eventout << "\n";
  }

  product.putMB1Digis(MB1SLayer, MB1Time, MB1Layer);
  product.putMB2Digis(MB2SLayer, MB2Time, MB2Layer);
  product.putMB3Digis(MB3SLayer, MB3Time, MB3Layer);
  product.putMB4Digis(MB4SLayer, MB4Time, MB4Layer);

  product.putCSCstripDigis(CSCStripADC);

  product.putCSCwireDigis(CSCWireTime);

  return;
}

void GlobalDigisProducer::clear() {
  std::string MsgLoggerCat = "GlobalDigisProducer_clear";

  if (verbosity > 0)
    edm::LogInfo(MsgLoggerCat) << "Clearing event holders";

  // reset electromagnetic info
  // EE info
  EECalmaxPos.clear();
  EECalAEE.clear();
  EECalSHE.clear();
  // EB info
  EBCalmaxPos.clear();
  EBCalAEE.clear();
  EBCalSHE.clear();
  // ES info
  ESCalADC0.clear();
  ESCalADC1.clear();
  ESCalADC2.clear();
  ESCalSHE.clear();

  // reset HCal Info
  HBCalAEE.clear();
  HBCalSHE.clear();
  HECalAEE.clear();
  HECalSHE.clear();
  HOCalAEE.clear();
  HOCalSHE.clear();
  HFCalAEE.clear();
  HFCalSHE.clear();

  // reset Track Info
  TIBL1ADC.clear();
  TIBL2ADC.clear();
  TIBL3ADC.clear();
  TIBL4ADC.clear();
  TIBL1Strip.clear();
  TIBL2Strip.clear();
  TIBL3Strip.clear();
  TIBL4Strip.clear();

  TOBL1ADC.clear();
  TOBL2ADC.clear();
  TOBL3ADC.clear();
  TOBL4ADC.clear();
  TOBL1Strip.clear();
  TOBL2Strip.clear();
  TOBL3Strip.clear();
  TOBL4Strip.clear();

  TIDW1ADC.clear();
  TIDW2ADC.clear();
  TIDW3ADC.clear();
  TIDW1Strip.clear();
  TIDW2Strip.clear();
  TIDW3Strip.clear();

  TECW1ADC.clear();
  TECW2ADC.clear();
  TECW3ADC.clear();
  TECW4ADC.clear();
  TECW5ADC.clear();
  TECW6ADC.clear();
  TECW7ADC.clear();
  TECW8ADC.clear();
  TECW1Strip.clear();
  TECW2Strip.clear();
  TECW3Strip.clear();
  TECW4Strip.clear();
  TECW5Strip.clear();
  TECW6Strip.clear();
  TECW7Strip.clear();
  TECW8Strip.clear();

  BRL1ADC.clear();
  BRL1Row.clear();
  BRL1Col.clear();
  BRL2ADC.clear();
  BRL2Row.clear();
  BRL2Col.clear();
  BRL3ADC.clear();
  BRL3Row.clear();
  BRL3Col.clear();

  FWD1pADC.clear();
  FWD1pRow.clear();
  FWD1pCol.clear();
  FWD1nADC.clear();
  FWD1nRow.clear();
  FWD1nCol.clear();
  FWD2pADC.clear();
  FWD2pRow.clear();
  FWD2pCol.clear();
  FWD2nADC.clear();
  FWD2nRow.clear();
  FWD2nCol.clear();

  // muon clear
  MB1SLayer.clear();
  MB1Time.clear();
  MB1Layer.clear();
  MB2SLayer.clear();
  MB2Time.clear();
  MB2Layer.clear();
  MB3SLayer.clear();
  MB3Time.clear();
  MB3Layer.clear();
  MB4SLayer.clear();
  MB4Time.clear();
  MB4Layer.clear();

  CSCStripADC.clear();

  CSCWireTime.clear();

  return;
}

// define this as a plug-in
// DEFINE_FWK_MODULE(GlobalDigisProducer);