HcalHitReconstructor.cc

Go to the documentation of this file.

#include "HcalHitReconstructor.h"
#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
#include "DataFormats/Common/interface/EDCollection.h"
#include "DataFormats/Common/interface/Handle.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "CalibFormats/HcalObjects/interface/HcalCoderDb.h"
#include "CalibFormats/HcalObjects/interface/HcalDbRecord.h"
#include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputer.h"
#include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputerRcd.h"
#include "CalibCalorimetry/HcalAlgos/interface/HcalDbASCIIIO.h"
#include "Geometry/CaloTopology/interface/HcalTopology.h"
#include "Geometry/Records/interface/HcalRecNumberingRecord.h"
#include "CondFormats/DataRecord/interface/HcalFrontEndMapRcd.h"
#include "CondFormats/DataRecord/interface/HcalOOTPileupCorrectionRcd.h"
#include "CondFormats/DataRecord/interface/HcalOOTPileupCompatibilityRcd.h"
#include "CondFormats/DataRecord/interface/HBHENegativeEFilterRcd.h"
#include "CondFormats/HcalObjects/interface/HcalFrontEndMap.h"
#include "CondFormats/HcalObjects/interface/OOTPileupCorrectionColl.h"
#include "CondFormats/HcalObjects/interface/OOTPileupCorrData.h"
#include <iostream>
#include <fstream>


/*  Hcal Hit reconstructor allows for CaloRecHits with status words */

HcalHitReconstructor::HcalHitReconstructor(edm::ParameterSet const& conf):
  reco_(conf.getParameter<bool>("correctForTimeslew"),
    conf.getParameter<bool>("correctForPhaseContainment"),
    conf.getParameter<double>("correctionPhaseNS")),
  det_(DetId::Hcal),
  inputLabel_(conf.getParameter<edm::InputTag>("digiLabel")),
  correctTiming_(conf.getParameter<bool>("correctTiming")),
  setNoiseFlags_(conf.getParameter<bool>("setNoiseFlags")),
  setHSCPFlags_(conf.getParameter<bool>("setHSCPFlags")),
  setSaturationFlags_(conf.getParameter<bool>("setSaturationFlags")),
  setTimingTrustFlags_(conf.getParameter<bool>("setTimingTrustFlags")),
  setPulseShapeFlags_(conf.getParameter<bool>("setPulseShapeFlags")),
  setNegativeFlags_(false),
  dropZSmarkedPassed_(conf.getParameter<bool>("dropZSmarkedPassed")),
  firstAuxTS_(conf.getParameter<int>("firstAuxTS")),
  firstSample_(conf.getParameter<int>("firstSample")),
  samplesToAdd_(conf.getParameter<int>("samplesToAdd")),
  tsFromDB_(conf.getParameter<bool>("tsFromDB")),
  useLeakCorrection_(conf.getParameter<bool>("useLeakCorrection")),
  dataOOTCorrectionName_(""),
  dataOOTCorrectionCategory_("Data"),
  mcOOTCorrectionName_(""),
  mcOOTCorrectionCategory_("MC"),
  setPileupCorrection_(nullptr),
  paramTS(nullptr)

