|
|
|
#include <EcalTBWeightUncalibRecHitProducer.h>
Public Types | |
| typedef std::vector< double > | EcalRecoAmplitudes |
Public Member Functions | |
| EcalTBWeightUncalibRecHitProducer (const edm::ParameterSet &ps) | |
| virtual void | produce (edm::Event &evt, const edm::EventSetup &es) |
| ~EcalTBWeightUncalibRecHitProducer () | |
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_ |
produce ECAL uncalibrated rechits from dataframes
$Alex Zabi$
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.
| typedef std::vector<double> EcalTBWeightUncalibRecHitProducer::EcalRecoAmplitudes |
Definition at line 23 of file EcalTBWeightUncalibRecHitProducer.h.
| 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.
{
}
| 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(), EBDataFrame::id(), EcalXtalGroupId::id(), EEDataFrame::id(), edm::InputTag::instance(), edm::InputTag::label(), LogDebug, EcalUncalibRecHitRecWeightsAlgo< C >::makeRecHit(), nbTimeBin_, edm::Event::put(), EcalDataFrame::sample(), EcalDataFrame::size(), tdcRecInfoCollection_, testbeamEBShape, testbeamEEShape, use2004OffsetConvention_, and ExpressReco_HICollisions_FallBack::weights.
{
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_ );
}
Definition at line 37 of file EcalTBWeightUncalibRecHitProducer.h.
Referenced by produce().
Definition at line 30 of file EcalTBWeightUncalibRecHitProducer.h.
Referenced by EcalTBWeightUncalibRecHitProducer(), and produce().
std::string EcalTBWeightUncalibRecHitProducer::EBhitCollection_ [private] |
Definition at line 34 of file EcalTBWeightUncalibRecHitProducer.h.
Referenced by EcalTBWeightUncalibRecHitProducer(), and produce().
Definition at line 38 of file EcalTBWeightUncalibRecHitProducer.h.
Referenced by produce().
Definition at line 31 of file EcalTBWeightUncalibRecHitProducer.h.
Referenced by EcalTBWeightUncalibRecHitProducer(), and produce().
std::string EcalTBWeightUncalibRecHitProducer::EEhitCollection_ [private] |
Definition at line 35 of file EcalTBWeightUncalibRecHitProducer.h.
Referenced by EcalTBWeightUncalibRecHitProducer(), and produce().
int EcalTBWeightUncalibRecHitProducer::nbTimeBin_ [private] |
Definition at line 46 of file EcalTBWeightUncalibRecHitProducer.h.
Referenced by EcalTBWeightUncalibRecHitProducer(), and produce().
Definition at line 32 of file EcalTBWeightUncalibRecHitProducer.h.
Referenced by EcalTBWeightUncalibRecHitProducer(), and produce().
const EBShape EcalTBWeightUncalibRecHitProducer::testbeamEBShape [private] |
Definition at line 41 of file EcalTBWeightUncalibRecHitProducer.h.
Referenced by produce().
const EEShape EcalTBWeightUncalibRecHitProducer::testbeamEEShape [private] |
Definition at line 40 of file EcalTBWeightUncalibRecHitProducer.h.
Referenced by produce().
bool EcalTBWeightUncalibRecHitProducer::use2004OffsetConvention_ [private] |
Definition at line 49 of file EcalTBWeightUncalibRecHitProducer.h.
Referenced by EcalTBWeightUncalibRecHitProducer(), and produce().
1.7.3