#include <EcalTBWeightUncalibRecHitProducer.h>

Inheritance diagram for EcalTBWeightUncalibRecHitProducer:

Public Types
typedef std::vector< double >	EcalRecoAmplitudes

Public Types inherited from edm::EDProducer
typedef EDProducer	ModuleType

typedef WorkerT< EDProducer >	WorkerType

Public Types inherited from edm::ProducerBase
typedef ProductRegistryHelper::TypeLabelList	TypeLabelList

Public Member Functions
	EcalTBWeightUncalibRecHitProducer (const edm::ParameterSet &ps)

virtual void	produce (edm::Event &evt, const edm::EventSetup &es)

	~EcalTBWeightUncalibRecHitProducer ()

Public Member Functions inherited from edm::EDProducer
	EDProducer ()

virtual	~EDProducer ()

Public Member Functions inherited from edm::ProducerBase
	ProducerBase ()

void	registerProducts (ProducerBase , ProductRegistry , ModuleDescription const &)

boost::function< void(const BranchDescription &)>	registrationCallback () const
	used by the fwk to register list of products More...

virtual	~ProducerBase ()

Private Attributes
EcalUncalibRecHitRecWeightsAlgo < EBDataFrame >	EBalgo_

edm::InputTag	EBdigiCollection_

std::string	EBhitCollection_

EcalUncalibRecHitRecWeightsAlgo < EEDataFrame >	EEalgo_

edm::InputTag	EEdigiCollection_

std::string	EEhitCollection_

int	nbTimeBin_

edm::InputTag	tdcRecInfoCollection_

const EBShape	testbeamEBShape

const EEShape	testbeamEEShape

bool	use2004OffsetConvention_

Additional Inherited Members
Static Public Member Functions inherited from edm::EDProducer
static const std::string &	baseType ()

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

Protected Member Functions inherited from edm::EDProducer
CurrentProcessingContext const *	currentContext () const

Protected Member Functions inherited from edm::ProducerBase
template<class TProducer , class TMethod >
void	callWhenNewProductsRegistered (TProducer *iProd, TMethod iMethod)

Detailed Description

produce ECAL uncalibrated rechits from dataframes

Id:: EcalTBWeightUncalibRecHitProducer.cc,v 1.13 2008/11/21 15:31:24 meridian Exp

Date:: 2008/11/21 15:31:24

Revision:: 1.13

$Alex Zabi$

Date:: 2008/11/21 15:31:24

Revision:: 1.13

Modification to detect first sample to switch gain. used for amplitude recontruction at high energy Add TDC convention option (P. Meridiani)

Definition at line 20 of file EcalTBWeightUncalibRecHitProducer.h.

Member Typedef Documentation

typedef std::vector<double> EcalTBWeightUncalibRecHitProducer::EcalRecoAmplitudes

Definition at line 23 of file EcalTBWeightUncalibRecHitProducer.h.

Constructor & Destructor Documentation

EcalTBWeightUncalibRecHitProducer::EcalTBWeightUncalibRecHitProducer ( const edm::ParameterSet & ps )

explicit

Definition at line 53 of file EcalTBWeightUncalibRecHitProducer.cc.

References EBdigiCollection_, EBhitCollection_, EEdigiCollection_, EEhitCollection_, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), nbTimeBin_, tdcRecInfoCollection_, and use2004OffsetConvention_.

                                                                                               {
 
    EBdigiCollection_ = ps.getParameter<edm::InputTag>("EBdigiCollection");
    EEdigiCollection_ = ps.getParameter<edm::InputTag>("EEdigiCollection");
    tdcRecInfoCollection_ = ps.getParameter<edm::InputTag>("tdcRecInfoCollection");
    EBhitCollection_  = ps.getParameter<std::string>("EBhitCollection");
    EEhitCollection_  = ps.getParameter<std::string>("EEhitCollection");
    nbTimeBin_  = ps.getParameter<int>("nbTimeBin");
    use2004OffsetConvention_ = ps.getUntrackedParameter< bool >("use2004OffsetConvention",false);
    produces< EBUncalibratedRecHitCollection >(EBhitCollection_);
    produces< EEUncalibratedRecHitCollection >(EEhitCollection_);
 }