{
  // register for data access
  tok_ho_ = consumes<HODigiCollection>(inputLabel_);
  tok_hf_ = consumes<HFDigiCollection>(inputLabel_);
  tok_calib_ = consumes<HcalCalibDigiCollection>(inputLabel_);

  std::string subd=conf.getParameter<std::string>("Subdetector");
  //Set all FlagSetters to 0
  /* Important to do this!  Otherwise, if the setters are turned off,
     the "if (XSetter_) delete XSetter_;" commands can crash
  */

  recoParamsFromDB_ = conf.getParameter<bool>("recoParamsFromDB");
  //  recoParamsFromDB_ = false ; //  trun off for now.

  // std::cout<<"  HcalHitReconstructor   recoParamsFromDB_ "<<recoParamsFromDB_<<std::endl;

  if (conf.existsAs<bool>("setNegativeFlags"))
      setNegativeFlags_ = conf.getParameter<bool>("setNegativeFlags");

  hfdigibit_                  = nullptr;

  hfS9S1_                     = nullptr;
  hfS8S1_                     = nullptr;
  hfPET_                      = nullptr;
  saturationFlagSetter_       = nullptr;
  HFTimingTrustFlagSetter_    = nullptr;
  digiTimeFromDB_             = false; // only need for HF
  
  if (setSaturationFlags_)
    {
      const edm::ParameterSet& pssat      = conf.getParameter<edm::ParameterSet>("saturationParameters");
      saturationFlagSetter_ = new HcalADCSaturationFlag(pssat.getParameter<int>("maxADCvalue"));
    }


  if (!strcasecmp(subd.c_str(),"HO")) {
    subdet_=HcalOuter;
    // setPileupCorrection_ = &HcalSimpleRecAlgo::setHOPileupCorrection;
    setPileupCorrection_ = nullptr;
    produces<HORecHitCollection>();
  } else if (!strcasecmp(subd.c_str(),"HF")) {
    subdet_=HcalForward;
    // setPileupCorrection_ = &HcalSimpleRecAlgo::setHFPileupCorrection;
    setPileupCorrection_ = nullptr;
    digiTimeFromDB_=conf.getParameter<bool>("digiTimeFromDB");

    if (setTimingTrustFlags_) {
      
      const edm::ParameterSet& pstrust      = conf.getParameter<edm::ParameterSet>("hfTimingTrustParameters");
      HFTimingTrustFlagSetter_=new HFTimingTrustFlag(pstrust.getParameter<int>("hfTimingTrustLevel1"),
                             pstrust.getParameter<int>("hfTimingTrustLevel2"));
    }

    if (setNoiseFlags_)
      {
    const edm::ParameterSet& psdigi    =conf.getParameter<edm::ParameterSet>("digistat");
    const edm::ParameterSet& psTimeWin =conf.getParameter<edm::ParameterSet>("HFInWindowStat");
    hfdigibit_=new HcalHFStatusBitFromDigis(psdigi,psTimeWin);

    const edm::ParameterSet& psS9S1   = conf.getParameter<edm::ParameterSet>("S9S1stat");
    hfS9S1_   = new HcalHF_S9S1algorithm(psS9S1.getParameter<std::vector<double> >("short_optimumSlope"),
                         psS9S1.getParameter<std::vector<double> >("shortEnergyParams"),
                         psS9S1.getParameter<std::vector<double> >("shortETParams"),
                         psS9S1.getParameter<std::vector<double> >("long_optimumSlope"),
                         psS9S1.getParameter<std::vector<double> >("longEnergyParams"),
                         psS9S1.getParameter<std::vector<double> >("longETParams"),
                         psS9S1.getParameter<int>("HcalAcceptSeverityLevel"),
                         psS9S1.getParameter<bool>("isS8S1")
                         );

    const edm::ParameterSet& psS8S1   = conf.getParameter<edm::ParameterSet>("S8S1stat");
    hfS8S1_   = new HcalHF_S9S1algorithm(psS8S1.getParameter<std::vector<double> >("short_optimumSlope"),
                         psS8S1.getParameter<std::vector<double> >("shortEnergyParams"),
                         psS8S1.getParameter<std::vector<double> >("shortETParams"),
                         psS8S1.getParameter<std::vector<double> >("long_optimumSlope"),
                         psS8S1.getParameter<std::vector<double> >("longEnergyParams"),
                         psS8S1.getParameter<std::vector<double> >("longETParams"),
                         psS8S1.getParameter<int>("HcalAcceptSeverityLevel"),
                         psS8S1.getParameter<bool>("isS8S1")
                         );

    const edm::ParameterSet& psPET    = conf.getParameter<edm::ParameterSet>("PETstat");
    hfPET_    = new HcalHF_PETalgorithm(psPET.getParameter<std::vector<double> >("short_R"),
                        psPET.getParameter<std::vector<double> >("shortEnergyParams"),
                        psPET.getParameter<std::vector<double> >("shortETParams"),
                        psPET.getParameter<std::vector<double> >("long_R"),
                        psPET.getParameter<std::vector<double> >("longEnergyParams"),
                        psPET.getParameter<std::vector<double> >("longETParams"),
                        psPET.getParameter<int>("HcalAcceptSeverityLevel"),
                        psPET.getParameter<std::vector<double> >("short_R_29"),
                        psPET.getParameter<std::vector<double> >("long_R_29")
                        );
      }
    produces<HFRecHitCollection>();
  } else if (!strcasecmp(subd.c_str(),"ZDC")) {
    det_=DetId::Calo;
    subdet_=HcalZDCDetId::SubdetectorId;
    produces<ZDCRecHitCollection>();
  } else if (!strcasecmp(subd.c_str(),"CALIB")) {
    subdet_=HcalOther;
    subdetOther_=HcalCalibration;
    produces<HcalCalibRecHitCollection>();
  } else {
    edm::LogWarning("Configuration") << "HcalHitReconstructor is not associated with a specific subdetector!" << std::endl;
  }

  // If no valid OOT pileup correction name specified,
  // disable the correction
  if (conf.existsAs<std::string>("dataOOTCorrectionName"))
      dataOOTCorrectionName_ = conf.getParameter<std::string>("dataOOTCorrectionName");
  if (conf.existsAs<std::string>("dataOOTCorrectionCategory"))
      dataOOTCorrectionCategory_ = conf.getParameter<std::string>("dataOOTCorrectionCategory");
  if (conf.existsAs<std::string>("mcOOTCorrectionName"))
      mcOOTCorrectionName_ = conf.getParameter<std::string>("mcOOTCorrectionName");
  if (conf.existsAs<std::string>("mcOOTCorrectionCategory"))
      mcOOTCorrectionCategory_ = conf.getParameter<std::string>("mcOOTCorrectionCategory");
  if (dataOOTCorrectionName_.empty() && mcOOTCorrectionName_.empty())
      setPileupCorrection_ = nullptr;

}

