#include <HcalRecHitsMaker.h>

Public Member Functions
	HcalRecHitsMaker (edm::ParameterSet const &p, int)

void	init (const edm::EventSetup &es, bool dodigis, bool domiscalib)

void	loadHcalRecHits (edm::Event &iEvent, const HcalTopology &, HBHERecHitCollection &hbheHits, HBHEDigiCollection &hbheDigis, RandomEngineAndDistribution const *)

void	loadHcalRecHits (edm::Event &iEvent, const HcalTopology &, HORecHitCollection &ho, HODigiCollection &hoDigis, RandomEngineAndDistribution const *)

void	loadHcalRecHits (edm::Event &iEvent, const HcalTopology &, HFRecHitCollection &hfHits, HFDigiCollection &hfDigis, RandomEngineAndDistribution const *)

	~HcalRecHitsMaker ()

Private Member Functions
void	clean ()

void	cleanSubDet (std::vector< float > &hits, std::vector< int > &cells)

unsigned	createVectorOfSubdetectorCells (const CaloGeometry &, const HcalTopology &, int subdetn, std::vector< int > &)

unsigned	createVectorsOfCells (const edm::EventSetup &es)

int	fCtoAdc (double fc) const

void	Fill (int id, float energy, std::vector< int > &myHits, float noise, float correctionfactor, RandomEngineAndDistribution const *)

double	fractionOOT (int time_slice)

void	loadPCaloHits (const edm::Event &iEvent, const HcalTopology &, RandomEngineAndDistribution const *)

double	noiseInfCfromDB (const HcalDbService *conditions, const HcalDetId &detId)

void	noisify (RandomEngineAndDistribution const *)

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, double correctionfactor, RandomEngineAndDistribution const )

Private Attributes
std::vector< double >	corrfac_

unsigned	det_

bool	doDigis_

bool	doMiscalib_

bool	doSaturation_

std::vector< int >	fctoadc_

std::vector< int >	firedCells_

std::vector< float >	gains_

std::vector< int >	hbhi_

std::string	hcalfileinpath_

std::vector< double >	hcalHotFraction_

std::vector< float >	hcalRecHits_

std::vector< int >	hehi_

std::vector< int >	hfhi_

std::vector< int >	hohi_

bool	initialized_

bool	initializedHB_

bool	initializedHE_

bool	initializedHF_

bool	initializedHO_

edm::InputTag	inputCol_

unsigned	maxIndex_

std::vector< float >	miscalib_

std::vector< GaussianTail * >	myGaussianTailGenerators_

const HcalRespCorrs *	myRespCorr

unsigned	nhbcells_

unsigned	nhecells_

unsigned	nhfcells_

unsigned	nhocells_

unsigned	nnoise_

std::vector< double >	noise_

std::vector< bool >	noiseFromDb_

std::vector< float >	noisesigma_

std::vector< float >	peds_

double	refactor_

double	refactor_mean_

std::vector< float >	sat_

std::vector< HcalDetId >	theDetIds_

std::vector< double >	threshold_

std::vector< float >	TPGFactor_

Detailed Description

Definition at line 24 of file HcalRecHitsMaker.h.

Constructor & Destructor Documentation

HcalRecHitsMaker::HcalRecHitsMaker	(	edm::ParameterSet const &	p,
		int	det
	)

Definition at line 34 of file HcalRecHitsMaker.cc.

References corrfac_, det_, doSaturation_, edm::ParameterSet::getParameter(), hbhi_, hcalfileinpath_, hcalHotFraction_, hehi_, hfhi_, hohi_, initialized_, initializedHB_, initializedHE_, initializedHF_, initializedHO_, inputCol_, maxIndex_, myGaussianTailGenerators_, nhbcells_, nhecells_, nhfcells_, nhocells_, nnoise_, noise_, noiseFromDb_, refactor_, refactor_mean_, mathSSE::sqrt(), AlCaHLTBitMon_QueryRunRegistry::string, and threshold_.

   :
   det_(det),
   doDigis_(false)
 {
   edm::ParameterSet RecHitsParameters=p.getParameter<edm::ParameterSet>("HCAL");
   noise_ = RecHitsParameters.getParameter<std::vector<double> >("Noise");
   corrfac_ = RecHitsParameters.getParameter<std::vector<double> >("NoiseCorrectionFactor");
   threshold_ = RecHitsParameters.getParameter<std::vector<double> >("Threshold");
   doSaturation_ = RecHitsParameters.getParameter<bool>("EnableSaturation");
     
   refactor_ = RecHitsParameters.getParameter<double> ("Refactor");
   refactor_mean_ = RecHitsParameters.getParameter<double> ("Refactor_mean");
   hcalfileinpath_= RecHitsParameters.getParameter<std::string> ("fileNameHcal");  
   inputCol_=RecHitsParameters.getParameter<edm::InputTag>("MixedSimHits");
   nhbcells_=nhecells_=nhocells_=nhfcells_=0;
 
   initialized_ = false;
   initializedHB_ = false;
   initializedHE_ = false;
   initializedHO_ = false;
   initializedHF_ = false;
   maxIndex_  = 0 ;
 
   if(det_==4)
     {
       hbhi_.reserve(2600);
       hehi_.reserve(2600);
     }
   else if (det_==5)
     hohi_.reserve(2200);
   else if (det_==6)
     hfhi_.reserve(1800);
 
   nnoise_=noise_.size();
   if(threshold_.size()!=noise_.size())
     {
       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;
       noise_.resize(nnoise_,0);
     }
 
   //  edm::ParameterSet hcalparam = p2.getParameter<edm::ParameterSet>("hcalSimParam"); 
   //  myHcalSimParameterMap_ = new HcalSimParameterMap(hcalparam);
 
   // Computes the fraction of HCAL above the threshold
   Genfun::Erf myErf;
   noiseFromDb_.resize(nnoise_,false);
   hcalHotFraction_.resize(nnoise_,0.);
   myGaussianTailGenerators_.resize(nnoise_,0);
   if(noise_.size()>0) 
     {
       for(unsigned inoise=0;inoise<nnoise_;++inoise)
     {
       if(noise_[inoise]==0) {
         hcalHotFraction_[inoise]=0.;
         continue;
       } else if(noise_[inoise]==-1) {
         noiseFromDb_[inoise]=true;
         continue;
       } else {
         hcalHotFraction_[inoise] = 0.5-0.5*myErf(threshold_[inoise]/noise_[inoise]/sqrt(2.));
         myGaussianTailGenerators_[inoise]=new GaussianTail(noise_[inoise],threshold_[inoise]);
       }
     }   
     }  
 }

