15 #include "CLHEP/GenericFunctions/Erf.hh"
79 if(threshold_.size()!=noise_.size())
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;
97 for(
unsigned inoise=0;inoise<
nnoise_;++inoise)
104 else if(noise_[inoise]==-1) {
153 edm::LogInfo(
"CaloRecHitsProducer") <<
"Total number of cells in HCAL " << ncells << std::endl;
174 std::string hcalfile=hcalfiletmp.
fullPath();
180 std::map<uint32_t,float>::const_iterator it=mapHcal.
get().begin();
181 std::map<uint32_t,float>::const_iterator itend=mapHcal.
get().end();
185 float icalconst=it->second;
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();
202 for(
unsigned ibin=0;ibin<
size;++ibin)
205 myGraf->GetPoint(ibin,x,y);
206 int index=(int)floor(x);
207 if(index<0||index>=10000)
continue;
211 for(
unsigned ivec=p_index;ivec<(unsigned)index;++ivec)
219 edm::FileInPath myTPGFilePath(
"CalibCalorimetry/HcalTPGAlgos/data/RecHit-TPG-calib.dat");
226 for(
unsigned i=0;
i<86;++
i)
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;
242 for(
unsigned ig=0;ig<4;++ig)
250 float tpgfactor=
TPGFactor_[(ieta>0)?ieta+43:-ieta];
251 mgain*=2500./0.92*tpgfactor ;
266 for(
unsigned ig=0;ig<4;++ig)
272 float tpgfactor=
TPGFactor_[(ieta>0)?ieta+44:-ieta+1];
273 mgain*=2500./0.92*tpgfactor;
287 for(
unsigned ig=0;ig<4;++ig)
293 float tpgfactor=
TPGFactor_[(ieta>0)?ieta+43:-ieta];
294 mgain*=2500./0.92*tpgfactor;
298 peds_[hohi_[ic]]=ped;
308 for(
unsigned ig=0;ig<4;++ig)
322 peds_[hfhi_[ic]]=ped;
350 int hashedindex=
detid.hashed_index();
351 switch(
detid.subdet())
399 for(
unsigned ihit=0;ihit<nhits;++ihit)
411 if(energy>
sat_[cellhashedindex])
414 energy=
sat_[cellhashedindex];
448 for(
unsigned ihit=0;ihit<nhits;++ihit)
460 if(energy>
sat_[cellhashedindex]) energy=
sat_[cellhashedindex];
481 for(
unsigned ihit=0;ihit<nhits;++ihit)
493 if(energy>
sat_[cellhashedindex]) energy=
sat_[cellhashedindex];
532 for (std::vector<DetId>::const_iterator
i=ids.begin();
i!=ids.end();
i++)
540 cellsvec.push_back(hi);
545 return cellsvec.size();
564 theHits.push_back(
id);
604 edm::LogInfo(
"CaloRecHitsProducer") <<
"CaloRecHitsProducer : added noise in "<< total <<
" HCAL cells " << std::endl;
613 double mean = (double)(ncells-theHits.size())*hcalHotFraction;
617 unsigned cellindex=0;
618 uint32_t cellhashedindex=0;
623 cellhashedindex = thecells[cellindex];
627 theHits.push_back(cellhashedindex);
635 uint32_t cellhashedindex=0;
636 unsigned nhcal=thecells.size();
639 for(
unsigned ncell=0;ncell<nhcal;++ncell)
641 cellhashedindex = thecells[ncell];
651 theHits.push_back(cellhashedindex);
667 unsigned size=cells.size();
669 for(
unsigned ic=0;ic<
size;++ic)
671 hits[cells[ic]] = 0.;
681 if(fc>9985.)
return 127;
682 return fctoadc_[(unsigned)floor(fc)];
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);
697 static float corrfac[4]={1.39,1.32,1.53,3.76};
703 double sig_sq_mean = 0.25 * ( ssqq_1 + ssqq_2 + ssqq_3 + ssqq_4);
708 double term =
sqrt (1. +
sqrt(1. - 2.*f*f));
710 double beta =
sqrt (0.5) * f / term;
713 double RMS4 =
sqrt(sig_sq_mean) *
sqrt (2.*beta*beta + 2.*(alpha-beta)*(alpha-beta) + 2.*(alpha-2.*beta)*(alpha-2.*beta)) ;
721 noise_rms_fC *= corrfac[sub-1];
T getParameter(std::string const &) const
const CaloSubdetectorGeometry * getSubdetectorGeometry(const DetId &id) const
access the subdetector geometry for the given subdetector directly
const std::map< uint32_t, float > & get()
static std::vector< float > peds_
HcalRecHitsMaker(edm::ParameterSet const &p, int, const RandomEngine *random)
HcalSubdetector subdet() const
get the subdetector
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)
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
unsigned createVectorsOfCells(const edm::EventSetup &es)
double gaussShoot(double mean=0.0, double sigma=1.0) const
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_
std::vector< double > threshold_
static std::vector< int > hohi_
static std::vector< int > fctoadc_
void setSample(int i, const HcalQIESample &sam)
unsigned int poissonShoot(double mean) const
int ieta() const
get the cell ieta
void loadPCaloHits(const edm::Event &iEvent)
float getValue(int fCapId) const
get value for capId = 0..3
std::string hcalfileinpath_
static std::vector< HcalDetId > theDetIds_
void setSample(int i, const HcalQIESample &sam)
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
T const * product() const
const HcalGain * getGain(const HcalGenericDetId &fId) const
T const * product() const
double flatShoot(double xmin=0.0, double xmax=1.0) const
static unsigned maxIndex_
const HcalRespCorrs * myRespCorr
static std::vector< float > gains_
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_
std::string fullPath() const
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.