HcalHitReconstructor::~HcalHitReconstructor() {
  delete hfdigibit_;
  
  delete hfS9S1_;
  delete hfS8S1_;
  delete hfPET_;
  delete saturationFlagSetter_;
  delete HFTimingTrustFlagSetter_;

  delete paramTS;
}

void HcalHitReconstructor::beginRun(edm::Run const&r, edm::EventSetup const & es){

  edm::ESHandle<HcalTopology> htopo;
  es.get<HcalRecNumberingRecord>().get(htopo);

  if ( tsFromDB_== true || recoParamsFromDB_ == true )
    {
      edm::ESHandle<HcalRecoParams> p;
      es.get<HcalRecoParamsRcd>().get(p);
      paramTS = new HcalRecoParams(*p.product());
      paramTS->setTopo(htopo.product());

      


      // std::cout<<" skdump in HcalHitReconstructor::beginRun   dupm RecoParams "<<std::endl;
      // std::ofstream skfile("skdumpRecoParamsNewFormat.txt");
      // HcalDbASCIIIO::dumpObject(skfile, (*paramTS) );
    }

  if (digiTimeFromDB_==true)
    {
      edm::ESHandle<HcalFlagHFDigiTimeParams> p;
      es.get<HcalFlagHFDigiTimeParamsRcd>().get(p);
      HFDigiTimeParams.reset( new HcalFlagHFDigiTimeParams( *p ) );
      HFDigiTimeParams->setTopo(htopo.product());
    }

  reco_.beginRun(es);
}

void HcalHitReconstructor::endRun(edm::Run const&r, edm::EventSetup const & es){
  if (tsFromDB_==true)
    {
      delete paramTS; paramTS=nullptr;
    }
  if (digiTimeFromDB_==true)
    {
      //DL delete HFDigiTimeParams; HFDigiTimeParams = 0;
    }
  reco_.endRun();
}