HcalRecHitsMaker::~HcalRecHitsMaker ( )

Definition at line 101 of file HcalRecHitsMaker.cc.

References clean(), hbhi_, hehi_, hfhi_, hohi_, myGaussianTailGenerators_, and theDetIds_.

 {
   clean();  
   if(myGaussianTailGenerators_.size()) 
     {
       for(unsigned igt=0; igt<myGaussianTailGenerators_.size();++igt)
         delete myGaussianTailGenerators_[igt];
     }
   myGaussianTailGenerators_.clear();
   theDetIds_.clear();
   hbhi_.clear();
   hehi_.clear();
   hohi_.clear();
   hfhi_.clear();
     
 }

Member Function Documentation

void HcalRecHitsMaker::clean ( )

private

Definition at line 685 of file HcalRecHitsMaker.cc.

References cleanSubDet(), firedCells_, and hcalRecHits_.

Referenced by loadPCaloHits(), and ~HcalRecHitsMaker().

 {
   cleanSubDet(hcalRecHits_,firedCells_);
 }

void HcalRecHitsMaker::cleanSubDet	(	std::vector< float > &	hits,
		std::vector< int > &	cells
	)

private

Definition at line 690 of file HcalRecHitsMaker.cc.

References findQualityFiles::size.

Referenced by clean().

 {
   unsigned size=cells.size();
   // Reset the energies
   for(unsigned ic=0;ic<size;++ic)
     {
       hits[cells[ic]] = 0.;
     }
   // Clear the list of fired cells 
   cells.clear();
 }

unsigned HcalRecHitsMaker::createVectorOfSubdetectorCells	(	const CaloGeometry &	cg,
		const HcalTopology &	topo,
		int	subdetn,
		std::vector< int > &	cellsvec
	)

private

Definition at line 547 of file HcalRecHitsMaker.cc.

References HcalTopology::detId2denseId(), relativeConstraints::geom, CaloGeometry::getSubdetectorGeometry(), CaloSubdetectorGeometry::getValidDetIds(), DetId::Hcal, i, maxIndex_, and theDetIds_.

Referenced by createVectorsOfCells().

 {
 
   const CaloSubdetectorGeometry* geom=cg.getSubdetectorGeometry(DetId::Hcal,subdetn);  
   const std::vector<DetId>& ids=geom->getValidDetIds(DetId::Hcal,subdetn);  
 
   for (std::vector<DetId>::const_iterator i=ids.begin(); i!=ids.end(); i++) 
     {
       HcalDetId myDetId(*i);
       //      std::cout << myDetId << myHcalSimParameterMap_->simParameters(myDetId).simHitToPhotoelectrons() << std::endl;;
       //      std::cout << " hi " << hi << " " theDetIds_.size() << std::endl;
       unsigned int hi=topo.detId2denseId(myDetId);
       theDetIds_[hi]=myDetId;
       //      std::cout << myDetId << " " << hi <<  std::endl;
       cellsvec.push_back(hi);      
 
       if(hi>maxIndex_)
     maxIndex_=hi;
     }
   return cellsvec.size();
 }

unsigned HcalRecHitsMaker::createVectorsOfCells ( const edm::EventSetup & es )

private

Definition at line 531 of file HcalRecHitsMaker.cc.

References createVectorOfSubdetectorCells(), edm::EventSetup::get(), hbhi_, HcalBarrel, HcalEndcap, HcalForward, HcalOuter, hehi_, hfhi_, hohi_, nhbcells_, nhecells_, nhfcells_, and nhocells_.

