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  float energy=hcalRecHits_[cellhashedindex];
456  // apply RespCorr
457  energy *= myRespCorr->getValues(theDetIds_[cellhashedindex])->getValue();
458 
459  // poor man saturation
460  if(energy>sat_[cellhashedindex]) energy=sat_[cellhashedindex];
461 
462  const HcalDetId& detid=theDetIds_[cellhashedindex];
463 
464  hoHits.push_back(HORecHit(detid,energy,0));
465  }
466 }
467 
468 // Fills the collections.
470 {
471  loadPCaloHits(iEvent);
472  noisify();
473  hfHits.reserve(firedCells_.size());
474  if(doDigis_)
475  {
476  hfDigis.reserve(firedCells_.size());
477  }
478  static HcalQIESample zeroSample(0,0,0,0);
479 
480  unsigned nhits=firedCells_.size();
481  for(unsigned ihit=0;ihit<nhits;++ihit)
482  {
483  unsigned cellhashedindex=firedCells_[ihit];
484  // Check if it is above the threshold
485  if(hcalRecHits_[cellhashedindex]<threshold_[0]) continue;
486 
487  float energy=hcalRecHits_[cellhashedindex];
488 
489  // apply RespCorr
490  energy *= myRespCorr->getValues(theDetIds_[cellhashedindex])->getValue();
491 
492  // poor man saturation
493  if(energy>sat_[cellhashedindex]) energy=sat_[cellhashedindex];
494 
495  const HcalDetId & detid=theDetIds_[cellhashedindex];
496  hfHits.push_back(HFRecHit(detid,energy,0.));
497  if(doDigis_)
498  {
499  HFDataFrame myDataFrame(detid);
500  myDataFrame.setSize(1);
501 
502  double nfc= hcalRecHits_[cellhashedindex]/gains_[cellhashedindex]+peds_[cellhashedindex];
503  int nadc=fCtoAdc(nfc/2.6);
504  HcalQIESample qie(nadc, 0, 0, 0) ;
505  myDataFrame.setSample(0,qie);
506  hfDigis.push_back(myDataFrame);
507  }
508  }
509 }
510 
511 
512 // For a fast injection of the noise: the list of cell ids is stored
514 {
516  es.get<CaloGeometryRecord>().get(pG);
521 
523 }
524 
525 // list of the cellids for a given subdetector
526 unsigned HcalRecHitsMaker::createVectorOfSubdetectorCells(const CaloGeometry& cg,int subdetn,std::vector<int>& cellsvec )
527 {
528 
530  const std::vector<DetId>& ids=geom->getValidDetIds(DetId::Hcal,subdetn);
531 
532  for (std::vector<DetId>::const_iterator i=ids.begin(); i!=ids.end(); i++)
533  {
534  HcalDetId myDetId(*i);
535  // std::cout << myDetId << myHcalSimParameterMap_->simParameters(myDetId).simHitToPhotoelectrons() << std::endl;;
536  // std::cout << " hi " << hi << " " theDetIds_.size() << std::endl;
537  unsigned hi=myDetId.hashed_index();
538  theDetIds_[hi]=myDetId;
539  // std::cout << myDetId << " " << hi << std::endl;
540  cellsvec.push_back(hi);
541 
542  if(hi>maxIndex_)
543  maxIndex_=hi;
544  }
545  return cellsvec.size();
546 }
547 
548 // Takes a hit (from a PSimHit) and fills a map
549 void HcalRecHitsMaker::Fill(int id, float energy, std::vector<int>& theHits,float noise)
550 {
551  if(doMiscalib_)
552  energy*=miscalib_[id];
553 
554  if(noiseFromDb_)
555  noise=noisesigma_[id];
556 
557  // Check if the RecHit exists
558  if(hcalRecHits_[id]>0.)
560  else
561  {
562  // the noise is injected only the first time
563  hcalRecHits_[id]=energy + random_->gaussShoot(0.,noise);
564  theHits.push_back(id);
565  }
566 }
567 
569 {
570  unsigned total=0;
571  switch(det_)
572  {
573  case 4:
574  {
575  // do the HB
576  if(noise_[0] != 0.) {
578  }
579  // do the HE
580  if(noise_[1] != 0.) {
582  }
583  }
584  break;
585  case 5:
586  {
587  // do the HO
588  if(noise_[0] != 0.) {
590  }
591  }
592  break;
593  case 6:
594  {
595  // do the HF
596  if(noise_[0] != 0.) {
598  }
599  }
600  break;
601  default:
602  break;
603  }
604  edm::LogInfo("CaloRecHitsProducer") << "CaloRecHitsProducer : added noise in "<< total << " HCAL cells " << std::endl;
605 }
606 
607 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)
608 {
609  // If the fraction of "hot " is small, use an optimized method to inject noise only in noisy cells. The 30% has not been tuned
610  if(!noiseFromDb_ && hcalHotFraction==0.) return 0;
611  if(hcalHotFraction<0.3 && !noiseFromDb_)
612  {
613  double mean = (double)(ncells-theHits.size())*hcalHotFraction;
614  unsigned nhcal = random_->poissonShoot(mean);
615 
616  unsigned ncell=0;
617  unsigned cellindex=0;
618  uint32_t cellhashedindex=0;
619 
620  while(ncell < nhcal)
621  {
622  cellindex = (unsigned)(random_->flatShoot()*ncells);
623  cellhashedindex = thecells[cellindex];
624  if(hcalRecHits_[cellhashedindex]==0.) // new cell
625  {
626  hcalRecHits_[cellhashedindex]=myGT->shoot();
627  theHits.push_back(cellhashedindex);
628  ++ncell;
629  }
630  }
631  return ncell;
632  }
633  else // otherwise, inject noise everywhere
634  {
635  uint32_t cellhashedindex=0;
636  unsigned nhcal=thecells.size();
637 
638 
639  for(unsigned ncell=0;ncell<nhcal;++ncell)
640  {
641  cellhashedindex = thecells[ncell];
642  if(hcalRecHits_[cellhashedindex]==0.) // new cell
643  {
644 
645  sigma=noisesigma_[cellhashedindex];
646 
647  double noise =random_->gaussShoot(0.,sigma);
648  if(noise>threshold)
649  {
650  hcalRecHits_[cellhashedindex]=noise;
651  theHits.push_back(cellhashedindex);
652  }
653  }
654  }
655  return nhcal;
656  }
657  return 0;
658 }
659 
661 {
663 }
664 
665 void HcalRecHitsMaker::cleanSubDet(std::vector<float>& hits,std::vector<int>& cells)
666 {
667  unsigned size=cells.size();
668  // Reset the energies
669  for(unsigned ic=0;ic<size;++ic)
670  {
671  hits[cells[ic]] = 0.;
672  }
673  // Clear the list of fired cells
674  cells.clear();
675 }
676 
677 // fC to ADC conversion
678 int HcalRecHitsMaker::fCtoAdc(double fc) const
679 {
680  if(fc<0.) return 0;
681  if(fc>9985.) return 127;
682  return fctoadc_[(unsigned)floor(fc)];
683 }
684 
685 double HcalRecHitsMaker::noiseInfCfromDB(const HcalDbService * conditions,const HcalDetId & detId)
686 {
687  // method from Salavat
688  // fetch pedestal widths (noise sigmas for all 4 CapID)
689  const HcalPedestalWidth* pedWidth =
690  conditions-> getPedestalWidth(detId);
691  double ssqq_1 = pedWidth->getSigma(0,0);
692  double ssqq_2 = pedWidth->getSigma(1,1);
693  double ssqq_3 = pedWidth->getSigma(2,2);
694  double ssqq_4 = pedWidth->getSigma(3,3);
695 
696  // correction factors (hb,he,ho,hf)
697  static float corrfac[4]={1.39,1.32,1.53,3.76};
698 
699  int sub = detId.subdet();
700 
701  // effective RecHits (for this particular detId) noise calculation :
702 
703  double sig_sq_mean = 0.25 * ( ssqq_1 + ssqq_2 + ssqq_3 + ssqq_4);
704 
705  // f - external parameter (currently set to 0.5 in the FullSim) !!!
706  double f=0.5;
707 
708  double term = sqrt (1. + sqrt(1. - 2.*f*f));
709  double alpha = sqrt (0.5) * term;
710  double beta = sqrt (0.5) * f / term;
711 
712  // double RMS1 = sqrt(sig_sq_mean) * sqrt (2.*beta*beta + alpha*alpha) ;
713  double RMS4 = sqrt(sig_sq_mean) * sqrt (2.*beta*beta + 2.*(alpha-beta)*(alpha-beta) + 2.*(alpha-2.*beta)*(alpha-2.*beta)) ;
714 
715  double noise_rms_fC;
716 
717  // if(sub == 4) noise_rms_fC = RMS1;
718  // else noise_rms_fC = RMS4;
719  noise_rms_fC = RMS4;
720 
721  noise_rms_fC *= corrfac[sub-1];
722 
723  // to convert from above fC to GeV - multiply by gain (GeV/fC)
724  // const HcalGain* gain = conditions->getGain(detId);
725  // double noise_rms_GeV = noise_rms_fC * gain->getValue(0); // Noise RMS (GeV) for detId channel
726  return noise_rms_fC;
727 }
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)
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
float threshold
Definition: crabWrap.py:319
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:359
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
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