CMS 3D CMS Logo

EcalTBWeightUncalibRecHitProducer.cc
Go to the documentation of this file.
1 
15 
16 #include <iostream>
17 #include <iomanip>
18 #include <cmath>
19 
21 
23 
26 
29 
33 
38 
42 
43 #include "CLHEP/Matrix/Matrix.h"
44 #include "CLHEP/Matrix/SymMatrix.h"
45 #include <vector>
46 
47 #define DEBUG
49 
50  EBdigiCollection_ = ps.getParameter<edm::InputTag>("EBdigiCollection");
51  EEdigiCollection_ = ps.getParameter<edm::InputTag>("EEdigiCollection");
52  tdcRecInfoCollection_ = ps.getParameter<edm::InputTag>("tdcRecInfoCollection");
53  EBhitCollection_ = ps.getParameter<std::string>("EBhitCollection");
54  EEhitCollection_ = ps.getParameter<std::string>("EEhitCollection");
55  nbTimeBin_ = ps.getParameter<int>("nbTimeBin");
56  use2004OffsetConvention_ = ps.getUntrackedParameter< bool >("use2004OffsetConvention",false);
57  produces< EBUncalibratedRecHitCollection >(EBhitCollection_);
58  produces< EEUncalibratedRecHitCollection >(EEhitCollection_);
59 }
60 
62 }
63 
64 void
66 
67  using namespace edm;
68 
70  const EBDigiCollection* EBdigis =nullptr;
71  if (EBdigiCollection_.label() != "" || EBdigiCollection_.instance() != "")
72  {
73  // evt.getByLabel( digiProducer_, EBdigiCollection_, pEBDigis);
74  evt.getByLabel( EBdigiCollection_, pEBDigis);
75  if (!pEBDigis.isValid()) {
76  edm::LogError("EcalUncalibRecHitError") << "Error! can't get the product " << EBdigiCollection_ ;
77  } else {
78  EBdigis = pEBDigis.product(); // get a ptr to the produc
79 #ifdef DEBUG
80  LogDebug("EcalUncalibRecHitInfo") << "total # EBdigis: " << EBdigis->size() ;
81 #endif
82  }
83  }
84 
86  const EEDigiCollection* EEdigis =nullptr;
87 
88  if (EEdigiCollection_.label() != "" || EEdigiCollection_.instance() != "")
89  {
90  // evt.getByLabel( digiProducer_, EEdigiCollection_, pEEDigis);
91  evt.getByLabel( EEdigiCollection_, pEEDigis);
92  if (!pEEDigis.isValid()) {
93  edm::LogError("EcalUncalibRecHitError") << "Error! can't get the product " << EEdigiCollection_ ;
94  } else {
95  EEdigis = pEEDigis.product(); // get a ptr to the produc
96 #ifdef DEBUG
97  LogDebug("EcalUncalibRecHitInfo") << "total # EEdigis: " << EEdigis->size() ;
98  // std::cout << "total # EEdigis: " << EEdigis->size() << std::endl;
99 #endif
100  }
101  }
102 
103 
104  if (!EBdigis && !EEdigis)
105  return;
106 
108  const EcalTBTDCRecInfo* recTDC =nullptr;
109 
110  // evt.getByLabel( digiProducer_, EBdigiCollection_, pEBDigis);
111  evt.getByLabel( tdcRecInfoCollection_, pRecTDC);
112  if (pRecTDC.isValid()) {
113  recTDC = pRecTDC.product(); // get a ptr to the product
114  }
115 
116  // fetch map of groups of xtals
118  es.get<EcalWeightXtalGroupsRcd>().get(pGrp);
119  const EcalWeightXtalGroups* grp = pGrp.product();
120  if (!grp)
121  return;
122 
123  const EcalXtalGroupsMap& grpMap = grp->getMap();
124 
125 
126  // Gain Ratios
128  es.get<EcalGainRatiosRcd>().get(pRatio);
129  const EcalGainRatioMap& gainMap = pRatio.product()->getMap(); // map of gain ratios
130 
131 
132  // fetch TB weights
133 #ifdef DEBUG
134  LogDebug("EcalUncalibRecHitDebug") <<"Fetching EcalTBWeights from DB " ;
135  // std::cout <<"Fetching EcalTBWeights from DB " ;
136 #endif
138  es.get<EcalTBWeightsRcd>().get(pWgts);
139  const EcalTBWeights* wgts = pWgts.product();
140 
141  if (!wgts)
142  return;
143 
144 #ifdef DEBUG
145  LogDebug("EcalUncalibRecHitDebug") << "EcalTBWeightMap.size(): " << std::setprecision(3) << wgts->getMap().size() ;
146  // std::cout << "EcalTBWeightMap.size(): " << std::setprecision(3) << wgts->getMap().size() ;
147 #endif
148 
149  // fetch the pedestals from the cond DB via EventSetup
150 #ifdef DEBUG
151  LogDebug("EcalUncalibRecHitDebug") << "fetching pedestals....";
152  // std::cout << "fetching pedestals....";
153 #endif
155  es.get<EcalPedestalsRcd>().get( pedHandle );
156  const EcalPedestalsMap& pedMap = pedHandle.product()->getMap(); // map of pedestals
157 #ifdef DEBUG
158  LogDebug("EcalUncalibRecHitDebug") << "done." ;
159  // std::cout << "done." ;
160 #endif
161  // collection of reco'ed ampltudes to put in the event
162 
163  auto EBuncalibRechits = std::make_unique<EBUncalibratedRecHitCollection>();
164  auto EEuncalibRechits = std::make_unique<EEUncalibratedRecHitCollection>();
165 
166  EcalPedestalsMapIterator pedIter; // pedestal iterator
167 
168  EcalGainRatioMap::const_iterator gainIter; // gain iterator
169 
170  EcalXtalGroupsMap::const_iterator git; // group iterator
171 
172  EcalTBWeights::EcalTBWeightMap::const_iterator wit; //weights iterator
173  // loop over EB digis
174  //Getting the TDC bin
175  EcalTBWeights::EcalTDCId tdcid(int(nbTimeBin_/2)+1);
176 
177  if (recTDC)
178  if (recTDC->offset() == -999.)
179  {
180  edm::LogError("EcalUncalibRecHitError") << "TDC bin completely out of range. Returning" ;
181  return;
182  }
183 
184  if (EBdigis)
185  {
186  for(unsigned int idig = 0; idig < EBdigis->size(); ++idig) {
187 
188  EBDataFrame itdg = (*EBdigis)[idig];
189 
190 
191  // counter_++; // verbosity counter
192 
193  // find pedestals for this channel
194 #ifdef DEBUG
195  LogDebug("EcalUncalibRecHitDebug") << "looking up pedestal for crystal: " << EBDetId(itdg.id()) ;
196 #endif
197  pedIter = pedMap.find(itdg.id().rawId());
198  if( pedIter == pedMap.end() ) {
199  edm::LogError("EcalUncalibRecHitError") << "error!! could not find pedestals for channel: " << EBDetId(itdg.id())
200  << "\n no uncalib rechit will be made for this digi!"
201  ;
202  continue;
203  }
204  const EcalPedestals::Item& aped = (*pedIter);
205  double pedVec[3];
206  double pedRMSVec[3];
207  pedVec[0]=aped.mean_x12;pedVec[1]=aped.mean_x6;pedVec[2]=aped.mean_x1;
208  pedRMSVec[0]=aped.rms_x12;pedRMSVec[1]=aped.rms_x6;pedRMSVec[2]=aped.rms_x1;
209 
210  // find gain ratios
211 #ifdef DEBUG
212  LogDebug("EcalUncalibRecHitDebug") << "looking up gainRatios for crystal: " << EBDetId(itdg.id()) ;
213 #endif
214  gainIter = gainMap.find(itdg.id().rawId());
215  if( gainIter == gainMap.end() ) {
216  edm::LogError("EcalUncalibRecHitError") << "error!! could not find gain ratios for channel: " << EBDetId(itdg.id())
217  << "\n no uncalib rechit will be made for this digi!"
218  ;
219  continue;
220  }
221  const EcalMGPAGainRatio& aGain = (*gainIter);
222  double gainRatios[3];
223  gainRatios[0]=1.;gainRatios[1]=aGain.gain12Over6();gainRatios[2]=aGain.gain6Over1()*aGain.gain12Over6();
224 
225 
226 
227  // lookup group ID for this channel
228  git = grpMap.find( itdg.id().rawId() );
229  if( git == grpMap.end() ) {
230  edm::LogError("EcalUncalibRecHitError") << "No group id found for this crystal. something wrong with EcalWeightXtalGroups in your DB?"
231  << "\n no uncalib rechit will be made for digi with id: " << EBDetId(itdg.id())
232  ;
233  continue;
234  }
235  const EcalXtalGroupId& gid = (*git);
236 
237 
238  //GAIN SWITCHING DETECTION ///////////////////////////////////////////////////////////////////////////////////////////////////
239  double sampleGainRef = 1;
240  int sampleSwitch = 999;
241  for (int sample = 0; sample < itdg.size(); ++sample)
242  {
243  double gainSample = itdg.sample(sample).gainId();
244  if(gainSample != sampleGainRef) {sampleGainRef = gainSample; sampleSwitch = sample;}
245  }//loop sample
247 
248  if (recTDC)
249  {
250  int tdcBin=0;
251  if (recTDC->offset() <= 0.)
252  tdcBin = 1;
253  if (recTDC->offset() >= 1.)
254  tdcBin = nbTimeBin_;
255  else
256  tdcBin = int(recTDC->offset()*float(nbTimeBin_))+1;
257 
258  if (tdcBin < 1 || tdcBin > nbTimeBin_ )
259  {
260  edm::LogError("EcalUncalibRecHitError") << "TDC bin out of range " << tdcBin << " offset " << recTDC->offset();
261  continue;
262  }
263 
264  // In case gain switching happens at the sample 4 (5th sample)
265  // (sample 5 (6th sample) in 2004 TDC convention) an extra
266  // set of weights has to be used. This set of weights is assigned to
267  // TDC values going from 25 and up.
268  if (use2004OffsetConvention_ && sampleSwitch == 5)
269  tdcid=EcalTBWeights::EcalTDCId(tdcBin+25);
270  else if (!use2004OffsetConvention_ && sampleSwitch == 4)
271  tdcid=EcalTBWeights::EcalTDCId(tdcBin+25);
272  else
273  tdcid=EcalTBWeights::EcalTDCId(tdcBin);
274  }//check TDC
275 
276  // now lookup the correct weights in the map
277  wit = wgts->getMap().find( std::make_pair(gid,tdcid) );
278  if( wit == wgts->getMap().end() ) { // no weights found for this group ID
279  edm::LogError("EcalUncalibRecHitError") << "No weights found for EcalGroupId: " << gid.id() << " and EcalTDCId: " << tdcid
280  << "\n skipping digi with id: " << EBDetId(itdg.id())
281  ;
282  continue;
283  }
284  const EcalWeightSet& wset = wit->second; // this is the EcalWeightSet
285 
286 
287  // EcalWeightMatrix is vec<vec:double>>
288 
289 #ifdef DEBUG
290  LogDebug("EcalUncalibRecHitDebug") << "accessing matrices of weights...";
291 #endif
294  //const EcalWeightSet::EcalChi2WeightMatrix& mat3 = wset.getChi2WeightsBeforeGainSwitch();
295  //const EcalWeightSet::EcalChi2WeightMatrix& mat4 = wset.getChi2WeightsAfterGainSwitch();
297  weights[0]=&mat1;
298  weights[1]=&mat2;
299  // weights.push_back(clmat1);
300  // weights.push_back(clmat2);
301  // LogDebug("EcalUncalibRecHitDebug") << "weights before switch:\n" << clmat1 ;
302  // LogDebug("EcalUncalibRecHitDebug") << "weights after switch:\n" << clmat2 ;
303 
304 
305  // build CLHEP chi2 matrices
306  //const EcalWeightSet::EcalChi2WeightMatrix* chi2mat[2];
307  // chi2mat[0]=&mat3;
308  // chi2mat[1]=&mat4;
309 
311  EBalgo_.makeRecHit(itdg, pedVec, pedRMSVec, gainRatios, weights, testbeamEBShape);
312  //EBalgo_.makeRecHit(itdg, pedVec, gainRatios, weights, chi2mat);
313  EBuncalibRechits->push_back( aHit );
314 #ifdef DEBUG
315  if(aHit.amplitude()>0.) {
316  LogDebug("EcalUncalibRecHitDebug") << "processed EBDataFrame with id: "
317  << EBDetId(itdg.id()) << "\n"
318  << "uncalib rechit amplitude: " << aHit.amplitude()
319  ;
320  }
321 #endif
322  }
323  }
324  // put the collection of reconstructed hits in the event
325  evt.put(std::move(EBuncalibRechits), EBhitCollection_);
326 
327 
328  if (EEdigis)
329  {
330  for(unsigned int idig = 0; idig < EEdigis->size(); ++idig) {
331 
332  EEDataFrame itdg = (*EEdigis)[idig];
333 
334  // counter_++; // verbosity counter
335 
336  // find pedestals for this channel
337 #ifdef DEBUG
338  LogDebug("EcalUncalibRecHitDebug") << "looking up pedestal for crystal: " << EEDetId(itdg.id()) ;
339  // std::cout << "looking up pedestal for crystal: " << EEDetId(itdg.id()) ;
340 #endif
341  pedIter = pedMap.find(itdg.id().rawId());
342  if( pedIter == pedMap.end() ) {
343  edm::LogError("EcalUncalibRecHitError") << "error!! could not find pedestals for channel: " << EEDetId(itdg.id())
344  << "\n no uncalib rechit will be made for this digi!"
345  ;
346  continue;
347  }
348  const EcalPedestals::Item& aped = (*pedIter);
349  double pedVec[3];
350  double pedRMSVec[3];
351  pedVec[0]=aped.mean_x12;pedVec[1]=aped.mean_x6;pedVec[2]=aped.mean_x1;
352  pedRMSVec[0]=aped.rms_x12;pedRMSVec[1]=aped.rms_x6;pedRMSVec[2]=aped.rms_x1;
353 
354  // find gain ratios
355 #ifdef DEBUG
356  LogDebug("EcalUncalibRecHitDebug") << "looking up gainRatios for crystal: " << EEDetId(itdg.id()) ;
357  // std::cout << "looking up gainRatios for crystal: " << EEDetId(itdg.id()) ;
358 #endif
359  gainIter = gainMap.find(itdg.id().rawId());
360  if( gainIter == gainMap.end() ) {
361  edm::LogError("EcalUncalibRecHitError") << "error!! could not find gain ratios for channel: " << EEDetId(itdg.id())
362  << "\n no uncalib rechit will be made for this digi!"
363  ;
364  continue;
365  }
366  const EcalMGPAGainRatio& aGain = (*gainIter);
367  double gainRatios[3];
368  gainRatios[0]=1.;gainRatios[1]=aGain.gain12Over6();gainRatios[2]=aGain.gain6Over1()*aGain.gain12Over6();
369 
370 
371 
372  // lookup group ID for this channel
373  git = grpMap.find( itdg.id().rawId() );
374  if( git == grpMap.end() ) {
375  edm::LogError("EcalUncalibRecHitError") << "No group id found for this crystal. something wrong with EcalWeightXtalGroups in your DB?"
376  << "\n no uncalib rechit will be made for digi with id: " << EEDetId(itdg.id())
377  ;
378  continue;
379  }
380  const EcalXtalGroupId& gid = (*git);
381 
382 
383  //GAIN SWITCHING DETECTION ///////////////////////////////////////////////////////////////////////////////////////////////////
384  double sampleGainRef = 1;
385  int sampleSwitch = 999;
386  for (int sample = 0; sample < itdg.size(); ++sample)
387  {
388  double gainSample = itdg.sample(sample).gainId();
389  if(gainSample != sampleGainRef) {sampleGainRef = gainSample; sampleSwitch = sample;}
390  }//loop sample
392 
393  if (recTDC)
394  {
395  int tdcBin=0;
396  if (recTDC->offset() <= 0.)
397  tdcBin = 1;
398  if (recTDC->offset() >= 1.)
399  tdcBin = nbTimeBin_;
400  else
401  tdcBin = int(recTDC->offset()*float(nbTimeBin_))+1;
402 
403  if (tdcBin < 1 || tdcBin > nbTimeBin_ )
404  {
405  edm::LogError("EcalUncalibRecHitError") << "TDC bin out of range " << tdcBin << " offset " << recTDC->offset();
406  continue;
407  }
408 
409  // In case gain switching happens at the sample 4 (5th sample)
410  // (sample 5 (6th sample) in 2004 TDC convention) an extra
411  // set of weights has to be used. This set of weights is assigned to
412  // TDC values going from 25 and up.
413  if (use2004OffsetConvention_ && sampleSwitch == 5)
414  tdcid=EcalTBWeights::EcalTDCId(tdcBin+25);
415  else if (!use2004OffsetConvention_ && sampleSwitch == 4)
416  tdcid=EcalTBWeights::EcalTDCId(tdcBin+25);
417  else
418  tdcid=EcalTBWeights::EcalTDCId(tdcBin);
419  }//check TDC
420 
421  // now lookup the correct weights in the map
422  wit = wgts->getMap().find( std::make_pair(gid,tdcid) );
423  if( wit == wgts->getMap().end() ) { // no weights found for this group ID
424  edm::LogError("EcalUncalibRecHitError") << "No weights found for EcalGroupId: " << gid.id() << " and EcalTDCId: " << tdcid
425  << "\n skipping digi with id: " << EEDetId(itdg.id())
426  ;
427  continue;
428  }
429  const EcalWeightSet& wset = wit->second; // this is the EcalWeightSet
430 
431 
432  // EcalWeightMatrix is vec<vec:double>>
433 
434 #ifdef DEBUG
435  LogDebug("EcalUncalibRecHitDebug") << "accessing matrices of weights...";
436  // std::cout << "accessing matrices of weights...";
437 #endif
440  //const EcalWeightSet::EcalChi2WeightMatrix& mat3 = wset.getChi2WeightsBeforeGainSwitch();
441  //const EcalWeightSet::EcalChi2WeightMatrix& mat4 = wset.getChi2WeightsAfterGainSwitch();
443  weights[0]=&mat1;
444  weights[1]=&mat2;
445  // weights.push_back(clmat1);
446  // weights.push_back(clmat2);
447  // LogDebug("EcalUncalibRecHitDebug") << "weights before switch:\n" << clmat1 ;
448  // LogDebug("EcalUncalibRecHitDebug") << "weights after switch:\n" << clmat2 ;
449 
450 
451  // build CLHEP chi2 matrices
452  //const EcalWeightSet::EcalChi2WeightMatrix* chi2mat[2];
453  //chi2mat[0]=&mat3;
454  //chi2mat[1]=&mat4;
455 
457  EEalgo_.makeRecHit(itdg, pedVec, pedRMSVec, gainRatios, weights, testbeamEEShape);
458  //EEalgo_.makeRecHit(itdg, pedVec, gainRatios, weights, chi2mat);
459  EEuncalibRechits->push_back( aHit );
460 #ifdef DEBUG
461  if(aHit.amplitude()>0.) {
462  LogDebug("EcalUncalibRecHitDebug") << "processed EEDataFrame with id: "
463  << EEDetId(itdg.id()) << "\n"
464  << "uncalib rechit amplitude: " << aHit.amplitude()
465 
466 // std::cout << "processed EEDataFrame with id: "
467 // << EEDetId(itdg.id()) << "\n"
468 // << "uncalib rechit amplitude: " << aHit.amplitude() << std::endl;
469  ;
470  }
471 #endif
472  }
473  }
474  // put the collection of reconstructed hits in the event
475  evt.put(std::move(EEuncalibRechits), EEhitCollection_);
476 }
477 
478 // HepMatrix
479 // EcalTBWeightUncalibRecHitProducer::makeMatrixFromVectors(const std::vector< std::vector<EcalWeight> >& vecvec) {
480 // int nrow = vecvec.size();
481 // int ncol = (vecvec[0]).size();
482 // HepMatrix clmat(nrow,ncol);
483 // //LogDebug("EcalUncalibRecHitDebug") << "created HepMatrix(" << nrow << "," << ncol << ")" ;
484 // for(int irow=0;irow<nrow;++irow) {
485 // for(int icol=0;icol<ncol;++icol) {
486 // clmat[irow][icol] = ((vecvec[irow])[icol]).value();
487 // }
488 // }
489 // return clmat;
490 // }
491 
492 // HepMatrix
493 // EcalTBWeightUncalibRecHitProducer::makeDummySymMatrix(int size)
494 // {
495 // HepMatrix clmat(10,10);
496 // //LogDebug("EcalUncalibRecHitDebug") << "created HepMatrix(" << nrow << "," << ncol << ")" ;
497 // for(int irow=0; irow<size; ++irow) {
498 // for(int icol=0 ; icol<size; ++icol) {
499 // clmat[irow][icol] = irow+icol;
500 // }
501 // }
502 // return clmat;
503 // }
#define LogDebug(id)
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
Definition: Event.h:137
key_type id() const
Definition: EBDataFrame.h:31
const self & getMap() const
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:47
EcalMGPASample sample(int i) const
Definition: EcalDataFrame.h:29
int gainId() const
get the gainId (2 bits)
int size() const
Definition: EcalDataFrame.h:26
const unsigned int id() const
EcalPedestalsMap::const_iterator EcalPedestalsMapIterator
Definition: EcalPedestals.h:52
EcalTBWeightUncalibRecHitProducer(const edm::ParameterSet &ps)
virtual EcalUncalibratedRecHit makeRecHit(const C &dataFrame, const double *pedestals, const double *pedestalsRMS, const double *gainRatios, const EcalWeightSet::EcalWeightMatrix **weights, const EcalShapeBase &testbeamPulseShape)
Compute parameters.
float gain6Over1() const
bool isValid() const
Definition: HandleBase.h:74
EcalWeightMatrix & getWeightsAfterGainSwitch()
Definition: EcalWeightSet.h:30
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:535
key_type id() const
Definition: EEDataFrame.h:28
EcalUncalibRecHitRecWeightsAlgo< EEDataFrame > EEalgo_
EcalWeightMatrix & getWeightsBeforeGainSwitch()
Definition: EcalWeightSet.h:29
T const * product() const
Definition: Handle.h:81
std::vector< Item >::const_iterator const_iterator
float gain12Over6() const
std::string const & label() const
Definition: InputTag.h:36
HLT enums.
void produce(edm::Event &evt, const edm::EventSetup &es) override
T get() const
Definition: EventSetup.h:63
const_iterator find(uint32_t rawId) const
const EcalTBWeightMap & getMap() const
Definition: EcalTBWeights.h:31
math::Matrix< 3, 10 >::type EcalWeightMatrix
Definition: EcalWeightSet.h:22
const_iterator end() const
EcalUncalibRecHitRecWeightsAlgo< EBDataFrame > EBalgo_
T const * product() const
Definition: ESHandle.h:86
std::string const & instance() const
Definition: InputTag.h:37
def move(src, dest)
Definition: eostools.py:510
float offset() const