Referenced by init().

 {
     edm::ESHandle<CaloGeometry> pG;
     es.get<CaloGeometryRecord>().get(pG);     
     edm::ESHandle<HcalTopology> topo;
     es.get<IdealGeometryRecord>().get( topo );
 
     nhbcells_ = createVectorOfSubdetectorCells(*pG, *topo, HcalBarrel,  hbhi_);    
     nhecells_ = createVectorOfSubdetectorCells(*pG, *topo, HcalEndcap,  hehi_);
     nhocells_ = createVectorOfSubdetectorCells(*pG, *topo, HcalOuter,   hohi_);
     nhfcells_ = createVectorOfSubdetectorCells(*pG, *topo, HcalForward, hfhi_);    
 
     return nhbcells_+nhecells_+nhocells_+nhfcells_;
 }

int HcalRecHitsMaker::fCtoAdc ( double fc ) const

private

Definition at line 703 of file HcalRecHitsMaker.cc.

References fctoadc_.

Referenced by loadHcalRecHits().

 {
   if(fc<0.) return 0;
   if(fc>9985.) return 127;
   return fctoadc_[(unsigned)floor(fc)];
 }

void HcalRecHitsMaker::Fill	(	int	id,
		float	energy,
		std::vector< int > &	myHits,
		float	noise,
		float	correctionfactor,
		RandomEngineAndDistribution const *	random
	)

private

Definition at line 570 of file HcalRecHitsMaker.cc.

References doMiscalib_, relval_parameters_module::energy, RandomEngineAndDistribution::gaussShoot(), hcalRecHits_, miscalib_, and noisesigma_.

Referenced by loadPCaloHits().

 {
   if(doMiscalib_) 
     energy*=miscalib_[id];
 
   if(noise==-1)
     noise=noisesigma_[id]*correctionfactor;
 
   // Check if the RecHit exists
   if(hcalRecHits_[id]>0.)
     hcalRecHits_[id]+=energy;
   else
     {
       // the noise is injected only the first time
       hcalRecHits_[id]=energy + random->gaussShoot(0.,noise);
       theHits.push_back(id);
     }
 }

double HcalRecHitsMaker::fractionOOT ( int time_slice )

private

Definition at line 759 of file HcalRecHitsMaker.cc.

Referenced by loadPCaloHits().

 {
   double SF = 100./88.; // to normalize to in-time signal (88% is in the reco window)
   if (time_slice==-4) {
     return 0.02*SF;
   } else if (time_slice==-3) {
     return 0.06*SF;
   } else if (time_slice==-2) {
     return 0.19*SF;
   } else if (time_slice==-1) {
     return 0.24*SF;
   } else if (time_slice==0) {
     return 1.;
   } else if (time_slice==1) {
     return 0.70*SF;
   } else {
     return 0.;
   }
 
 }

void HcalRecHitsMaker::init	(	const edm::EventSetup &	es,
		bool	dodigis,
		bool	domiscalib
	)

Definition at line 118 of file HcalRecHitsMaker.cc.

References gather_cfg::cout, createVectorsOfCells(), det_, doDigis_, doMiscalib_, doSaturation_, fctoadc_, edm::FileInPath::fullPath(), gains_, edm::EventSetup::get(), CaloMiscalibMapHcal::get(), HcalDbService::getGain(), HcalDbService::getPedestal(), HcalPedestal::getValue(), HcalGain::getValue(), hbhi_, hcalfileinpath_, hcalRecHits_, hehi_, hfhi_, hohi_, i, HcalDetId::ieta(), recoMuon::in, cmsHarvester::index, initialized_, initializedHB_, initializedHE_, initializedHF_, initializedHO_, maxIndex_, miscalib_, myRespCorr, nhbcells_, nhecells_, nhfcells_, nhocells_, noiseFromDb_, noiseInfCfromDB(), noisesigma_, MiscalibReaderFromXML::parseXMLMiscalibFile(), peds_, CaloMiscalibMapHcal::prefillMap(), edm::ESHandle< class >::product(), refactor_, refactor_mean_, sat_, findQualityFiles::size, AlCaHLTBitMon_QueryRunRegistry::string, theDetIds_, TPGFactor_, x, and detailsBasic3DVector::y.

