Skip to content

Commit

Permalink
Merge pull request #23894 from hvanhaev/from-CMSSW_8_0_31
Browse files Browse the repository at this point in the history 
backport of CASTOR updates on rechitcorrector and noise simulation
  • Loading branch information
@cmsbuild
cmsbuild committed Sep 11, 2018
2 parents 36b143d + 33c5a78 commit 2c691cc
Show file tree
Hide file tree
Showing 5 changed files with 74 additions and 109 deletions.
73 changes: 24 additions & 49 deletions RecoLocalCalo/Castor/src/RecHitCorrector.cc
View file
Expand Up @@ -4,9 +4,9 @@
// Class: RecHitCorrector
//
/**\class RecHitCorrector RecHitCorrector.cc RecoLocalCalo/Castor/src/RecHitCorrector.cc
Description: [one line class summary]
Implementation:
[Notes on implementation]
*/
Expand All @@ -22,7 +22,7 @@

// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EDProducer.h"
#include "FWCore/Framework/interface/stream/EDProducer.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
Expand All @@ -44,15 +44,13 @@
// class declaration
//

class RecHitCorrector : public edm::EDProducer {
class RecHitCorrector : public edm::stream::EDProducer<> {
public:
explicit RecHitCorrector(const edm::ParameterSet&);
~RecHitCorrector();
~RecHitCorrector() override;

private:
virtual void beginJob() override ;
virtual void produce(edm::Event&, const edm::EventSetup&) override;
virtual void endJob() override ;
void produce(edm::Event&, const edm::EventSetup&) override;

// ----------member data ---------------------------
edm::EDGetTokenT<CastorRecHitCollection> tok_input_;
Expand Down Expand Up @@ -114,72 +112,49 @@ RecHitCorrector::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
iSetup.get<CastorChannelQualityRcd>().get(p);
CastorChannelQuality* myqual = new CastorChannelQuality(*p.product());

if (!rechits.isValid()) std::cout << "No valid CastorRecHitCollection found, please check the InputLabel..." << std::endl;
if (!rechits.isValid()) edm::LogWarning("CastorRecHitCorrector") << "No valid CastorRecHitCollection found, please check the InputLabel...";

CastorCalibrations calibrations;

std::auto_ptr<CastorRecHitCollection> rec(new CastorRecHitCollection);
auto rec = std::make_unique<CastorRecHitCollection>();

for (unsigned int i=0;i<rechits->size();i++) {
CastorRecHit rechit = (*rechits)[i];
//std::cout << "rechit energy = " << rechit.energy() << std::endl;
double fC = factor_*rechit.energy();
double time = rechit.time();
//std::cout << "rechit energy(fC) = " << fC << " time = " << time << std::endl;

// do proper gain calibration reading the latest entries in the condDB
const CastorCalibrations& calibrations=conditions->getCastorCalibrations(rechit.id());
int capid = 0; // take some capid, gains are the same for all capid's
double correctedenergy = factor_*rechit.energy();

double correctedenergy = 0;
if (doInterCalib_) {
if (rechit.id().module() <= 2) {
correctedenergy = 0.5*fC*calibrations.gain(capid);
//std::cout << " correctedenergy = " << correctedenergy << " gain = " << calibrations.gain(capid) << std::endl;
} else {
correctedenergy = fC*calibrations.gain(capid);
}
} else {
if (rechit.id().module() <= 2) {
correctedenergy = 0.5*fC;
} else {
correctedenergy = fC;
}
// do proper gain calibration reading the latest entries in the condDB
const CastorCalibrations& calibrations=conditions->getCastorCalibrations(rechit.id());
int capid = 0; // take some capid, gains are the same for all capid's
correctedenergy *= calibrations.gain(capid);
}

// now check the channelquality of this rechit
bool ok = true;
DetId detcell=(DetId)rechit.id();
std::vector<DetId> channels = myqual->getAllChannels();
//std::cout << "number of specified quality flags = " << channels.size() << std::endl;
for (std::vector<DetId>::iterator channel = channels.begin();channel != channels.end();channel++) {
if (channel->rawId() == detcell.rawId()) {
const CastorChannelStatus* mydigistatus=myqual->getValues(*channel);
//std::cout << "CastorChannelStatus = " << mydigistatus->getValue() << std::endl;
if (mydigistatus->getValue() == 2989) ok = false; // 2989 = BAD
for (auto channel : channels) {
if (channel.rawId() == detcell.rawId()) {
const CastorChannelStatus* mydigistatus=myqual->getValues(channel);
if (mydigistatus->getValue() == 2989) {
ok = false; // 2989 = BAD
break;
}
}
}

if (ok) {
CastorRecHit *correctedhit = new CastorRecHit(rechit.id(),correctedenergy,time);
rec->push_back(*correctedhit);
rec->emplace_back(rechit.id(),correctedenergy,time);
}
}

iEvent.put(rec);
iEvent.put(std::move(rec));

delete myqual;

}

// ------------ method called once each job just before starting event loop ------------
void
RecHitCorrector::beginJob()
{
}

// ------------ method called once each job just after ending the event loop ------------
void
RecHitCorrector::endJob() {
}

//define this as a plug-in
DEFINE_FWK_MODULE(RecHitCorrector);
37 changes: 22 additions & 15 deletions SimCalorimetry/CastorSim/src/CastorAmplifier.cc
View file
Expand Up @@ -14,7 +14,7 @@
#include<iostream>

CastorAmplifier::CastorAmplifier(const CastorSimParameterMap * parameters, bool addNoise) :
theDbService(0),
theDbService(nullptr),
theParameterMap(parameters),
theStartingCapId(0),
addNoise_(addNoise)
Expand All @@ -23,30 +23,37 @@ CastorAmplifier::CastorAmplifier(const CastorSimParameterMap * parameters, bool

void CastorAmplifier::amplify(CaloSamples & frame, CLHEP::HepRandomEngine* engine) const {
const CastorSimParameters & parameters = theParameterMap->castorParameters();
assert(theDbService != 0);
assert(theDbService != nullptr);
HcalGenericDetId hcalGenDetId(frame.id());
const CastorPedestal* peds = theDbService->getPedestal(hcalGenDetId);
const CastorPedestalWidth* pwidths = theDbService->getPedestalWidth(hcalGenDetId);
if (!peds || !pwidths )
{
edm::LogError("CastorAmplifier") << "Could not fetch HCAL/CASTOR conditions for channel " << hcalGenDetId;
}
else
{
double gauss [32]; //big enough
double noise [32]; //big enough
double fCperPE = parameters.photoelectronsToAnalog(frame.id());

double gauss [32]; //big enough
double noise [32]; //big enough
double fCperPE = parameters.photoelectronsToAnalog(frame.id());

for (int i = 0; i < frame.size(); i++) gauss[i] = CLHEP::RandGaussQ::shoot(engine, 0., 1.);
pwidths->makeNoise (frame.size(), gauss, noise);
for(int tbin = 0; tbin < frame.size(); ++tbin) {
int capId = (theStartingCapId + tbin)%4;
double pedestal = peds->getValue (capId);
for (int i = 0; i < frame.size(); i++) { gauss[i] = CLHEP::RandGaussQ::shoot(engine, 0., 1.); }
if(addNoise_) {
pedestal += noise [tbin];
pwidths->makeNoise (frame.size(), gauss, noise);
}
for(int tbin = 0; tbin < frame.size(); ++tbin) {
int capId = (theStartingCapId + tbin)%4;
double pedestal = peds->getValue(capId);
if(addNoise_) {
if (parameters.doDynamicNoise()) {
pedestal += noise [tbin]*(fCperPE/parameters.getNominalfCperPE());
} else {
pedestal += noise [tbin];
}
}
frame[tbin] *= fCperPE;
frame[tbin] += pedestal;
}
frame[tbin] *= fCperPE;
frame[tbin] += pedestal;
}
LogDebug("CastorAmplifier") << frame;
}

51 changes: 22 additions & 29 deletions SimCalorimetry/CastorSim/src/CastorSimParameters.cc
View file
Expand Up @@ -11,64 +11,57 @@
CastorSimParameters::CastorSimParameters(double simHitToPhotoelectrons, double photoelectronsToAnalog,double samplingFactor, double timePhase, bool syncPhase)
: CaloSimParameters(simHitToPhotoelectrons, photoelectronsToAnalog, samplingFactor, timePhase, 6, 4, false, syncPhase),
theDbService(0),
theSamplingFactor( samplingFactor )
theSamplingFactor( samplingFactor ),
nominalfCperPE( 1),
dynamicNoise(false)
{
}


CastorSimParameters::CastorSimParameters(const edm::ParameterSet & p)
: CaloSimParameters(p),
theDbService(0),
theSamplingFactor( p.getParameter<double>("samplingFactor") )
theSamplingFactor( p.getParameter<double>("samplingFactor") ),
nominalfCperPE( p.getParameter<double>("photoelectronsToAnalog") ),
dynamicNoise(p.getParameter<bool>("doDynamicNoise") )
{
}

/*
double CastorSimParameters::simHitToPhotoelectrons(const DetId & detId) const
double CastorSimParameters::getNominalfCperPE() const
{
// the gain is in units of GeV/fC. We want a constant with pe/dGeV
// pe/dGeV = (GeV/dGeV) / (GeV/fC) / (fC/pe)
double result = CaloSimParameters::simHitToPhotoelectrons(detId);
if(HcalGenericDetId(detId).genericSubdet() != HcalGenericDetId::HcalGenForward
|| HcalGenericDetId(detId).genericSubdet() != HcalGenericDetId::HcalGenCastor)
{
result = samplingFactor(detId) / fCtoGeV(detId) / photoelectronsToAnalog();
}
return result;
// return the nominal PMT gain value of CASTOR from the config file.
return nominalfCperPE;
}

bool CastorSimParameters::doDynamicNoise() const
{
// activate proper noise treatment depending on used gain values.
return dynamicNoise;
}
*/

double CastorSimParameters::photoelectronsToAnalog(const DetId & detId) const
{
// calculate factor (PMT gain) using sampling factor value & available electron gain
//std::cout << " sampling factor = " << theSamplingFactor << " fCtoGeV = " << fCtoGeV(detId) << " and photoelectronsToAnalog = " << theSamplingFactor/fCtoGeV(detId) << std::endl;
return theSamplingFactor/fCtoGeV(detId);
}



double CastorSimParameters::fCtoGeV(const DetId & detId) const
{
assert(theDbService != 0);
assert(theDbService != nullptr);
HcalGenericDetId hcalGenDetId(detId);
const CastorGain* gains = theDbService->getGain(hcalGenDetId);
const CastorGainWidth* gwidths = theDbService->getGainWidth(hcalGenDetId);
double result = 0.0;
if (!gains || !gwidths )
{
edm::LogError("CastorAmplifier") << "Could not fetch HCAL conditions for channel " << hcalGenDetId;
}
// only one gain will be recorded per channel, so just use capID 0 for now

double result = gains->getValue(0);
// if(doNoise_)
/// {
// result += CLHEP::RandGaussQ::shoot(0., gwidths->getValue(0));
// }
else
{
// only one gain will be recorded per channel, so just use capID 0 for now
result = gains->getValue(0);
}
return result;
}
/*
double CastorSimParameters::samplingFactor(const DetId & detId) const {
HcalDetId hcalDetId(detId);
return theSamplingFactors[hcalDetId.ietaAbs()-theFirstRing];
}
*/
19 changes: 4 additions & 15 deletions SimCalorimetry/CastorSim/src/CastorSimParameters.h
View file
Expand Up @@ -12,34 +12,23 @@ class CastorSimParameters : public CaloSimParameters
CastorSimParameters(double simHitToPhotoelectrons, double photoelectronsToAnalog, double samplingFactor, double timePhase, bool syncPhase);
CastorSimParameters(const edm::ParameterSet & p);

/*
CastorSimParameters(double simHitToPhotoelectrons, double photoelectronsToAnalog,
double samplingFactor, double timePhase,
int readoutFrameSize, int binOfMaximum,
bool doPhotostatistics, bool syncPhase,
int firstRing, const std::vector<double> & samplingFactors);
CastorSimParameters(const edm::ParameterSet & p);
*/

virtual ~CastorSimParameters() {}

void setDbService(const CastorDbService * service) {theDbService = service;}

//virtual double simHitToPhotoelectrons(const DetId & detId) const;
double getNominalfCperPE() const;
bool doDynamicNoise() const;

virtual double photoelectronsToAnalog(const DetId & detId) const;

double fCtoGeV(const DetId & detId) const;

/// the ratio of actual incident energy to deposited energy
/// in the SimHit
// virtual double samplingFactor(const DetId & detId) const;


private:
const CastorDbService * theDbService;
double theSamplingFactor;
//std::vector<double> theSamplingFactors;
double nominalfCperPE;
bool dynamicNoise;
};

#endif
Expand Down
3 changes: 2 additions & 1 deletion SimGeneral/MixingModule/python/castorDigitizer_cfi.py
View file
Expand Up @@ -15,6 +15,7 @@
photoelectronsToAnalog = cms.double(4.009),
simHitToPhotoelectrons = cms.double(1000.0),
syncPhase = cms.bool(True),
timePhase = cms.double(-4.0)
timePhase = cms.double(-4.0),
doDynamicNoise = cms.bool(False)
)
)

0 comments on commit 2c691cc

Please sign in to comment.