void HcalHitReconstructor::produce(edm::Event& e, const edm::EventSetup& eventSetup)
{

  // get conditions
  edm::ESHandle<HcalTopology> topo;
  eventSetup.get<HcalRecNumberingRecord>().get(topo);

  edm::ESHandle<HcalDbService> conditions;
  eventSetup.get<HcalDbRecord>().get(conditions);

  if (useLeakCorrection_) reco_.setLeakCorrection();

  edm::ESHandle<HcalChannelQuality> p;
  eventSetup.get<HcalChannelQualityRcd>().get("withTopo",p);
  const HcalChannelQuality* myqual = p.product();

  edm::ESHandle<HcalSeverityLevelComputer> mycomputer;
  eventSetup.get<HcalSeverityLevelComputerRcd>().get(mycomputer);
  const HcalSeverityLevelComputer* mySeverity = mycomputer.product();


  // GET THE BEAM CROSSING INFO HERE, WHEN WE UNDERSTAND HOW THINGS WORK.
  // Then, call "setBXInfo" method of the reco_ object.
  // Also remember to call SetBXInfo in the negative energy flag setter.

  if (det_==DetId::Hcal) {

      //  HO ------------------------------------------------------------------
    if (subdet_==HcalOuter) {
      edm::Handle<HODigiCollection> digi;
      e.getByToken(tok_ho_,digi);
      
      // create empty output
      auto rec = std::make_unique<HORecHitCollection>();
      rec->reserve(digi->size());
      // run the algorithm
      HODigiCollection::const_iterator i;

      // Vote on majority TS0 CapId
      int favorite_capid = 0; 
      if (correctTiming_) {
        long capid_votes[4] = {0,0,0,0};
        for (i=digi->begin(); i!=digi->end(); i++) {
          capid_votes[(*i)[0].capid()]++;
        }
        for (int k = 0; k < 4; k++)
          if (capid_votes[k] > capid_votes[favorite_capid])
            favorite_capid = k;
      }

      for (i=digi->begin(); i!=digi->end(); i++) {
    HcalDetId cell = i->id();
    DetId detcell=(DetId)cell;
        // firstSample & samplesToAdd
        if(tsFromDB_ || recoParamsFromDB_) {
          const HcalRecoParam* param_ts = paramTS->getValues(detcell.rawId());
      if(tsFromDB_) {
        firstSample_  = param_ts->firstSample();
        samplesToAdd_ = param_ts->samplesToAdd();
      }
          if(recoParamsFromDB_) {
             bool correctForTimeslew=param_ts->correctForTimeslew();
             bool correctForPhaseContainment= param_ts->correctForPhaseContainment();
             float phaseNS=param_ts->correctionPhaseNS();
             useLeakCorrection_= param_ts->useLeakCorrection();
             correctTiming_ = param_ts->correctTiming();
             firstAuxTS_ = param_ts->firstAuxTS();
             int pileupCleaningID = param_ts->pileupCleaningID();
             reco_.setRecoParams(correctForTimeslew,correctForPhaseContainment,useLeakCorrection_,pileupCleaningID,phaseNS);
          }
        }

        int first = firstSample_;
        int toadd = samplesToAdd_;

    // check on cells to be ignored and dropped: (rof,20.Feb.09)
    const HcalChannelStatus* mydigistatus=myqual->getValues(detcell.rawId());
    if (mySeverity->dropChannel(mydigistatus->getValue() ) ) continue;
    if (dropZSmarkedPassed_)
      if (i->zsMarkAndPass()) continue;

    const HcalCalibrations& calibrations=conditions->getHcalCalibrations(cell);
    const HcalQIECoder* channelCoder = conditions->getHcalCoder (cell);
    const HcalQIEShape* shape = conditions->getHcalShape (channelCoder);
    HcalCoderDb coder (*channelCoder, *shape);

    rec->push_back(reco_.reconstruct(*i,first,toadd,coder,calibrations));

    // Set auxiliary flag
    int auxflag=0;
        int fTS = firstAuxTS_;
    if (fTS<0) fTS=0; //silly protection against negative time slice values
    for (int xx=fTS; xx<fTS+4 && xx<i->size();++xx)
      auxflag+=(i->sample(xx).adc())<<(7*(xx-fTS)); // store the time slices in the first 28 bits of aux, a set of 4 7-bit adc values
    // bits 28 and 29 are reserved for capid of the first time slice saved in aux
    auxflag+=((i->sample(fTS).capid())<<28);
    (rec->back()).setAux(auxflag);
    // (rec->back()).setFlags(0);  Don't want to do this because the algorithm
        //                             can already set some flags
    // Fill Presample ADC flag
    if (fTS>0)
      (rec->back()).setFlagField((i->sample(fTS-1).adc()), HcalCaloFlagLabels::PresampleADC,7);

    if (setSaturationFlags_)
      saturationFlagSetter_->setSaturationFlag(rec->back(),*i);
    if (correctTiming_)
      HcalTimingCorrector::Correct(rec->back(), *i, favorite_capid);
      }
      // return result
      e.put(std::move(rec));    

      // HF -------------------------------------------------------------------
    } else if (subdet_==HcalForward) {
      edm::Handle<HFDigiCollection> digi;
      e.getByToken(tok_hf_,digi);


      // create empty output
      auto rec = std::make_unique<HFRecHitCollection>();
      rec->reserve(digi->size());
      // run the algorithm
      HFDigiCollection::const_iterator i;

      // Vote on majority TS0 CapId
      int favorite_capid = 0; 
      if (correctTiming_) {
        long capid_votes[4] = {0,0,0,0};
        for (i=digi->begin(); i!=digi->end(); i++) {
          capid_votes[(*i)[0].capid()]++;
        }
        for (int k = 0; k < 4; k++)
          if (capid_votes[k] > capid_votes[favorite_capid])
            favorite_capid = k;
      }

      for (i=digi->begin(); i!=digi->end(); i++) {
    HcalDetId cell = i->id();
    DetId detcell=(DetId)cell;

        if(tsFromDB_ || recoParamsFromDB_) {
          const HcalRecoParam* param_ts = paramTS->getValues(detcell.rawId());
      if(tsFromDB_) {
        firstSample_  = param_ts->firstSample();
        samplesToAdd_ = param_ts->samplesToAdd();
      }
          if(recoParamsFromDB_) {
             bool correctForTimeslew=param_ts->correctForTimeslew();
             bool correctForPhaseContainment= param_ts->correctForPhaseContainment();
             float phaseNS=param_ts->correctionPhaseNS();
             useLeakCorrection_= param_ts->useLeakCorrection();
             correctTiming_ = param_ts->correctTiming();
             firstAuxTS_ = param_ts->firstAuxTS();
             int pileupCleaningID = param_ts->pileupCleaningID();
             reco_.setRecoParams(correctForTimeslew,correctForPhaseContainment,useLeakCorrection_,pileupCleaningID,phaseNS);
          }
        }

        int first = firstSample_;
        int toadd = samplesToAdd_;

    // check on cells to be ignored and dropped: (rof,20.Feb.09)
    const HcalChannelStatus* mydigistatus=myqual->getValues(detcell.rawId());
    if (mySeverity->dropChannel(mydigistatus->getValue() ) ) continue;
    if (dropZSmarkedPassed_)
      if (i->zsMarkAndPass()) continue;

    const HcalCalibrations& calibrations=conditions->getHcalCalibrations(cell);
    const HcalQIECoder* channelCoder = conditions->getHcalCoder (cell);
    const HcalQIEShape* shape = conditions->getHcalShape (channelCoder);
    HcalCoderDb coder (*channelCoder, *shape);

    // Set HFDigiTime flag values from digiTimeFromDB_
    if (digiTimeFromDB_==true && hfdigibit_!=nullptr)
      {
        const HcalFlagHFDigiTimeParam* hfDTparam = HFDigiTimeParams->getValues(detcell.rawId());
        hfdigibit_->resetParamsFromDB(hfDTparam->HFdigiflagFirstSample(),
                      hfDTparam->HFdigiflagSamplesToAdd(),
                      hfDTparam->HFdigiflagExpectedPeak(),
                      hfDTparam->HFdigiflagMinEThreshold(),
                      hfDTparam->HFdigiflagCoefficients()
                      );
      }

    //std::cout << "TOADDHF " << toadd << " " << first << " " << std::endl;
    rec->push_back(reco_.reconstruct(*i,first,toadd,coder,calibrations));

    // Set auxiliary flag
    int auxflag=0;
        int fTS = firstAuxTS_;
    if (fTS<0) fTS=0; // silly protection against negative time slice values
    for (int xx=fTS; xx<fTS+4 && xx<i->size();++xx)
      auxflag+=(i->sample(xx).adc())<<(7*(xx-fTS)); // store the time slices in the first 28 bits of aux, a set of 4 7-bit adc values
    // bits 28 and 29 are reserved for capid of the first time slice saved in aux
    auxflag+=((i->sample(fTS).capid())<<28);
    (rec->back()).setAux(auxflag);

    // (rec->back()).setFlags(0);  Don't want to do this because the algorithm
        //                             can already set some flags

    // Fill Presample ADC flag
    if (fTS>0)
      (rec->back()).setFlagField((i->sample(fTS-1).adc()), HcalCaloFlagLabels::PresampleADC,7);

    // This calls the code for setting the HF noise bit determined from digi shape
    if (setNoiseFlags_) 
      hfdigibit_->hfSetFlagFromDigi(rec->back(),*i,coder,calibrations);
    if (setSaturationFlags_)
      saturationFlagSetter_->setSaturationFlag(rec->back(),*i);
    if (setTimingTrustFlags_)
      HFTimingTrustFlagSetter_->setHFTimingTrustFlag(rec->back(),*i);
    if (correctTiming_)
      HcalTimingCorrector::Correct(rec->back(), *i, favorite_capid);
      } // for (i=digi->begin(); i!=digi->end(); i++) -- loop on all HF digis

      // The following flags require the full set of rechits
      // These need to be set consecutively, so an energy check should be the first 
      // test performed on these hits (to minimize the loop time)
      if (setNoiseFlags_) 
    {
      // Step 1:  Set PET flag  (short fibers of |ieta|==29)
      // Neighbor/partner channels that are flagged by Pulse Shape algorithm (HFDigiTime)
      // won't be considered in these calculations
      for (HFRecHitCollection::iterator i = rec->begin();i!=rec->end();++i)
        {
          int depth=i->id().depth();
          int ieta=i->id().ieta();
          // Short fibers and all channels at |ieta|=29 use PET settings in Algo 3
          if (depth==2 || abs(ieta)==29 ) 
        hfPET_->HFSetFlagFromPET(*i,*rec,myqual,mySeverity);
        }

      // Step 2:  Set S8S1 flag (short fibers or |ieta|==29)
      for (HFRecHitCollection::iterator i = rec->begin();i!=rec->end();++i)
        {
          int depth=i->id().depth();
          int ieta=i->id().ieta();
          // Short fibers and all channels at |ieta|=29 use PET settings in Algo 3
          if (depth==2 || abs(ieta)==29 ) 
        hfS8S1_->HFSetFlagFromS9S1(*i,*rec,myqual,mySeverity);
        }

      // Set 3:  Set S9S1 flag (long fibers)
      for (HFRecHitCollection::iterator i = rec->begin();i!=rec->end();++i)
        {
          int depth=i->id().depth();
          int ieta=i->id().ieta();
          // Short fibers and all channels at |ieta|=29 use PET settings in Algo 3
          if (depth==1 && abs(ieta)!=29 ) 
        hfS9S1_->HFSetFlagFromS9S1(*i,*rec,myqual, mySeverity);
        }
    }

      // return result
      e.put(std::move(rec));     
    } else if (subdet_==HcalOther && subdetOther_==HcalCalibration) {
      edm::Handle<HcalCalibDigiCollection> digi;
      e.getByToken(tok_calib_,digi);
      
      // create empty output
      auto rec = std::make_unique<HcalCalibRecHitCollection>();
      rec->reserve(digi->size());
      // run the algorithm
      int first = firstSample_;
      int toadd = samplesToAdd_;

      HcalCalibDigiCollection::const_iterator i;
      for (i=digi->begin(); i!=digi->end(); i++) {
    HcalCalibDetId cell = i->id();
    //  HcalDetId cellh = i->id();
    DetId detcell=(DetId)cell;
    // check on cells to be ignored and dropped: (rof,20.Feb.09)
    const HcalChannelStatus* mydigistatus=myqual->getValues(detcell.rawId());
    if (mySeverity->dropChannel(mydigistatus->getValue() ) ) continue;
    if (dropZSmarkedPassed_)
      if (i->zsMarkAndPass()) continue;

    const HcalCalibrations& calibrations=conditions->getHcalCalibrations(cell);
    const HcalQIECoder* channelCoder = conditions->getHcalCoder (cell);
    const HcalQIEShape* shape = conditions->getHcalShape (channelCoder);
    HcalCoderDb coder (*channelCoder, *shape);

    // firstSample & samplesToAdd
        if(tsFromDB_) {
      const HcalRecoParam* param_ts = paramTS->getValues(detcell.rawId());
      first = param_ts->firstSample();    
      toadd = param_ts->samplesToAdd();    
    }
    rec->push_back(reco_.reconstruct(*i,first,toadd,coder,calibrations));

    /*
      // Flag setting not available for calibration rechits
    // Set auxiliary flag
    int auxflag=0;
        int fTS = firstAuxTS_;
    for (int xx=fTS; xx<fTS+4 && xx<i->size();++xx)
      auxflag+=(i->sample(xx).adc())<<(7*(xx-fTS)); // store the time slices in the first 28 bits of aux, a set of 4 7-bit adc values
    // bits 28 and 29 are reserved for capid of the first time slice saved in aux
    auxflag+=((i->sample(fTS).capid())<<28);
    (rec->back()).setAux(auxflag);

    (rec->back()).setFlags(0); // Not yet implemented for HcalCalibRecHit
    */
      }
      // return result
      e.put(std::move(rec));     
    }
  } 
  //DL  delete myqual;
} // void HcalHitReconstructor::produce(...)