Referenced by CaloRecHitsProducer::beginRun().

 {
   doDigis_=doDigis;
   doMiscalib_=doMiscalib;
 // needed for the noise simulation
   edm::ESHandle<HcalDbService> conditions;
   es.get<HcalDbRecord>().get(conditions);
   const HcalDbService * theDbService=conditions.product();
 
   // get the correction factors
   edm::ESHandle<HcalRespCorrs> rchandle;
   es.get<HcalRespCorrsRcd>().get(rchandle);
   myRespCorr= rchandle.product();
   
   //general initialization
   if(!initialized_) 
     {
       theDetIds_.resize(9201);
       unsigned ncells=createVectorsOfCells(es);
       edm::LogInfo("CaloRecHitsProducer") << "Total number of cells in HCAL " << ncells << std::endl;
       hcalRecHits_.resize(maxIndex_+1,0.);
       edm::LogInfo("CaloRecHitsProducer") << "Largest HCAL hashedindex" << maxIndex_ << std::endl;
 
       peds_.resize(9201);
       gains_.resize(9201);
       sat_.resize(9201);
       noisesigma_.resize(9201);
       
       miscalib_.resize(maxIndex_+1,1.);
       // Read from file ( a la HcalRecHitsRecalib.cc)
       // here read them from xml (particular to HCAL)
       CaloMiscalibMapHcal mapHcal;
   
       edm::ESHandle<HcalTopology> topo;
       es.get<IdealGeometryRecord>().get( topo );
       mapHcal.prefillMap(*topo);
   
       edm::FileInPath hcalfiletmp("CalibCalorimetry/CaloMiscalibTools/data/"+hcalfileinpath_);      
       std::string hcalfile=hcalfiletmp.fullPath();            
       if(doMiscalib_ && !hcalfile.empty()) 
         {
           MiscalibReaderFromXMLHcal hcalreader_(mapHcal);
         
           hcalreader_.parseXMLMiscalibFile(hcalfile);
           //      mapHcal.print();
           std::map<uint32_t,float>::const_iterator it=mapHcal.get().begin();
           std::map<uint32_t,float>::const_iterator itend=mapHcal.get().end();
           for(;it!=itend;++it)
             {
               HcalDetId myDetId(it->first);
               float icalconst=it->second;
               miscalib_[topo->detId2denseId(myDetId)]=refactor_mean_+(icalconst-refactor_mean_)*refactor_;
             }
         }
       
       // Open the histogram for the fC to ADC conversion
       gROOT->cd();
       edm::FileInPath myDataFile("FastSimulation/CaloRecHitsProducer/data/adcvsfc.root");
       TFile * myFile = new TFile(myDataFile.fullPath().c_str(),"READ");
       TGraph * myGraf = (TGraph*)myFile->Get("adcvsfc");
       unsigned size=myGraf->GetN();
       fctoadc_.resize(10000);
       unsigned p_index=0;
       fctoadc_[0]=0;
       int prev_nadc=0;
       int nadc=0;
       for(unsigned ibin=0;ibin<size;++ibin)
         {
           double x,y;
           myGraf->GetPoint(ibin,x,y);
           int index=(int)floor(x);
           if(index<0||index>=10000) continue;
           prev_nadc=nadc;
           nadc=(int)y;
           // Now fills the vector
           for(unsigned ivec=p_index;ivec<(unsigned)index;++ivec)
             {
               fctoadc_[ivec] = prev_nadc;
             }
           p_index = index;
         }
       myFile->Close();
       gROOT->cd();
       edm::FileInPath myTPGFilePath("CalibCalorimetry/HcalTPGAlgos/data/RecHit-TPG-calib.dat");
       TPGFactor_.resize(87,1.2);
       std::ifstream  myTPGFile(myTPGFilePath.fullPath().c_str(),std::ifstream::in);
       if(myTPGFile)
         {
           float gain;
           myTPGFile >> gain;
           for(unsigned i=0;i<86;++i)
             {
               myTPGFile >> TPGFactor_[i] ;
               //std::cout << TPGFactor_[i] << std::endl;
             }
         }
       else
         {
           std::cout << " Unable to open CalibCalorimetry/HcalTPGAlgos/data/RecHit-TPG-calib.dat" << std::endl;
           std::cout <<  " Using a constant 1.2 factor " << std::endl;
         }
           
       initialized_=true; 
     }
     
     if(!initializedHB_ && det_==4){
       nhbcells_ = hbhi_.size();
       //HB
       for(unsigned ic=0;ic<nhbcells_;++ic)
         {
           float gain = theDbService->getGain(theDetIds_[hbhi_[ic]])->getValue(0);
           float mgain=0.;
           for(unsigned ig=0;ig<4;++ig)
             mgain+=theDbService->getGain(theDetIds_[hbhi_[ic]])->getValue(ig);
           if(noiseFromDb_[0])
             noisesigma_[hbhi_[ic]]=noiseInfCfromDB(theDbService,theDetIds_[hbhi_[ic]])*mgain*0.25;
           // std::cout << " NOISEHB " << theDetIds_[hbhi_[ic]].ieta() << " " << noisesigma_[hbhi_[ic]] << "  "<< std::endl;
           // 10000 (ADC scale) / 4. (to compute the mean) / 0.92  ADC/fC
           // *1.25 (only ~80% in 1ts Digi, while saturation applied to 4ts RecHit) 
           int ieta=theDetIds_[hbhi_[ic]].ieta();
           float tpgfactor=TPGFactor_[(ieta>0)?ieta+43:-ieta];
           mgain*=2500./0.92*tpgfactor ;// 10000 (ADC scale) / 4. (to compute the mean)
           sat_[hbhi_[ic]]=(doSaturation_)?mgain:99999.;
           
           peds_[hbhi_[ic]]=theDbService->getPedestal(theDetIds_[hbhi_[ic]])->getValue(0);
         
           gain*=tpgfactor;
           gains_[hbhi_[ic]]=gain;
         }
       
       initializedHB_=true;
     }
     
     if(!initializedHE_ && det_==4){
       nhecells_ = hehi_.size();
      //HE
       for(unsigned ic=0;ic<nhecells_;++ic)
         {
           float gain= theDbService->getGain(theDetIds_[hehi_[ic]])->getValue(0);
           int ieta=theDetIds_[hehi_[ic]].ieta();
           float mgain=0.;
           for(unsigned ig=0;ig<4;++ig)
             mgain+=theDbService->getGain(theDetIds_[hehi_[ic]])->getValue(ig);
           if(noiseFromDb_[1])
             noisesigma_[hehi_[ic]]=noiseInfCfromDB(theDbService,theDetIds_[hehi_[ic]])*mgain*0.25;
           //      std::cout << " NOISEHE " << theDetIds_[hehi_[ic]].ieta() << " " << noisesigma_[hehi_[ic]] << "  "<< std::endl;
           float tpgfactor=TPGFactor_[(ieta>0)?ieta+44:-ieta+1];
           mgain*=2500./0.92*tpgfactor;
           sat_[hehi_[ic]]=(doSaturation_)?mgain:99999.;
           
           gain*=tpgfactor;
           peds_[hehi_[ic]]=theDbService->getPedestal(theDetIds_[hehi_[ic]])->getValue(0);
           gains_[hehi_[ic]]=gain;
         }
 
       initializedHE_=true;
     }
     
     if(!initializedHO_ && det_==5){
       nhocells_ = hohi_.size();
      //HO
       for(unsigned ic=0;ic<nhocells_;++ic)
         {
           float ped=theDbService->getPedestal(theDetIds_[hohi_[ic]])->getValue(0);
           float gain=theDbService->getGain(theDetIds_[hohi_[ic]])->getValue(0);
           float mgain=0.;
           for(unsigned ig=0;ig<4;++ig)
             mgain+=theDbService->getGain(theDetIds_[hohi_[ic]])->getValue(ig);
           if(noiseFromDb_[0])
             noisesigma_[hohi_[ic]]=noiseInfCfromDB(theDbService,theDetIds_[hohi_[ic]])*mgain*0.25;
           //      std::cout << " NOISEHO " << theDetIds_[hohi_[ic]].ieta() << " " << noisesigma_[hohi_[ic]] << "  "<< std::endl;
           int ieta=HcalDetId(hohi_[ic]).ieta();
           float tpgfactor=TPGFactor_[(ieta>0)?ieta+43:-ieta];
           mgain*=2500./0.92*tpgfactor;
           sat_[hohi_[ic]]=(doSaturation_)?mgain:99999.;
         
           gain*=tpgfactor;
           peds_[hohi_[ic]]=ped;
           gains_[hohi_[ic]]=gain;
         }
     
       initializedHO_=true;
     }
 
     if(!initializedHF_ && det_==6){
       nhfcells_ = hfhi_.size();
       //HF
       for(unsigned ic=0;ic<nhfcells_;++ic)
         {
           float ped=theDbService->getPedestal(theDetIds_[hfhi_[ic]])->getValue(0);
           float gain=theDbService->getGain(theDetIds_[hfhi_[ic]])->getValue(0);
           float mgain=0.;
           for(unsigned ig=0;ig<4;++ig)
             mgain+=theDbService->getGain(theDetIds_[hfhi_[ic]])->getValue(ig);
           // additional 1/2 factor for the HF (Salavat)
           if(noiseFromDb_[0])
             {
               // computation from when the noise was taken 
               noisesigma_[hfhi_[ic]]=noiseInfCfromDB(theDbService,theDetIds_[hfhi_[ic]])*mgain*0.25;
             }
           //      std::cout << " NOISEHF " << theDetIds_[hfhi_[ic]].ieta() << " " << noisesigma_[hfhi_[ic]] << "  "<< std::endl;
           
           mgain*=2500./0.36;
           sat_[hfhi_[ic]]=(doSaturation_)?mgain:99999.;
           int ieta=theDetIds_[hfhi_[ic]].ieta();
           gain*=TPGFactor_[(ieta>0)?ieta+45:-ieta+2];
           peds_[hfhi_[ic]]=ped;
           gains_[hfhi_[ic]]=gain;
         }
 
       initializedHF_=true;
     }
   
   // clear the vector we don't need. It is a bit stupid 
   hcalRecHits_.resize(maxIndex_+1,0.);
 
 }

