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)
456 int ieta = detid.
ieta();
459 if(ieta > 4 || ieta < -4 )
continue;
466 if(energy>
sat_[cellhashedindex]) energy=
sat_[cellhashedindex];
485 for(
unsigned ihit=0;ihit<nhits;++ihit)
497 if(energy>
sat_[cellhashedindex]) energy=
sat_[cellhashedindex];
536 for (std::vector<DetId>::const_iterator
i=ids.begin();
i!=ids.end();
i++)
544 cellsvec.push_back(hi);
549 return cellsvec.size();
568 theHits.push_back(
id);
608 edm::LogInfo(
"CaloRecHitsProducer") <<
"CaloRecHitsProducer : added noise in "<< total <<
" HCAL cells " << std::endl;
617 double mean = (double)(ncells-theHits.size())*hcalHotFraction;
621 unsigned cellindex=0;
622 uint32_t cellhashedindex=0;
627 cellhashedindex = thecells[cellindex];
632 theHits.push_back(cellhashedindex);
640 uint32_t cellhashedindex=0;
641 unsigned nhcal=thecells.size();
644 for(
unsigned ncell=0;ncell<nhcal;++ncell)
646 cellhashedindex = thecells[ncell];
656 theHits.push_back(cellhashedindex);
672 unsigned size=cells.size();
674 for(
unsigned ic=0;ic<
size;++ic)
676 hits[cells[ic]] = 0.;
686 if(fc>9985.)
return 127;
687 return fctoadc_[(unsigned)floor(fc)];
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);
702 static float corrfac[4]={1.39,1.32,1.17,3.76};
707 int ieta = detId.
ieta();
708 if (sub == 3 &&
abs (ieta) > 4)
return 0.;
712 double sig_sq_mean = 0.25 * ( ssqq_1 + ssqq_2 + ssqq_3 + ssqq_4);
717 double term =
sqrt (1. +
sqrt(1. - 2.*f*f));
719 double beta =
sqrt (0.5) * f / term;
722 double RMS4 =
sqrt(sig_sq_mean) *
sqrt (2.*beta*beta + 2.*(alpha-beta)*(alpha-beta) + 2.*(alpha-2.*beta)*(alpha-2.*beta)) ;
730 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.