Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

backport of CASTOR updates on rechitcorrector and noise simulation #23894

Merged
merged 4 commits into from Sep 11, 2018
Merged

backport of CASTOR updates on rechitcorrector and noise simulation #23894

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
0a1daca
backport of CASTOR updates on rechitcorrector and noise simulation
hvanhaev Jul 9, 2018
96cc6ae
make the noise bug fix available as parameter to restore default beha…
hvanhaev Aug 27, 2018
49ca4f1
forgot to set parameter to False as default
hvanhaev Aug 27, 2018
33c5a78
revert photoelectronsToAnalog parameter to old version for safe backp…
hvanhaev Sep 10, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
73 changes: 24 additions & 49 deletions RecoLocalCalo/Castor/src/RecHitCorrector.cc
View file
Open in desktop
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);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@hvanhaev
please remind me the reason for hits with the id().module() <= 2 to be not handled in a special way anymore.
Are hits with these module values present only for some period of data taking. If so, when?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@hvanhaev
In case you missed the question, please clarify on the hits with id().module() <= 2

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@hvanhaev: (maybe I missed some part of the discussion about it, sorry for that) shouldn't this 0.5 factor for id().module() <= 2 be restored to fully reproduce the previous 80X behaviour?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

#23894 (comment)
this code is not used in the standard reco workflows.
So, there is no requirement for unchanged default behavior.

//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);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@hvanhaev , is it correct? I would think about "rec->emplace_back(new CastorRecHit(....))"

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There is no need for new. One could think of explicitly calling a constructor emplace_back(CastorRecHit(..., but that's the same as what it is now.

the old implementation collected copies of the hits (type CastorRecHit). The updated version does exactly the same. and also avoids a memory leak that was there before.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No, using new here is totally unnecessary and in fact it's a memory leak. rec is a CastorRecHitCollection and it can construct the CastorRecHit objects from arguments using emplace_back().

}
}

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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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)
)
)