void HcalRecHitsMaker::loadHcalRecHits	(	edm::Event &	iEvent,
		const HcalTopology &	topo,
		HBHERecHitCollection &	hbheHits,
		HBHEDigiCollection &	hbheDigis,
		RandomEngineAndDistribution const *	random
	)

Definition at line 397 of file HcalRecHitsMaker.cc.

References cond::rpcobgas::detid, doDigis_, relval_parameters_module::energy, fCtoAdc(), firedCells_, gains_, HcalRespCorr::getValue(), HcalCondObjectContainer< Item >::getValues(), HcalBarrel, hcalRecHits_, loadPCaloHits(), myRespCorr, noisify(), peds_, edm::SortedCollection< T, SORT >::push_back(), edm::SortedCollection< T, SORT >::reserve(), sat_, HBHEDataFrame::setSample(), HBHEDataFrame::setSize(), HcalDetId::subdet(), theDetIds_, and threshold_.

Referenced by CaloRecHitsProducer::produce().

 {
   loadPCaloHits(iEvent,topo, random);
   noisify(random);
   hbheHits.reserve(firedCells_.size());
   if(doDigis_)
     {
       hbheDigis.reserve(firedCells_.size());
     }
   static HcalQIESample zeroSample(0,0,0,0);
   unsigned nhits=firedCells_.size();
   // HB and HE
 
   for(unsigned ihit=0;ihit<nhits;++ihit)
     {
       unsigned cellhashedindex=firedCells_[ihit];
       const HcalDetId& detid  = theDetIds_[cellhashedindex];
       unsigned subdet=(detid.subdet()==HcalBarrel) ? 0: 1;  
 
       float energy=hcalRecHits_[cellhashedindex];
       // Check if it is above the threshold
       if(energy<threshold_[subdet]) continue; 
       // apply RespCorr only to the RecHit
       energy *= myRespCorr->getValues(theDetIds_[cellhashedindex])->getValue();
       // poor man saturation
       if(energy>sat_[cellhashedindex]) 
     {
       //      std::cout << " Saturation " << energy << " " << sat_[cellhashedindex] << " HB " << std::endl;
       energy=sat_[cellhashedindex];
     }
       
 
       hbheHits.push_back(HBHERecHit(detid,energy,0.));      
       if(doDigis_)
     {
       HBHEDataFrame myDataFrame(detid);
       myDataFrame.setSize(2);
       double nfc=hcalRecHits_[cellhashedindex]/gains_[cellhashedindex]+peds_[cellhashedindex];
       int nadc=fCtoAdc(nfc);
       HcalQIESample qie(nadc, 0, 0, 0) ;
       myDataFrame.setSample(0,qie);
       myDataFrame.setSample(1,zeroSample);
       hbheDigis.push_back(myDataFrame);
     }
     }
 }