EcalTBWeightUncalibRecHitProducer::~EcalTBWeightUncalibRecHitProducer ( )

Definition at line 66 of file EcalTBWeightUncalibRecHitProducer.cc.

66 {

67 }

Member Function Documentation

void EcalTBWeightUncalibRecHitProducer::produce	(	edm::Event &	evt,
		const edm::EventSetup &	es
	)

virtual

Implements edm::EDProducer.

Definition at line 70 of file EcalTBWeightUncalibRecHitProducer.cc.

References EcalUncalibratedRecHit::amplitude(), EBalgo_, EBdigiCollection_, EBhitCollection_, EEalgo_, EEdigiCollection_, EEhitCollection_, EcalMGPAGainRatio::gain12Over6(), EcalMGPAGainRatio::gain6Over1(), EcalMGPASample::gainId(), edm::EventSetup::get(), edm::Event::getByLabel(), EcalWeightSet::getChi2WeightsAfterGainSwitch(), EcalWeightSet::getChi2WeightsBeforeGainSwitch(), EcalCondObjectContainer< T >::getMap(), EcalWeightSet::getWeightsAfterGainSwitch(), EcalWeightSet::getWeightsBeforeGainSwitch(), EcalXtalGroupId::id(), EEDataFrame::id(), EBDataFrame::id(), edm::InputTag::instance(), edm::InputTag::label(), LogDebug, EcalUncalibRecHitRecWeightsAlgo< C >::makeRecHit(), nbTimeBin_, edm::Event::put(), DetId::rawId(), EcalDataFrame::sample(), EcalDataFrame::size(), tdcRecInfoCollection_, testbeamEBShape, testbeamEEShape, and use2004OffsetConvention_.

                                                                                  {
 
    using namespace edm;
    
    Handle< EBDigiCollection > pEBDigis;
    const EBDigiCollection* EBdigis =0;
    if (EBdigiCollection_.label() != "" || EBdigiCollection_.instance() != "")
      {
        //     evt.getByLabel( digiProducer_, EBdigiCollection_, pEBDigis);
        evt.getByLabel( EBdigiCollection_, pEBDigis);
        if (!pEBDigis.isValid()) {
      edm::LogError("EcalUncalibRecHitError") << "Error! can't get the product " << EBdigiCollection_ ;
        } else {
      EBdigis = pEBDigis.product(); // get a ptr to the produc
 #ifdef DEBUG
      LogDebug("EcalUncalibRecHitInfo") << "total # EBdigis: " << EBdigis->size() ;
 #endif
        }
      }
 
    Handle< EEDigiCollection > pEEDigis;
    const EEDigiCollection* EEdigis =0;
 
    if (EEdigiCollection_.label() != "" || EEdigiCollection_.instance() != "")
      {
        //     evt.getByLabel( digiProducer_, EEdigiCollection_, pEEDigis);
        evt.getByLabel( EEdigiCollection_, pEEDigis);
        if (!pEEDigis.isValid()) {
      edm::LogError("EcalUncalibRecHitError") << "Error! can't get the product " << EEdigiCollection_ ;
        } else {
      EEdigis = pEEDigis.product(); // get a ptr to the produc
 #ifdef DEBUG
      LogDebug("EcalUncalibRecHitInfo") << "total # EEdigis: " << EEdigis->size() ;
      //  std::cout << "total # EEdigis: " << EEdigis->size() << std::endl;
 #endif
        }
      }
 
 
    if (!EBdigis && !EEdigis)
      return;
 
    Handle< EcalTBTDCRecInfo > pRecTDC;
    const EcalTBTDCRecInfo* recTDC =0;
 
    //     evt.getByLabel( digiProducer_, EBdigiCollection_, pEBDigis);
    evt.getByLabel( tdcRecInfoCollection_, pRecTDC);
    if (pRecTDC.isValid()) {
      recTDC = pRecTDC.product(); // get a ptr to the product
    } 
 
    // fetch map of groups of xtals
    edm::ESHandle<EcalWeightXtalGroups> pGrp;
    es.get<EcalWeightXtalGroupsRcd>().get(pGrp);
    const EcalWeightXtalGroups* grp = pGrp.product();
    if (!grp)
      return;
 
    const EcalXtalGroupsMap& grpMap = grp->getMap();
    
    
    // Gain Ratios
    edm::ESHandle<EcalGainRatios> pRatio;
    es.get<EcalGainRatiosRcd>().get(pRatio);
    const EcalGainRatioMap& gainMap = pRatio.product()->getMap(); // map of gain ratios
 
 
    // fetch TB weights
 #ifdef DEBUG
    LogDebug("EcalUncalibRecHitDebug") <<"Fetching EcalTBWeights from DB " ;
    //   std::cout  <<"Fetching EcalTBWeights from DB " ;
 #endif
    edm::ESHandle<EcalTBWeights> pWgts;
    es.get<EcalTBWeightsRcd>().get(pWgts);
    const EcalTBWeights* wgts = pWgts.product();
 
    if (!wgts)
      return;
 
 #ifdef DEBUG
    LogDebug("EcalUncalibRecHitDebug") << "EcalTBWeightMap.size(): " << std::setprecision(3) << wgts->getMap().size() ;
    //   std::cout << "EcalTBWeightMap.size(): " << std::setprecision(3) << wgts->getMap().size() ;
 #endif   
 
    // fetch the pedestals from the cond DB via EventSetup
 #ifdef DEBUG
    LogDebug("EcalUncalibRecHitDebug") << "fetching pedestals....";
    //   std::cout << "fetching pedestals....";
 #endif
    edm::ESHandle<EcalPedestals> pedHandle;
    es.get<EcalPedestalsRcd>().get( pedHandle );
    const EcalPedestalsMap& pedMap = pedHandle.product()->getMap(); // map of pedestals
 #ifdef DEBUG
    LogDebug("EcalUncalibRecHitDebug") << "done." ;
    //   std::cout << "done." ;
 #endif
    // collection of reco'ed ampltudes to put in the event
 
    std::auto_ptr< EBUncalibratedRecHitCollection > EBuncalibRechits( new EBUncalibratedRecHitCollection );
    std::auto_ptr< EEUncalibratedRecHitCollection > EEuncalibRechits( new EEUncalibratedRecHitCollection );
 
    EcalPedestalsMapIterator pedIter; // pedestal iterator
 
    EcalGainRatioMap::const_iterator gainIter; // gain iterator
 
    EcalXtalGroupsMap::const_iterator git; // group iterator
 
    EcalTBWeights::EcalTBWeightMap::const_iterator wit; //weights iterator 
    // loop over EB digis
    //Getting the TDC bin
    EcalTBWeights::EcalTDCId tdcid(int(nbTimeBin_/2)+1);
    
    if (recTDC)
      if (recTDC->offset() == -999.)
        {
      edm::LogError("EcalUncalibRecHitError") << "TDC bin completely out of range. Returning" ;
      return;
        }
 
    if (EBdigis)
      {
        for(unsigned int idig = 0; idig < EBdigis->size(); ++idig) {
      
      EBDataFrame itdg = (*EBdigis)[idig];
 
      
      //     counter_++; // verbosity counter
      
      // find pedestals for this channel
 #ifdef DEBUG
      LogDebug("EcalUncalibRecHitDebug") << "looking up pedestal for crystal: " << EBDetId(itdg.id()) ;
 #endif
      pedIter = pedMap.find(itdg.id().rawId());
      if( pedIter == pedMap.end() ) {
        edm::LogError("EcalUncalibRecHitError") << "error!! could not find pedestals for channel: " << EBDetId(itdg.id()) 
                            << "\n  no uncalib rechit will be made for this digi!"
          ;
        continue;
      }
      const EcalPedestals::Item& aped = (*pedIter);
      double pedVec[3];
      double pedRMSVec[3];
      pedVec[0]=aped.mean_x12;pedVec[1]=aped.mean_x6;pedVec[2]=aped.mean_x1;
      pedRMSVec[0]=aped.rms_x12;pedRMSVec[1]=aped.rms_x6;pedRMSVec[2]=aped.rms_x1;
 
      // find gain ratios
 #ifdef DEBUG
      LogDebug("EcalUncalibRecHitDebug") << "looking up gainRatios for crystal: " << EBDetId(itdg.id()) ;
 #endif
      gainIter = gainMap.find(itdg.id().rawId());
      if( gainIter == gainMap.end() ) {
        edm::LogError("EcalUncalibRecHitError") << "error!! could not find gain ratios for channel: " << EBDetId(itdg.id()) 
                            << "\n  no uncalib rechit will be made for this digi!"
          ;
        continue;
      }
      const EcalMGPAGainRatio& aGain = (*gainIter);
      double gainRatios[3];
      gainRatios[0]=1.;gainRatios[1]=aGain.gain12Over6();gainRatios[2]=aGain.gain6Over1()*aGain.gain12Over6();
      
      
      
      // lookup group ID for this channel
      git = grpMap.find( itdg.id().rawId() );
      if( git == grpMap.end() ) {
        edm::LogError("EcalUncalibRecHitError") << "No group id found for this crystal. something wrong with EcalWeightXtalGroups in your DB?"
                            << "\n  no uncalib rechit will be made for digi with id: " << EBDetId(itdg.id())
          ;
        continue;
      }
      const EcalXtalGroupId& gid = (*git);
 
 
      //GAIN SWITCHING DETECTION ///////////////////////////////////////////////////////////////////////////////////////////////////     
      double sampleGainRef = 1;
      int    sampleSwitch  = 999;
      for (int sample = 0; sample < itdg.size(); ++sample)
        {
          double gainSample = itdg.sample(sample).gainId();
          if(gainSample != sampleGainRef) {sampleGainRef = gainSample; sampleSwitch = sample;}
        }//loop sample
      
      if (recTDC)
        {
          int tdcBin=0;
          if (recTDC->offset() <= 0.)
            tdcBin = 1;
          if (recTDC->offset() >= 1.)
            tdcBin = nbTimeBin_;
          else
            tdcBin = int(recTDC->offset()*float(nbTimeBin_))+1;
        
          if (tdcBin < 1 || tdcBin > nbTimeBin_ )
            {
          edm::LogError("EcalUncalibRecHitError") << "TDC bin out of range " << tdcBin << " offset " << recTDC->offset();
          continue;
            }
 
          // In case gain switching happens at the sample 4 (5th sample) 
          // (sample 5 (6th sample) in 2004 TDC convention) an extra
          // set of weights has to be used. This set of weights is assigned to 
          // TDC values going from 25 and up.
          if (use2004OffsetConvention_ && sampleSwitch == 5)
            tdcid=EcalTBWeights::EcalTDCId(tdcBin+25);
          else if (!use2004OffsetConvention_ && sampleSwitch == 4)
            tdcid=EcalTBWeights::EcalTDCId(tdcBin+25);
          else 
            tdcid=EcalTBWeights::EcalTDCId(tdcBin);
        }//check TDC     
      
      // now lookup the correct weights in the map
      wit = wgts->getMap().find( std::make_pair(gid,tdcid) );
      if( wit == wgts->getMap().end() ) {  // no weights found for this group ID
        edm::LogError("EcalUncalibRecHitError") << "No weights found for EcalGroupId: " << gid.id() << " and  EcalTDCId: " << tdcid
                            << "\n  skipping digi with id: " << EBDetId(itdg.id())
          ;
        continue;
      }
      const EcalWeightSet& wset = wit->second; // this is the EcalWeightSet
 
 
      // EcalWeightMatrix is vec<vec:double>>
 
 #ifdef DEBUG
      LogDebug("EcalUncalibRecHitDebug") << "accessing matrices of weights...";
 #endif
      const EcalWeightSet::EcalWeightMatrix& mat1 = wset.getWeightsBeforeGainSwitch();
      const EcalWeightSet::EcalWeightMatrix& mat2 = wset.getWeightsAfterGainSwitch();
      const EcalWeightSet::EcalChi2WeightMatrix& mat3 = wset.getChi2WeightsBeforeGainSwitch();
      const EcalWeightSet::EcalChi2WeightMatrix& mat4 = wset.getChi2WeightsAfterGainSwitch();
      const EcalWeightSet::EcalWeightMatrix* weights[2];
      weights[0]=&mat1;
      weights[1]=&mat2;
      //      weights.push_back(clmat1);
      //      weights.push_back(clmat2);
      //      LogDebug("EcalUncalibRecHitDebug") << "weights before switch:\n" << clmat1 ;
      //      LogDebug("EcalUncalibRecHitDebug") << "weights after switch:\n" << clmat2 ;
      
      
      // build CLHEP chi2  matrices
      const EcalWeightSet::EcalChi2WeightMatrix* chi2mat[2];
      chi2mat[0]=&mat3;
      chi2mat[1]=&mat4;
 
      EcalUncalibratedRecHit aHit =
        EBalgo_.makeRecHit(itdg, pedVec, pedRMSVec, gainRatios, weights, testbeamEBShape);
        //EBalgo_.makeRecHit(itdg, pedVec, gainRatios, weights, chi2mat);
      EBuncalibRechits->push_back( aHit );
 #ifdef DEBUG
      if(aHit.amplitude()>0.) {
        LogDebug("EcalUncalibRecHitDebug") << "processed EBDataFrame with id: "
                           << EBDetId(itdg.id()) << "\n"
                           << "uncalib rechit amplitude: " << aHit.amplitude()
          ;
      }
 #endif
        }
      }
    // put the collection of recunstructed hits in the event
    evt.put( EBuncalibRechits, EBhitCollection_ );
 
 
    if (EEdigis)
      {
        for(unsigned int idig = 0; idig < EEdigis->size(); ++idig) {
 
      EEDataFrame itdg = (*EEdigis)[idig];
      
      //     counter_++; // verbosity counter
      
      // find pedestals for this channel
 #ifdef DEBUG
      LogDebug("EcalUncalibRecHitDebug") << "looking up pedestal for crystal: " << EEDetId(itdg.id()) ;
      //  std::cout << "looking up pedestal for crystal: " << EEDetId(itdg.id()) ;
 #endif
      pedIter = pedMap.find(itdg.id().rawId());
      if( pedIter == pedMap.end() ) {
        edm::LogError("EcalUncalibRecHitError") << "error!! could not find pedestals for channel: " << EEDetId(itdg.id()) 
                            << "\n  no uncalib rechit will be made for this digi!"
          ;
        continue;
      }
      const EcalPedestals::Item& aped = (*pedIter);
      double pedVec[3];
      double pedRMSVec[3];
      pedVec[0]=aped.mean_x12;pedVec[1]=aped.mean_x6;pedVec[2]=aped.mean_x1;
      pedRMSVec[0]=aped.rms_x12;pedRMSVec[1]=aped.rms_x6;pedRMSVec[2]=aped.rms_x1;
 
      // find gain ratios
 #ifdef DEBUG
      LogDebug("EcalUncalibRecHitDebug") << "looking up gainRatios for crystal: " << EEDetId(itdg.id()) ;
      //  std::cout << "looking up gainRatios for crystal: " << EEDetId(itdg.id()) ;
 #endif
      gainIter = gainMap.find(itdg.id().rawId());
      if( gainIter == gainMap.end() ) {
        edm::LogError("EcalUncalibRecHitError") << "error!! could not find gain ratios for channel: " << EEDetId(itdg.id()) 
                            << "\n  no uncalib rechit will be made for this digi!"
          ;
        continue;
      }
      const EcalMGPAGainRatio& aGain = (*gainIter);
      double gainRatios[3];
      gainRatios[0]=1.;gainRatios[1]=aGain.gain12Over6();gainRatios[2]=aGain.gain6Over1()*aGain.gain12Over6();
      
      
      
      // lookup group ID for this channel
      git = grpMap.find( itdg.id().rawId() );
      if( git == grpMap.end() ) {
        edm::LogError("EcalUncalibRecHitError") << "No group id found for this crystal. something wrong with EcalWeightXtalGroups in your DB?"
                            << "\n  no uncalib rechit will be made for digi with id: " << EEDetId(itdg.id())
          ;
        continue;
      }
      const EcalXtalGroupId& gid = (*git);
 
 
      //GAIN SWITCHING DETECTION ///////////////////////////////////////////////////////////////////////////////////////////////////     
      double sampleGainRef = 1;
      int    sampleSwitch  = 999;
      for (int sample = 0; sample < itdg.size(); ++sample)
        {
          double gainSample = itdg.sample(sample).gainId();
          if(gainSample != sampleGainRef) {sampleGainRef = gainSample; sampleSwitch = sample;}
        }//loop sample
      
      if (recTDC)
        {
          int tdcBin=0;
          if (recTDC->offset() <= 0.)
            tdcBin = 1;
          if (recTDC->offset() >= 1.)
            tdcBin = nbTimeBin_;
          else
            tdcBin = int(recTDC->offset()*float(nbTimeBin_))+1;
        
          if (tdcBin < 1 || tdcBin > nbTimeBin_ )
            {
          edm::LogError("EcalUncalibRecHitError") << "TDC bin out of range " << tdcBin << " offset " << recTDC->offset();
          continue;
            }
 
          // In case gain switching happens at the sample 4 (5th sample) 
          // (sample 5 (6th sample) in 2004 TDC convention) an extra
          // set of weights has to be used. This set of weights is assigned to 
          // TDC values going from 25 and up.
          if (use2004OffsetConvention_ && sampleSwitch == 5)
            tdcid=EcalTBWeights::EcalTDCId(tdcBin+25);
          else if (!use2004OffsetConvention_ && sampleSwitch == 4)
            tdcid=EcalTBWeights::EcalTDCId(tdcBin+25);
          else 
            tdcid=EcalTBWeights::EcalTDCId(tdcBin);
        }//check TDC     
      
      // now lookup the correct weights in the map
      wit = wgts->getMap().find( std::make_pair(gid,tdcid) );
      if( wit == wgts->getMap().end() ) {  // no weights found for this group ID
        edm::LogError("EcalUncalibRecHitError") << "No weights found for EcalGroupId: " << gid.id() << " and  EcalTDCId: " << tdcid
                            << "\n  skipping digi with id: " << EEDetId(itdg.id())
          ;
        continue;
      }
      const EcalWeightSet& wset = wit->second; // this is the EcalWeightSet
 
 
      // EcalWeightMatrix is vec<vec:double>>
 
 #ifdef DEBUG
      LogDebug("EcalUncalibRecHitDebug") << "accessing matrices of weights...";
      //  std::cout << "accessing matrices of weights...";
 #endif
      const EcalWeightSet::EcalWeightMatrix& mat1 = wset.getWeightsBeforeGainSwitch();
      const EcalWeightSet::EcalWeightMatrix& mat2 = wset.getWeightsAfterGainSwitch();
      const EcalWeightSet::EcalChi2WeightMatrix& mat3 = wset.getChi2WeightsBeforeGainSwitch();
      const EcalWeightSet::EcalChi2WeightMatrix& mat4 = wset.getChi2WeightsAfterGainSwitch();
      const EcalWeightSet::EcalWeightMatrix* weights[2];
      weights[0]=&mat1;
      weights[1]=&mat2;
      //      weights.push_back(clmat1);
      //      weights.push_back(clmat2);
      //      LogDebug("EcalUncalibRecHitDebug") << "weights before switch:\n" << clmat1 ;
      //      LogDebug("EcalUncalibRecHitDebug") << "weights after switch:\n" << clmat2 ;
      
      
      // build CLHEP chi2  matrices
      const EcalWeightSet::EcalChi2WeightMatrix* chi2mat[2];
      chi2mat[0]=&mat3;
      chi2mat[1]=&mat4;
 
      EcalUncalibratedRecHit aHit =
        EEalgo_.makeRecHit(itdg, pedVec, pedRMSVec, gainRatios, weights, testbeamEEShape);
        //EEalgo_.makeRecHit(itdg, pedVec, gainRatios, weights, chi2mat);
      EEuncalibRechits->push_back( aHit );
 #ifdef DEBUG
      if(aHit.amplitude()>0.) {
        LogDebug("EcalUncalibRecHitDebug") << "processed EEDataFrame with id: "
                           << EEDetId(itdg.id()) << "\n"
                           << "uncalib rechit amplitude: " << aHit.amplitude()
 
 //     std::cout << "processed EEDataFrame with id: "
 //                        << EEDetId(itdg.id()) << "\n"
 //                        << "uncalib rechit amplitude: " << aHit.amplitude() << std::endl;
          ;
      }
 #endif
        }
      }
    // put the collection of recunstructed hits in the event
    evt.put( EEuncalibRechits, EEhitCollection_ );
 }