CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HcalRecHitsMaker.cc
Go to the documentation of this file.
3 
15 #include "CLHEP/GenericFunctions/Erf.hh"
26 
27 #include "TFile.h"
28 #include "TGraph.h"
29 #include "TROOT.h"
30 #include <fstream>
31 
32 class RandomEngine;
33 
34 bool HcalRecHitsMaker::initialized_ = false;
35 std::vector<float> HcalRecHitsMaker::peds_;
36 std::vector<int> HcalRecHitsMaker::fctoadc_;
37 std::vector<float> HcalRecHitsMaker::sat_;
38 std::vector<float> HcalRecHitsMaker::gains_;
39 std::vector<float> HcalRecHitsMaker::noisesigma_;
40 std::vector<float> HcalRecHitsMaker::TPGFactor_;
41 std::vector<float> HcalRecHitsMaker::miscalib_;
42 std::vector<HcalDetId> HcalRecHitsMaker::theDetIds_;
43 std::vector<int> HcalRecHitsMaker::hbhi_;
44 std::vector<int> HcalRecHitsMaker::hehi_;
45 std::vector<int> HcalRecHitsMaker::hohi_;
46 std::vector<int> HcalRecHitsMaker::hfhi_;
47 unsigned HcalRecHitsMaker::maxIndex_ = 0 ;
48 
50  const RandomEngine * myrandom)
51  :
52  det_(det),
53  doDigis_(false),
54  noiseFromDb_(false),
55  random_(myrandom)
56  //,myHcalSimParameterMap_(0)
57 {
58  edm::ParameterSet RecHitsParameters=p.getParameter<edm::ParameterSet>("HCAL");
59  noise_ = RecHitsParameters.getParameter<std::vector<double> >("Noise");
60  threshold_ = RecHitsParameters.getParameter<std::vector<double> >("Threshold");
61  doSaturation_ = RecHitsParameters.getParameter<bool>("EnableSaturation");
62 
63  refactor_ = RecHitsParameters.getParameter<double> ("Refactor");
64  refactor_mean_ = RecHitsParameters.getParameter<double> ("Refactor_mean");
65  hcalfileinpath_= RecHitsParameters.getParameter<std::string> ("fileNameHcal");
66  inputCol_=RecHitsParameters.getParameter<edm::InputTag>("MixedSimHits");
68 
69  if(det_==4)
70  {
71  hbhi_.reserve(2600);
72  hehi_.reserve(2600);
73  }
74  else if (det_==5)
75  hohi_.reserve(2200);
76  else if (det_==6)
77  hfhi_.reserve(1800);
78 
79  if(threshold_.size()!=noise_.size())
80  {
81  edm::LogWarning("CaloRecHitsProducer") << " WARNING : HCAL Noise simulation, the number of parameters should be the same for the noise and the thresholds. Disabling the noise simulation " << std::endl;
82  noise_.clear();
83  noise_.push_back(0.);
84  }
85  else
86  nnoise_=noise_.size();
87 
88  // edm::ParameterSet hcalparam = p2.getParameter<edm::ParameterSet>("hcalSimParam");
89  // myHcalSimParameterMap_ = new HcalSimParameterMap(hcalparam);
90 
91  // Computes the fraction of HCAL above the threshold
92  Genfun::Erf myErf;
93  hcalHotFraction_.resize(nnoise_,0.);
95  if(noise_.size()>0)
96  {
97  for(unsigned inoise=0;inoise<nnoise_;++inoise)
98  {
99  if(noise_[inoise]==0)
100  {
101  hcalHotFraction_[inoise]=0.;
102  continue;
103  }
104  else if(noise_[inoise]==-1) {
105  noiseFromDb_=true;
106  continue;
107  }
108  else
109  {
110  hcalHotFraction_.push_back(0.5-0.5*myErf(threshold_[inoise]/noise_[inoise]/sqrt(2.)));
111  myGaussianTailGenerators_[inoise]=new GaussianTail(random_,noise_[inoise],threshold_[inoise]);
112  }
113  }
114  }
115 }
116 
118 {
119  clean();
120  if(myGaussianTailGenerators_.size())
121  {
122  for(unsigned igt=0; igt<myGaussianTailGenerators_.size();++igt)
123  delete myGaussianTailGenerators_[igt];
124  }
126  theDetIds_.clear();
127  hbhi_.clear();
128  hehi_.clear();
129  hohi_.clear();
130  hfhi_.clear();
131 
132 }
133 
134 void HcalRecHitsMaker::init(const edm::EventSetup &es,bool doDigis,bool doMiscalib)
135 {
136  doDigis_=doDigis;
137  doMiscalib_=doMiscalib;
138 // needed for the noise simulation
139  edm::ESHandle<HcalDbService> conditions;
140  es.get<HcalDbRecord>().get(conditions);
141  const HcalDbService * theDbService=conditions.product();
142 
143  // get the correction factors
145  es.get<HcalRespCorrsRcd>().get(rchandle);
146  myRespCorr= rchandle.product();
147 
148 
149  if(!initialized_)
150  {
151  theDetIds_.resize(9201);
152  unsigned ncells=createVectorsOfCells(es);
153  edm::LogInfo("CaloRecHitsProducer") << "Total number of cells in HCAL " << ncells << std::endl;
154  hcalRecHits_.resize(maxIndex_+1,0.);
155  edm::LogInfo("CaloRecHitsProducer") << "Largest HCAL hashedindex" << maxIndex_ << std::endl;
156 
157  peds_.resize(9201);
158  gains_.resize(9201);
159  if(doSaturation_)
160  sat_.resize(9201);
161  if(noiseFromDb_)
162  noisesigma_.resize(9201);
163 
164 
165 
166  miscalib_.resize(maxIndex_+1,1.);
167  // Read from file ( a la HcalRecHitsRecalib.cc)
168  // here read them from xml (particular to HCAL)
169  CaloMiscalibMapHcal mapHcal;
170  mapHcal.prefillMap();
171 
172 
173  edm::FileInPath hcalfiletmp("CalibCalorimetry/CaloMiscalibTools/data/"+hcalfileinpath_);
174  std::string hcalfile=hcalfiletmp.fullPath();
175  MiscalibReaderFromXMLHcal hcalreader_(mapHcal);
176  if(doMiscalib_ && !hcalfile.empty())
177  {
178  hcalreader_.parseXMLMiscalibFile(hcalfile);
179  // mapHcal.print();
180  std::map<uint32_t,float>::const_iterator it=mapHcal.get().begin();
181  std::map<uint32_t,float>::const_iterator itend=mapHcal.get().end();
182  for(;it!=itend;++it)
183  {
184  HcalDetId myDetId(it->first);
185  float icalconst=it->second;
187  }
188  }
189 
190 
191  // Open the histogram for the fC to ADC conversion
192  gROOT->cd();
193  edm::FileInPath myDataFile("FastSimulation/CaloRecHitsProducer/data/adcvsfc.root");
194  TFile * myFile = new TFile(myDataFile.fullPath().c_str(),"READ");
195  TGraph * myGraf = (TGraph*)myFile->Get("adcvsfc");
196  unsigned size=myGraf->GetN();
197  fctoadc_.resize(10000);
198  unsigned p_index=0;
199  fctoadc_[0]=0;
200  int prev_nadc=0;
201  int nadc=0;
202  for(unsigned ibin=0;ibin<size;++ibin)
203  {
204  double x,y;
205  myGraf->GetPoint(ibin,x,y);
206  int index=(int)floor(x);
207  if(index<0||index>=10000) continue;
208  prev_nadc=nadc;
209  nadc=(int)y;
210  // Now fills the vector
211  for(unsigned ivec=p_index;ivec<(unsigned)index;++ivec)
212  {
213  fctoadc_[ivec] = prev_nadc;
214  }
215  p_index = index;
216  }
217  myFile->Close();
218  gROOT->cd();
219  edm::FileInPath myTPGFilePath("CalibCalorimetry/HcalTPGAlgos/data/RecHit-TPG-calib.dat");
220  TPGFactor_.resize(87,1.2);
221  std::ifstream myTPGFile(myTPGFilePath.fullPath().c_str(),ifstream::in);
222  if(myTPGFile)
223  {
224  float gain;
225  myTPGFile >> gain;
226  for(unsigned i=0;i<86;++i)
227  {
228  myTPGFile >> TPGFactor_[i] ;
229  // std::cout << TPGFactor_[i] << std::endl;
230  }
231  }
232  else
233  {
234  std::cout << " Unable to open CalibCalorimetry/HcalTPGAlgos/data/RecHit-TPG-calib.dat" << std::endl;
235  std::cout << " Using a constant 1.2 factor " << std::endl;
236  }
237  //HB
238  for(unsigned ic=0;ic<nhbcells_;++ic)
239  {
240  float gain = theDbService->getGain(theDetIds_[hbhi_[ic]])->getValue(0);
241  float mgain=0.;
242  for(unsigned ig=0;ig<4;++ig)
243  mgain+=theDbService->getGain(theDetIds_[hbhi_[ic]])->getValue(ig);
244  if(noiseFromDb_)
245  noisesigma_[hbhi_[ic]]=noiseInfCfromDB(theDbService,theDetIds_[hbhi_[ic]])*mgain*0.25;
246  // std::cout << " NOISEHB " << theDetIds_[hbhi_[ic]].ieta() << " " << noisesigma_[hbhi_[ic]] << " "<< std::endl;
247  // 10000 (ADC scale) / 4. (to compute the mean) / 0.92 ADC/fC
248  // *1.25 (only ~80% in 1ts Digi, while saturation applied to 4ts RecHit)
249  int ieta=theDetIds_[hbhi_[ic]].ieta();
250  float tpgfactor=TPGFactor_[(ieta>0)?ieta+43:-ieta];
251  mgain*=2500./0.92*tpgfactor ;// 10000 (ADC scale) / 4. (to compute the mean)
252  sat_[hbhi_[ic]]=(doSaturation_)?mgain:99999.;
253 
254  peds_[hbhi_[ic]]=theDbService->getPedestal(theDetIds_[hbhi_[ic]])->getValue(0);
255 
256  gain*=tpgfactor;
257  gains_[hbhi_[ic]]=gain;
258  }
259  //HE
260 
261  for(unsigned ic=0;ic<nhecells_;++ic)
262  {
263  float gain= theDbService->getGain(theDetIds_[hehi_[ic]])->getValue(0);
264  int ieta=theDetIds_[hehi_[ic]].ieta();
265  float mgain=0.;
266  for(unsigned ig=0;ig<4;++ig)
267  mgain+=theDbService->getGain(theDetIds_[hehi_[ic]])->getValue(ig);
268  if(noiseFromDb_)
269  noisesigma_[hehi_[ic]]=noiseInfCfromDB(theDbService,theDetIds_[hehi_[ic]])*mgain*0.25;
270 
271  // std::cout << " NOISEHE " << theDetIds_[hehi_[ic]].ieta() << " " << noisesigma_[hehi_[ic]] << " "<< std::endl;
272  float tpgfactor=TPGFactor_[(ieta>0)?ieta+44:-ieta+1];
273  mgain*=2500./0.92*tpgfactor;
274  sat_[hehi_[ic]]=(doSaturation_)?mgain:99999.;
275 
276  gain*=tpgfactor;
277  peds_[hehi_[ic]]=theDbService->getPedestal(theDetIds_[hehi_[ic]])->getValue(0);
278  gains_[hehi_[ic]]=gain;
279  }
280  //HO
281 
282  for(unsigned ic=0;ic<nhocells_;++ic)
283  {
284  float ped=theDbService->getPedestal(theDetIds_[hohi_[ic]])->getValue(0);
285  float gain=theDbService->getGain(theDetIds_[hohi_[ic]])->getValue(0);
286  float mgain=0.;
287  for(unsigned ig=0;ig<4;++ig)
288  mgain+=theDbService->getGain(theDetIds_[hohi_[ic]])->getValue(ig);
289  if(noiseFromDb_)
290  noisesigma_[hohi_[ic]]=noiseInfCfromDB(theDbService,theDetIds_[hohi_[ic]])*mgain*0.25;
291  // std::cout << " NOISEHO " << theDetIds_[hohi_[ic]].ieta() << " " << noisesigma_[hohi_[ic]] << " "<< std::endl;
292  int ieta=HcalDetId(hohi_[ic]).ieta();
293  float tpgfactor=TPGFactor_[(ieta>0)?ieta+43:-ieta];
294  mgain*=2500./0.92*tpgfactor;
295  sat_[hohi_[ic]]=(doSaturation_)?mgain:99999.;
296 
297  gain*=tpgfactor;
298  peds_[hohi_[ic]]=ped;
299  gains_[hohi_[ic]]=gain;
300  }
301  //HF
302 
303  for(unsigned ic=0;ic<nhfcells_;++ic)
304  {
305  float ped=theDbService->getPedestal(theDetIds_[hfhi_[ic]])->getValue(0);
306  float gain=theDbService->getGain(theDetIds_[hfhi_[ic]])->getValue(0);
307  float mgain=0.;
308  for(unsigned ig=0;ig<4;++ig)
309  mgain+=theDbService->getGain(theDetIds_[hfhi_[ic]])->getValue(ig);
310  // additional 1/2 factor for the HF (Salavat)
311  if(noiseFromDb_)
312  {
313  // computation from when the noise was taken
314  noisesigma_[hfhi_[ic]]=noiseInfCfromDB(theDbService,theDetIds_[hfhi_[ic]])*mgain*0.25;
315  }
316  // std::cout << " NOISEHF " << theDetIds_[hfhi_[ic]].ieta() << " " << noisesigma_[hfhi_[ic]] << " "<< std::endl;
317 
318  mgain*=2500./0.36;
319  sat_[hfhi_[ic]]=(doSaturation_)?mgain:99999.;
320  int ieta=theDetIds_[hfhi_[ic]].ieta();
321  gain*=TPGFactor_[(ieta>0)?ieta+45:-ieta+2];
322  peds_[hfhi_[ic]]=ped;
323  gains_[hfhi_[ic]]=gain;
324  }
325  initialized_=true;
326  }
327 
328  // clear the vector we don't need. It is a bit stupid
329  hcalRecHits_.resize(maxIndex_+1,0.);
330 
331 }
332 
333 
334 // Get the PCaloHits from the event. They have to be stored in a map, because when
335 // the pile-up is added thanks to the Mixing Module, the same cell can be present several times
337 {
338  clean();
339 
341  iEvent.getByLabel(inputCol_,cf);
342  std::auto_ptr<MixCollection<PCaloHit> > colcalo(new MixCollection<PCaloHit>(cf.product(),std::pair<int,int>(0,0) ));
343 
344  MixCollection<PCaloHit>::iterator it=colcalo->begin();;
345  MixCollection<PCaloHit>::iterator itend=colcalo->end();
346  unsigned counter=0;
347  for(;it!=itend;++it)
348  {
349  HcalDetId detid(it->id());
350  int hashedindex=detid.hashed_index();
351  switch(detid.subdet())
352  {
353  case HcalBarrel:
354  {
355  if(det_==4)
356  Fill(hashedindex,it->energy(),firedCells_,noise_[0]);
357  }
358  break;
359  case HcalEndcap:
360  {
361  if(det_==4)
362  Fill(hashedindex,it->energy(),firedCells_,noise_[1]);
363  }
364  break;
365  case HcalOuter:
366  {
367  if(det_==5)
368  Fill(hashedindex,it->energy(),firedCells_,noise_[0]);
369  }
370  break;
371  case HcalForward:
372  {
373  if(det_==6)
374  Fill(hashedindex,it->energy(),firedCells_,noise_[0]);
375  }
376  break;
377  default:
378  edm::LogWarning("CaloRecHitsProducer") << "RecHit not registered\n";
379  ;
380  }
381  ++counter;
382  }
383 }
384 
385 // Fills the collections.
387 {
388  loadPCaloHits(iEvent);
389  noisify();
390  hbheHits.reserve(firedCells_.size());
391  if(doDigis_)
392  {
393  hbheDigis.reserve(firedCells_.size());
394  }
395  static HcalQIESample zeroSample(0,0,0,0);
396  unsigned nhits=firedCells_.size();
397  // HB and HE
398 
399  for(unsigned ihit=0;ihit<nhits;++ihit)
400  {
401  unsigned cellhashedindex=firedCells_[ihit];
402  const HcalDetId& detid = theDetIds_[cellhashedindex];
403  unsigned subdet=(detid.subdet()==HcalBarrel) ? 0: 1;
404 
405  // Check if it is above the threshold
406  if(hcalRecHits_[cellhashedindex]<threshold_[subdet]) continue;
407  float energy=hcalRecHits_[cellhashedindex];
408  // apply RespCorr only to the RecHit
409  energy *= myRespCorr->getValues(theDetIds_[cellhashedindex])->getValue();
410  // poor man saturation
411  if(energy>sat_[cellhashedindex])
412  {
413  // std::cout << " Saturation " << energy << " " << sat_[cellhashedindex] << " HB " << std::endl;
414  energy=sat_[cellhashedindex];
415  }
416 
417 
418  hbheHits.push_back(HBHERecHit(detid,energy,0.));
419  if(doDigis_)
420  {
421  HBHEDataFrame myDataFrame(detid);
422  myDataFrame.setSize(2);
423  double nfc=hcalRecHits_[cellhashedindex]/gains_[cellhashedindex]+peds_[cellhashedindex];
424  int nadc=fCtoAdc(nfc);
425  HcalQIESample qie(nadc, 0, 0, 0) ;
426  myDataFrame.setSample(0,qie);
427  myDataFrame.setSample(1,zeroSample);
428  hbheDigis.push_back(myDataFrame);
429  }
430  }
431 }
432 
433 
434 // Fills the collections.
436 {
437  loadPCaloHits(iEvent);
438  noisify();
439  hoHits.reserve(firedCells_.size());
440  if(doDigis_)
441  {
442  hoDigis.reserve(firedCells_.size());
443  }
444  static HcalQIESample zeroSample(0,0,0,0);
445 
446  // HO
447  unsigned nhits=firedCells_.size();
448  for(unsigned ihit=0;ihit<nhits;++ihit)
449  {
450 
451  unsigned cellhashedindex=firedCells_[ihit];
452  // Check if it is above the threshold
453  if(hcalRecHits_[cellhashedindex]<threshold_[0]) continue;
454 
455  const HcalDetId& detid=theDetIds_[cellhashedindex];
456  int ieta = detid.ieta();
457 
458  // Force only Ring#0 to remain
459  if(ieta > 4 || ieta < -4 ) continue;
460 
461  float energy=hcalRecHits_[cellhashedindex];
462  // apply RespCorr
463  energy *= myRespCorr->getValues(theDetIds_[cellhashedindex])->getValue();
464 
465  // poor man saturation
466  if(energy>sat_[cellhashedindex]) energy=sat_[cellhashedindex];
467 
468  hoHits.push_back(HORecHit(detid,energy,0));
469  }
470 }
471 
472 // Fills the collections.
474 {
475  loadPCaloHits(iEvent);
476  noisify();
477  hfHits.reserve(firedCells_.size());
478  if(doDigis_)
479  {
480  hfDigis.reserve(firedCells_.size());
481  }
482  static HcalQIESample zeroSample(0,0,0,0);
483 
484  unsigned nhits=firedCells_.size();
485  for(unsigned ihit=0;ihit<nhits;++ihit)
486  {
487  unsigned cellhashedindex=firedCells_[ihit];
488  // Check if it is above the threshold
489  if(hcalRecHits_[cellhashedindex]<threshold_[0]) continue;
490 
491  float energy=hcalRecHits_[cellhashedindex];
492 
493  // apply RespCorr
494  energy *= myRespCorr->getValues(theDetIds_[cellhashedindex])->getValue();
495 
496  // poor man saturation
497  if(energy>sat_[cellhashedindex]) energy=sat_[cellhashedindex];
498 
499  const HcalDetId & detid=theDetIds_[cellhashedindex];
500  hfHits.push_back(HFRecHit(detid,energy,0.));
501  if(doDigis_)
502  {
503  HFDataFrame myDataFrame(detid);
504  myDataFrame.setSize(1);
505 
506  double nfc= hcalRecHits_[cellhashedindex]/gains_[cellhashedindex]+peds_[cellhashedindex];
507  int nadc=fCtoAdc(nfc/2.6);
508  HcalQIESample qie(nadc, 0, 0, 0) ;
509  myDataFrame.setSample(0,qie);
510  hfDigis.push_back(myDataFrame);
511  }
512  }
513 }
514 
515 
516 // For a fast injection of the noise: the list of cell ids is stored
518 {
520  es.get<CaloGeometryRecord>().get(pG);
525 
527 }
528 
529 // list of the cellids for a given subdetector
530 unsigned HcalRecHitsMaker::createVectorOfSubdetectorCells(const CaloGeometry& cg,int subdetn,std::vector<int>& cellsvec )
531 {
532 
534  const std::vector<DetId>& ids=geom->getValidDetIds(DetId::Hcal,subdetn);
535 
536  for (std::vector<DetId>::const_iterator i=ids.begin(); i!=ids.end(); i++)
537  {
538  HcalDetId myDetId(*i);
539  // std::cout << myDetId << myHcalSimParameterMap_->simParameters(myDetId).simHitToPhotoelectrons() << std::endl;;
540  // std::cout << " hi " << hi << " " theDetIds_.size() << std::endl;
541  unsigned hi=myDetId.hashed_index();
542  theDetIds_[hi]=myDetId;
543  // std::cout << myDetId << " " << hi << std::endl;
544  cellsvec.push_back(hi);
545 
546  if(hi>maxIndex_)
547  maxIndex_=hi;
548  }
549  return cellsvec.size();
550 }
551 
552 // Takes a hit (from a PSimHit) and fills a map
553 void HcalRecHitsMaker::Fill(int id, float energy, std::vector<int>& theHits,float noise)
554 {
555  if(doMiscalib_)
556  energy*=miscalib_[id];
557 
558  if(noiseFromDb_)
559  noise=noisesigma_[id];
560 
561  // Check if the RecHit exists
562  if(hcalRecHits_[id]>0.)
563  hcalRecHits_[id]+=energy;
564  else
565  {
566  // the noise is injected only the first time
567  hcalRecHits_[id]=energy + random_->gaussShoot(0.,noise);
568  theHits.push_back(id);
569  }
570 }
571 
573 {
574  unsigned total=0;
575  switch(det_)
576  {
577  case 4:
578  {
579  // do the HB
580  if(noise_[0] != 0.) {
582  }
583  // do the HE
584  if(noise_[1] != 0.) {
586  }
587  }
588  break;
589  case 5:
590  {
591  // do the HO
592  if(noise_[0] != 0.) {
594  }
595  }
596  break;
597  case 6:
598  {
599  // do the HF
600  if(noise_[0] != 0.) {
602  }
603  }
604  break;
605  default:
606  break;
607  }
608  edm::LogInfo("CaloRecHitsProducer") << "CaloRecHitsProducer : added noise in "<< total << " HCAL cells " << std::endl;
609 }
610 
611 unsigned HcalRecHitsMaker::noisifySubdet(std::vector<float>& theMap, std::vector<int>& theHits, const std::vector<int>& thecells, unsigned ncells, double hcalHotFraction,const GaussianTail *myGT,double sigma,double threshold)
612 {
613  // If the fraction of "hot " is small, use an optimized method to inject noise only in noisy cells. The 30% has not been tuned
614  if(!noiseFromDb_ && hcalHotFraction==0.) return 0;
615  if(hcalHotFraction<0.3 && !noiseFromDb_)
616  {
617  double mean = (double)(ncells-theHits.size())*hcalHotFraction;
618  unsigned nhcal = random_->poissonShoot(mean);
619 
620  unsigned ncell=0;
621  unsigned cellindex=0;
622  uint32_t cellhashedindex=0;
623 
624  while(ncell < nhcal)
625  {
626  cellindex = (unsigned)(random_->flatShoot()*ncells);
627  cellhashedindex = thecells[cellindex];
628 
629  if(hcalRecHits_[cellhashedindex]==0.) // new cell
630  {
631  hcalRecHits_[cellhashedindex]=myGT->shoot();
632  theHits.push_back(cellhashedindex);
633  ++ncell;
634  }
635  }
636  return ncell;
637  }
638  else // otherwise, inject noise everywhere
639  {
640  uint32_t cellhashedindex=0;
641  unsigned nhcal=thecells.size();
642 
643 
644  for(unsigned ncell=0;ncell<nhcal;++ncell)
645  {
646  cellhashedindex = thecells[ncell];
647  if(hcalRecHits_[cellhashedindex]==0.) // new cell
648  {
649 
650  sigma=noisesigma_[cellhashedindex];
651 
652  double noise =random_->gaussShoot(0.,sigma);
653  if(noise>threshold)
654  {
655  hcalRecHits_[cellhashedindex]=noise;
656  theHits.push_back(cellhashedindex);
657  }
658  }
659  }
660  return nhcal;
661  }
662  return 0;
663 }
664 
666 {
668 }
669 
670 void HcalRecHitsMaker::cleanSubDet(std::vector<float>& hits,std::vector<int>& cells)
671 {
672  unsigned size=cells.size();
673  // Reset the energies
674  for(unsigned ic=0;ic<size;++ic)
675  {
676  hits[cells[ic]] = 0.;
677  }
678  // Clear the list of fired cells
679  cells.clear();
680 }
681 
682 // fC to ADC conversion
683 int HcalRecHitsMaker::fCtoAdc(double fc) const
684 {
685  if(fc<0.) return 0;
686  if(fc>9985.) return 127;
687  return fctoadc_[(unsigned)floor(fc)];
688 }
689 
690 double HcalRecHitsMaker::noiseInfCfromDB(const HcalDbService * conditions,const HcalDetId & detId)
691 {
692  // method from Salavat
693  // fetch pedestal widths (noise sigmas for all 4 CapID)
694  const HcalPedestalWidth* pedWidth =
695  conditions-> getPedestalWidth(detId);
696  double ssqq_1 = pedWidth->getSigma(0,0);
697  double ssqq_2 = pedWidth->getSigma(1,1);
698  double ssqq_3 = pedWidth->getSigma(2,2);
699  double ssqq_4 = pedWidth->getSigma(3,3);
700 
701  // correction factors (hb,he,ho,hf)
702  static float corrfac[4]={1.39,1.32,1.17,3.76};
703 
704  int sub = detId.subdet();
705 
706  // HO: only Ring#0 matters
707  int ieta = detId.ieta();
708  if (sub == 3 && abs (ieta) > 4) return 0.;
709 
710  // effective RecHits (for this particular detId) noise calculation :
711 
712  double sig_sq_mean = 0.25 * ( ssqq_1 + ssqq_2 + ssqq_3 + ssqq_4);
713 
714  // f - external parameter (currently set to 0.5 in the FullSim) !!!
715  double f=0.5;
716 
717  double term = sqrt (1. + sqrt(1. - 2.*f*f));
718  double alpha = sqrt (0.5) * term;
719  double beta = sqrt (0.5) * f / term;
720 
721  // double RMS1 = sqrt(sig_sq_mean) * sqrt (2.*beta*beta + alpha*alpha) ;
722  double RMS4 = sqrt(sig_sq_mean) * sqrt (2.*beta*beta + 2.*(alpha-beta)*(alpha-beta) + 2.*(alpha-2.*beta)*(alpha-2.*beta)) ;
723 
724  double noise_rms_fC;
725 
726  // if(sub == 4) noise_rms_fC = RMS1;
727  // else noise_rms_fC = RMS4;
728  noise_rms_fC = RMS4;
729 
730  noise_rms_fC *= corrfac[sub-1];
731 
732  // to convert from above fC to GeV - multiply by gain (GeV/fC)
733  // const HcalGain* gain = conditions->getGain(detId);
734  // double noise_rms_GeV = noise_rms_fC * gain->getValue(0); // Noise RMS (GeV) for detId channel
735  return noise_rms_fC;
736 }
const double beta
T getParameter(std::string const &) const
const CaloSubdetectorGeometry * getSubdetectorGeometry(const DetId &id) const
access the subdetector geometry for the given subdetector directly
Definition: CaloGeometry.cc:43
const std::map< uint32_t, float > & get()
int i
Definition: DBlmapReader.cc:9
static std::vector< float > peds_
float alpha
Definition: AMPTWrapper.h:95
HcalRecHitsMaker(edm::ParameterSet const &p, int, const RandomEngine *random)
HcalSubdetector subdet() const
get the subdetector
Definition: HcalDetId.h:32
static std::vector< float > sat_
std::vector< GaussianTail * > myGaussianTailGenerators_
float getSigma(int fCapId1, int fCapId2) const
get correlation element for capId1/2 = 0..3
int fCtoAdc(double fc) const
std::vector< double > hcalHotFraction_
unsigned createVectorOfSubdetectorCells(const CaloGeometry &, int subdetn, std::vector< int > &)
double noiseInfCfromDB(const HcalDbService *conditions, const HcalDetId &detId)
edm::InputTag inputCol_
std::vector< double > noise_
bool parseXMLMiscalibFile(std::string configFile)
std::vector< float > hcalRecHits_
static std::vector< float > TPGFactor_
void push_back(T const &t)
#define abs(x)
Definition: mlp_lapack.h:159
const RandomEngine * random_
static std::vector< int > hfhi_
float getValue(int fCapId) const
get value for capId = 0..3
Definition: HcalGain.h:20
unsigned createVectorsOfCells(const edm::EventSetup &es)
void setSize(int size)
double gaussShoot(double mean=0.0, double sigma=1.0) const
Definition: RandomEngine.h:37
virtual const std::vector< DetId > & getValidDetIds(DetId::Detector det=DetId::Detector(0), int subdet=0) const
Get a list of valid detector ids (for the given subdetector)
static std::vector< float > noisesigma_
int iEvent
Definition: GenABIO.cc:243
std::vector< double > threshold_
static std::vector< int > hohi_
static std::vector< int > fctoadc_
void setSample(int i, const HcalQIESample &sam)
Definition: HBHEDataFrame.h:50
T sqrt(T t)
Definition: SSEVec.h:28
unsigned int poissonShoot(double mean) const
Definition: RandomEngine.h:44
int ieta() const
get the cell ieta
Definition: HcalDetId.h:38
void loadPCaloHits(const edm::Event &iEvent)
static bool initialized_
float getValue(int fCapId) const
get value for capId = 0..3
Definition: HcalPedestal.h:19
double f[11][100]
std::string hcalfileinpath_
static std::vector< HcalDetId > theDetIds_
void setSample(int i, const HcalQIESample &sam)
Definition: HFDataFrame.h:50
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:355
const T & get() const
Definition: EventSetup.h:55
T const * product() const
Definition: ESHandle.h:62
const HcalGain * getGain(const HcalGenericDetId &fId) const
T const * product() const
Definition: Handle.h:74
double flatShoot(double xmin=0.0, double xmax=1.0) const
Definition: RandomEngine.h:30
static unsigned maxIndex_
void setSize(int size)
Definition: HFDataFrame.cc:17
const HcalRespCorrs * myRespCorr
double shoot() const
Definition: GaussianTail.cc:20
static std::vector< float > gains_
float getValue() const
Definition: HcalRespCorr.h:18
unsigned noisifySubdet(std::vector< float > &theMap, std::vector< int > &theHits, const std::vector< int > &thecells, unsigned ncells, double hcalHotFraction_, const GaussianTail *, double sigma, double threshold)
void Fill(int id, float energy, std::vector< int > &myHits, float noise)
void loadHcalRecHits(edm::Event &iEvent, HBHERecHitCollection &hbheHits, HBHEDigiCollection &hbheDigis)
void reserve(size_type n)
void cleanSubDet(std::vector< float > &hits, std::vector< int > &cells)
std::vector< int > firedCells_
static std::vector< int > hbhi_
tuple cout
Definition: gather_cfg.py:41
std::string fullPath() const
Definition: FileInPath.cc:170
Definition: DDAxes.h:10
static std::vector< float > miscalib_
const Item * getValues(DetId fId) const
const HcalPedestal * getPedestal(const HcalGenericDetId &fId) const
static std::vector< int > hehi_
void init(const edm::EventSetup &es, bool dodigis, bool domiscalib)
tuple size
Write out results.
int hashed_index() const
Definition: HcalDetId.cc:119