void HcalRecHitsMaker::loadHcalRecHits	(	edm::Event &	iEvent,
		const HcalTopology &	topo,
		HORecHitCollection &	ho,
		HODigiCollection &	hoDigis,
		RandomEngineAndDistribution const *	random
	)

Definition at line 447 of file HcalRecHitsMaker.cc.

References cond::rpcobgas::detid, doDigis_, relval_parameters_module::energy, firedCells_, HcalRespCorr::getValue(), HcalCondObjectContainer< Item >::getValues(), hcalRecHits_, HcalDetId::ieta(), loadPCaloHits(), myRespCorr, noisify(), edm::SortedCollection< T, SORT >::push_back(), edm::SortedCollection< T, SORT >::reserve(), sat_, theDetIds_, and threshold_.

 {
   loadPCaloHits(iEvent,topo, random);
   noisify(random);
   hoHits.reserve(firedCells_.size());
   if(doDigis_)
     {
       hoDigis.reserve(firedCells_.size());
     }
   static HcalQIESample zeroSample(0,0,0,0);
 
   // HO
   unsigned nhits=firedCells_.size();
   for(unsigned ihit=0;ihit<nhits;++ihit)
     {
 
       unsigned cellhashedindex=firedCells_[ihit];
       // Check if it is above the threshold
       if(hcalRecHits_[cellhashedindex]<threshold_[0]) continue; 
 
       const HcalDetId&  detid=theDetIds_[cellhashedindex];
       int ieta = detid.ieta();
       
       // Force  only Ring#0 to remain
       if(ieta > 4 || ieta < -4 ) continue;
 
       float energy=hcalRecHits_[cellhashedindex];
       // apply RespCorr
       energy *= myRespCorr->getValues(theDetIds_[cellhashedindex])->getValue();
 
       // poor man saturation
       if(energy>sat_[cellhashedindex]) energy=sat_[cellhashedindex];
 
       hoHits.push_back(HORecHit(detid,energy,0));
     }
 }

void HcalRecHitsMaker::loadHcalRecHits	(	edm::Event &	iEvent,
		const HcalTopology &	topo,
		HFRecHitCollection &	hfHits,
		HFDigiCollection &	hfDigis,
		RandomEngineAndDistribution const *	random
	)

Definition at line 486 of file HcalRecHitsMaker.cc.

References cond::rpcobgas::detid, doDigis_, relval_parameters_module::energy, fCtoAdc(), firedCells_, gains_, HcalRespCorr::getValue(), HcalCondObjectContainer< Item >::getValues(), hcalRecHits_, loadPCaloHits(), myRespCorr, noisify(), peds_, edm::SortedCollection< T, SORT >::push_back(), edm::SortedCollection< T, SORT >::reserve(), sat_, HFDataFrame::setSample(), HFDataFrame::setSize(), theDetIds_, and threshold_.

 {
   loadPCaloHits(iEvent,topo, random);
   noisify(random);
   hfHits.reserve(firedCells_.size());
   if(doDigis_)
     {
       hfDigis.reserve(firedCells_.size());
     }
   static HcalQIESample zeroSample(0,0,0,0);
 
   unsigned nhits=firedCells_.size();
   for(unsigned ihit=0;ihit<nhits;++ihit)
     {
       unsigned cellhashedindex=firedCells_[ihit];
       // Check if it is above the threshold
       if(hcalRecHits_[cellhashedindex]<threshold_[0]) continue; 
 
       float energy=hcalRecHits_[cellhashedindex];
 
       // apply RespCorr
       energy *= myRespCorr->getValues(theDetIds_[cellhashedindex])->getValue();
 
       // poor man saturation
       if(energy>sat_[cellhashedindex]) energy=sat_[cellhashedindex];
 
       const HcalDetId & detid=theDetIds_[cellhashedindex];
       hfHits.push_back(HFRecHit(detid,energy,0.));      
       if(doDigis_)
     {
       HFDataFrame myDataFrame(detid);
       myDataFrame.setSize(1);
 
       double nfc= hcalRecHits_[cellhashedindex]/gains_[cellhashedindex]+peds_[cellhashedindex];
       int nadc=fCtoAdc(nfc/2.6);
       HcalQIESample qie(nadc, 0, 0, 0) ;
       myDataFrame.setSample(0,qie);
       hfDigis.push_back(myDataFrame);
     }
     }
 }

void HcalRecHitsMaker::loadPCaloHits	(	const edm::Event &	iEvent,
		const HcalTopology &	topo,
		RandomEngineAndDistribution const *	random
	)

private

Definition at line 339 of file HcalRecHitsMaker.cc.

References clean(), corrfac_, counter, det_, cond::rpcobgas::detid, HcalTopology::detId2denseId(), Fill(), firedCells_, fractionOOT(), edm::Event::getByLabel(), HcalBarrel, HcalEndcap, HcalForward, HcalOuter, inputCol_, noise_, edm::Handle< T >::product(), and random.

Referenced by loadHcalRecHits().

 {
   clean();
 
   edm::Handle<CrossingFrame<PCaloHit> > cf;
   iEvent.getByLabel(inputCol_,cf);
   std::auto_ptr<MixCollection<PCaloHit> > colcalo(new MixCollection<PCaloHit>(cf.product(),std::pair<int,int>(0,0) ));
 
   MixCollection<PCaloHit>::iterator it=colcalo->begin();;
   MixCollection<PCaloHit>::iterator itend=colcalo->end();
   unsigned counter=0;
   for(;it!=itend;++it)
     {
       HcalDetId detid(it->id());
       int hashedindex=topo.detId2denseId(detid);
 
       // apply ToF correction
       int time_slice=0; // temporary
       double fTOF=1.;
       if (detid.subdet()==HcalForward) fTOF = (time_slice==0) ? 1. : 0.;    
       else fTOF = fractionOOT(time_slice);
 
       switch(detid.subdet())
     {
     case HcalBarrel: 
       {
         if(det_==4)
           Fill(hashedindex,fTOF*(it->energy()),firedCells_,noise_[0],corrfac_[0], random);
       }
       break;
     case HcalEndcap: 
       {   
         if(det_==4)
           Fill(hashedindex,fTOF*(it->energy()),firedCells_,noise_[1],corrfac_[1], random);
       }
       break;
     case HcalOuter: 
       {
         if(det_==5)
           Fill(hashedindex,fTOF*(it->energy()),firedCells_,noise_[0],corrfac_[0], random);
       }
       break;             
     case HcalForward: 
       {
         if(det_==6 && time_slice==0) // skip the HF hit if out-of-time
           Fill(hashedindex,it->energy(),firedCells_,noise_[0],corrfac_[0], random);
       }
       break;
     default:
       edm::LogWarning("CaloRecHitsProducer") << "RecHit not registered\n";
       ;
     }
       ++counter;
     }
 }

double HcalRecHitsMaker::noiseInfCfromDB	(	const HcalDbService *	conditions,
		const HcalDetId &	detId
	)

private

Definition at line 710 of file HcalRecHitsMaker.cc.

References funct::abs(), alpha, beta, f, HcalPedestalWidth::getSigma(), HcalDetId::ieta(), mathSSE::sqrt(), and HcalDetId::subdet().

Referenced by init().

 {
   // method from Salavat
   // fetch pedestal widths (noise sigmas for all 4 CapID)
   const HcalPedestalWidth* pedWidth =
     conditions-> getPedestalWidth(detId); 
   double ssqq_1 = pedWidth->getSigma(0,0);
   double ssqq_2 = pedWidth->getSigma(1,1);
   double ssqq_3 = pedWidth->getSigma(2,2);
   double ssqq_4 = pedWidth->getSigma(3,3);
 
   // correction factors (hb,he,ho,hf)
   //  static float corrfac[4]={1.39,1.32,1.17,3.76};
   //  static float corrfac[4]={0.93,0.88,1.17,2.51}; // divided HB, HE, HF by 1.5 (to take into account the halving of the number of time slices), HO did not change
 
   int sub   = detId.subdet();
 
   // HO: only Ring#0 matters 
   int ieta  = detId.ieta();
   if (sub == 3 && abs (ieta) > 4) return 0.;   
 
   // effective RecHits (for this particular detId) noise calculation :
   
   double sig_sq_mean =  0.25 * ( ssqq_1 + ssqq_2 + ssqq_3 + ssqq_4);
   
   // f - external parameter (currently set to 0.5 in the FullSim) !!!
   double f=0.5;
 
   double term  = sqrt (1. + sqrt(1. - 2.*f*f));
   double alpha = sqrt (0.5) * term;
   double beta  = sqrt (0.5) * f / term;
 
   //  double RMS1   = sqrt(sig_sq_mean) * sqrt (2.*beta*beta + alpha*alpha) ;
   double RMS4   = sqrt(sig_sq_mean) * sqrt (2.*beta*beta + 2.*(alpha-beta)*(alpha-beta) + 2.*(alpha-2.*beta)*(alpha-2.*beta)) ;
 
   double noise_rms_fC;
 
   //  if(sub == 4)  noise_rms_fC = RMS1;
   //  else          noise_rms_fC = RMS4;
   noise_rms_fC = RMS4;
 
 
   // to convert from above fC to GeV - multiply by gain (GeV/fC)        
   //  const HcalGain*  gain = conditions->getGain(detId); 
   //  double noise_rms_GeV = noise_rms_fC * gain->getValue(0); // Noise RMS (GeV) for detId channel
   return noise_rms_fC;
 }

void HcalRecHitsMaker::noisify ( RandomEngineAndDistribution const * random )

private

Definition at line 590 of file HcalRecHitsMaker.cc.

References corrfac_, det_, firedCells_, hbhi_, hcalHotFraction_, hcalRecHits_, hehi_, hfhi_, hohi_, myGaussianTailGenerators_, nhbcells_, nhecells_, nhfcells_, nhocells_, noise_, noisifySubdet(), threshold_, and pileupDistInMC::total.

Referenced by loadHcalRecHits().

 {
   unsigned total=0;
   switch(det_)
     {
     case 4:
       {
     // do the HB
     if(noise_[0] != 0.) {
       total+=noisifySubdet(hcalRecHits_,firedCells_,hbhi_,nhbcells_,hcalHotFraction_[0],myGaussianTailGenerators_[0],noise_[0],threshold_[0],corrfac_[0], random);
     }
     // do the HE
     if(noise_[1] != 0.) {
       total+=noisifySubdet(hcalRecHits_,firedCells_,hehi_,nhecells_,hcalHotFraction_[1],myGaussianTailGenerators_[1],noise_[1],threshold_[1],corrfac_[1], random);
     }
       }
       break;
     case 5:
       {
     // do the HO
     if(noise_[0] != 0.) {
       total+=noisifySubdet(hcalRecHits_,firedCells_,hohi_,nhocells_,hcalHotFraction_[0],myGaussianTailGenerators_[0],noise_[0],threshold_[0],corrfac_[0], random);
     }
       }
       break;
     case 6:
       {
     // do the HF
     if(noise_[0] != 0.) {
       total+=noisifySubdet(hcalRecHits_,firedCells_,hfhi_,nhfcells_,hcalHotFraction_[0],myGaussianTailGenerators_[0],noise_[0],threshold_[0],corrfac_[0], random);
     }
       }
       break;
     default:
       break;
     }
   edm::LogInfo("CaloRecHitsProducer") << "CaloRecHitsProducer : added noise in "<<  total << " HCAL cells "  << std::endl;
 }

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,
		double	correctionfactor,
		RandomEngineAndDistribution const *	random
	)

private

Definition at line 629 of file HcalRecHitsMaker.cc.

References RandomEngineAndDistribution::flatShoot(), RandomEngineAndDistribution::gaussShoot(), hcalRecHits_, timingPdfMaker::mean, noise_, noisesigma_, RandomEngineAndDistribution::poissonShoot(), and GaussianTail::shoot().

Referenced by noisify().

 {
  // If the fraction of "hot " is small, use an optimized method to inject noise only in noisy cells. The 30% has not been tuned
   if(noise_!=-1 && hcalHotFraction==0.) return 0;
   if(hcalHotFraction<0.3 && noise_!=-1)
     {
       double mean = (double)(ncells-theHits.size())*hcalHotFraction;
       unsigned nhcal = random->poissonShoot(mean);
   
       unsigned ncell=0;
       unsigned cellindex=0;
       uint32_t cellhashedindex=0;
       
       while(ncell < nhcal)
     {
       cellindex = (unsigned)(random->flatShoot()*ncells);
       cellhashedindex = thecells[cellindex];
 
       if(hcalRecHits_[cellhashedindex]==0.) // new cell
         {
           hcalRecHits_[cellhashedindex]=myGT->shoot(random);
           theHits.push_back(cellhashedindex);
           ++ncell;
         }
     }
       return ncell;
     }
   else // otherwise, inject noise everywhere
     {
       uint32_t cellhashedindex=0;
       unsigned nhcal=thecells.size();
 
       double sigma = noise_;
       for(unsigned ncell=0;ncell<nhcal;++ncell)
     {
       cellhashedindex = thecells[ncell];
       if(hcalRecHits_[cellhashedindex]==0.) // new cell
         { 
           //only use noisesigma if it has been filled from the DB
           //otherwise, leave sigma alone (default is user input value)
           if(noise_==-1) sigma=noisesigma_[cellhashedindex]*correctionfactor;
 
           double noise =random->gaussShoot(0.,sigma);
           if(noise>threshold)
         {
           hcalRecHits_[cellhashedindex]=noise;          
           theHits.push_back(cellhashedindex);
         }
         }
     }
       return nhcal;
     }
   return 0;
 }