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.

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

Update Phase 1 pixel gains in digitizer/clusterizer #19181

Merged
merged 8 commits into from Jun 15, 2017
21 changes: 21 additions & 0 deletions HLTrigger/Configuration/python/customizeHLTforCMSSW.py
Expand Up @@ -173,6 +173,22 @@ def customiseFor19098(process):
if hasattr(producer, "infinitesimalPt"): del producer.infinitesimalPt
return process

def customiseFor19181_pixel_phase0(process):
for producer in producers_by_type(process, "SiPixelClusterProducer"):
producer.VCaltoElectronGain_L1 = cms.int32(65)
producer.VCaltoElectronOffset_L1 = cms.int32(-414)
producer.ClusterThreshold = cms.int32(4000)
producer.ClusterThreshold_L1 = cms.int32(4000)

Copy link
Contributor

Choose a reason for hiding this comment

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

forgot return process

def customiseFor19181_pixel_phase1(process):
for producer in producers_by_type(process, "SiPixelClusterProducer"):
producer.VCaltoElectronGain = cms.int32(47)
producer.VCaltoElectronGain_L1 = cms.int32(50)
producer.VCaltoElectronOffset = cms.int32(-60)
producer.VCaltoElectronOffset_L1 = cms.int32(-670)
producer.ClusterThreshold = cms.int32(4000)
producer.ClusterThreshold_L1 = cms.int32(2000)

Copy link
Contributor

Choose a reason for hiding this comment

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

forgot return process

# CMSSW version specific customizations
def customizeHLTforCMSSW(process, menuType="GRun"):

Expand All @@ -189,4 +205,9 @@ def customizeHLTforCMSSW(process, menuType="GRun"):
process = customiseFor18559(process)
process = customiseFor19098(process)

if (menuType == "GRun2016"):
process = customiseFor19181_pixel_phase0(process)
else:
process = customiseFor19181_pixel_phase1(process)
Copy link
Contributor

Choose a reason for hiding this comment

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

what will happen to run1 workflows? Do they also have GRun2016?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I only found phase0 configurations in GRun2016 configs. An HLT expert should comment

Copy link
Contributor

Choose a reason for hiding this comment

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

Run1 workflows use a stripped HLT menu with nearly no physics.
From the HLT technical point this looks ok, but I can't judge the physics/detector...


return process
Expand Up @@ -47,9 +47,12 @@ PixelThresholdClusterizer::PixelThresholdClusterizer
// Get thresholds in electrons
thePixelThreshold( conf.getParameter<int>("ChannelThreshold") ),
theSeedThreshold( conf.getParameter<int>("SeedThreshold") ),
theClusterThreshold( conf.getParameter<double>("ClusterThreshold") ),
theClusterThreshold( conf.getParameter<int>("ClusterThreshold") ),
theClusterThreshold_L1( conf.getParameter<int>("ClusterThreshold_L1") ),
theConversionFactor( conf.getParameter<int>("VCaltoElectronGain") ),
theConversionFactor_L1( conf.getParameter<int>("VCaltoElectronGain_L1") ),
theOffset( conf.getParameter<int>("VCaltoElectronOffset") ),
theOffset_L1( conf.getParameter<int>("VCaltoElectronOffset_L1") ),
theStackADC_( conf.exists("AdcFullScaleStack") ? conf.getParameter<int>("AdcFullScaleStack") : 255 ),
theFirstStack_( conf.exists("FirstStackLayer") ? conf.getParameter<int>("FirstStackLayer") : 5 ),
theElectronPerADCGain_( conf.exists("ElectronPerADCGain") ? conf.getParameter<double>("ElectronPerADCGain") : 135. ),
Expand All @@ -63,6 +66,28 @@ PixelThresholdClusterizer::PixelThresholdClusterizer
/////////////////////////////////////////////////////////////////////////////
PixelThresholdClusterizer::~PixelThresholdClusterizer() {}


// Configuration descriptions
void
PixelThresholdClusterizer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
// siPixelClusters
edm::ParameterSetDescription desc;
desc.add<edm::InputTag>("src", edm::InputTag("siPixelDigis"));
desc.add<int>("ChannelThreshold", 1000);
desc.addUntracked<bool>("MissCalibrate", true);
desc.add<bool>("SplitClusters", false);
desc.add<int>("VCaltoElectronGain", 65);
desc.add<int>("VCaltoElectronGain_L1", 65);
desc.add<int>("VCaltoElectronOffset", -414);
desc.add<int>("VCaltoElectronOffset_L1", -414);
desc.add<std::string>("payloadType", "Offline");
desc.add<int>("SeedThreshold", 1000);
desc.add<int>("ClusterThreshold_L1", 4000);
desc.add<int>("ClusterThreshold", 4000);
descriptions.add("siPixelClusters", desc);
desc.add<int>("maxNumberOfClusters", -1);
}
Copy link
Contributor

Choose a reason for hiding this comment

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

descriptions.add must be last...

Copy link
Contributor

Choose a reason for hiding this comment

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

If I'm not mistaken, currently this method is a dead code.
Perhaps in the future refactoring this can become useful.


//----------------------------------------------------------------------------
//! Prepare the Clusterizer to work on a particular DetUnit. Re-init the
//! size of the panel/plaquette (so update nrows and ncols),
Expand Down Expand Up @@ -119,6 +144,11 @@ void PixelThresholdClusterizer::clusterizeDetUnitT( const T & input,

detid_ = input.detId();

// Set separate cluster threshold for L1 (needed for phase1)
auto clusterThreshold = theClusterThreshold;
int layer = (DetId(detid_).subdetId()==1) ? PXBDetId(detid_).layer() : 0;
Copy link
Contributor

Choose a reason for hiding this comment

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

@jkarancs are you sure you can use the PXBDetId class for the phase 1 detector? Between phase 0 and phase 1 the detid definition has changed and you should use TrackerTopology.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

@venturia
This class was already using PXBDetId to get the layer number elsewhere. TTopo gives the same value for layer. It indeed could be changed to use TrackerTopology, but EventSetup is not easily accessible in this class. The class which inherits from this could pass it I think.

Copy link
Contributor

Choose a reason for hiding this comment

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

This is strange. According to this page: https://github.com/cms-sw/cmssw/tree/master/Geometry/TrackerNumberingBuilder
in the old detector the layer number is encoded in the detid starting at bit 16 while in the phase 1 detector the layer number starts at bit 20. According to that it is not possible that you get the right layer number if you decode phase 1 detids with the PXBDetId class.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Uh, you are absolutely right, the bit coding is only same for phase 0. We must use TrackerTopology indeed.

Copy link
Contributor

Choose a reason for hiding this comment

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

Do you see any difference in the results after the fix of the DetId decoding? Is the problem also in the code used for the present data reconstruction? Do you know or can you check if these outdated class is used in another parts of the pixel reco code?

if (layer==1) clusterThreshold = theClusterThreshold_L1;

// Copy PixelDigis to the buffer array; select the seed pixels
// on the way, and store them in theSeeds.
copy_to_buffer(begin, end);
Expand All @@ -138,7 +168,7 @@ void PixelThresholdClusterizer::clusterizeDetUnitT( const T & input,

// Check if the cluster is above threshold
// (TO DO: one is signed, other unsigned, gcc warns...)
if ( cluster.charge() >= theClusterThreshold)
if ( cluster.charge() >= clusterThreshold)
Copy link
Contributor

Choose a reason for hiding this comment

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

actually given that cluster.charge() is a int please use a int for clusterThreshold!
moving to int the computation of cluster.charge() was a major speed-up improvement

{
// std::cout << "putting in this cluster " << i << " " << cluster.charge() << " " << cluster.pixelADC().size() << endl;
// sort by row (x)
Expand Down Expand Up @@ -202,10 +232,14 @@ void PixelThresholdClusterizer::copy_to_buffer( DigiIterator begin, DigiIterator
#endif
int electron[end-begin];
memset(electron, 0, sizeof(electron));
int layer = (DetId(detid_).subdetId()==1) ? PXBDetId(detid_).layer() : 0;
if ( doMissCalibrate ) {
(*theSiPixelGainCalibrationService_).calibrate(detid_,begin,end,theConversionFactor, theOffset,electron);
if (layer==1) {
(*theSiPixelGainCalibrationService_).calibrate(detid_,begin,end,theConversionFactor_L1, theOffset_L1,electron);
} else {
(*theSiPixelGainCalibrationService_).calibrate(detid_,begin,end,theConversionFactor, theOffset, electron);
}
} else {
int layer = (DetId(detid_).subdetId()==1) ? PXBDetId(detid_).layer() : 0;
int i=0;
for(DigiIterator di = begin; di != end; ++di) {
auto adc = di->adc();
Expand Down Expand Up @@ -315,7 +349,11 @@ int PixelThresholdClusterizer::calibrate(int adc, int col, int row)
//const float p3 = 113.0;
//float vcal = ( atanh( (adc-p3)/p2) + p1)/p0;

electrons = int( vcal * theConversionFactor + theOffset);
if (layer==1) {
electrons = int( vcal * theConversionFactor_L1 + theOffset_L1);
} else {
electrons = int( vcal * theConversionFactor + theOffset);
}
}
}
else
Expand Down Expand Up @@ -427,6 +465,11 @@ PixelThresholdClusterizer::make_cluster( const SiPixelCluster::PixelPos& pix,

if (dead_flag && doSplitClusters)
{
// Set separate cluster threshold for L1 (needed for phase1)
auto clusterThreshold = theClusterThreshold;
int layer = (DetId(detid_).subdetId()==1) ? PXBDetId(detid_).layer() : 0;
if (layer==1) clusterThreshold = theClusterThreshold_L1;

//Set the first cluster equal to the existing cluster.
SiPixelCluster first_cluster = cluster;
bool have_second_cluster = false;
Expand All @@ -440,8 +483,8 @@ PixelThresholdClusterizer::make_cluster( const SiPixelCluster::PixelPos& pix,
SiPixelCluster second_cluster = make_cluster(deadpix, output);

//If both clusters would normally have been found by the clusterizer, put them into output
if ( second_cluster.charge() >= theClusterThreshold &&
first_cluster.charge() >= theClusterThreshold )
if ( second_cluster.charge() >= clusterThreshold &&
first_cluster.charge() >= clusterThreshold )
{
output.push_back( second_cluster );
have_second_cluster = true;
Expand All @@ -460,7 +503,7 @@ PixelThresholdClusterizer::make_cluster( const SiPixelCluster::PixelPos& pix,
}

//Remember to also add the first cluster if we added the second one.
if ( first_cluster.charge() >= theClusterThreshold && have_second_cluster)
if ( first_cluster.charge() >= clusterThreshold && have_second_cluster)
{
output.push_back( first_cluster );
std::push_heap(output.begin(),output.end(),[](SiPixelCluster const & cl1,SiPixelCluster const & cl2) { return cl1.minPixelRow() < cl2.minPixelRow();});
Expand Down
Expand Up @@ -50,6 +50,9 @@
// Parameter Set:
#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"

#include <vector>


Expand All @@ -69,6 +72,8 @@ class dso_hidden PixelThresholdClusterizer final : public PixelClusterizerBase {
const std::vector<short>& badChannels,
edmNew::DetSetVector<SiPixelCluster>::FastFiller& output) { clusterizeDetUnitT(input, pixDet, badChannels, output); }

static void fillDescriptions(edm::ConfigurationDescriptions & descriptions);

private:

template<typename T>
Expand All @@ -88,11 +93,14 @@ class dso_hidden PixelThresholdClusterizer final : public PixelClusterizerBase {
float theSeedThresholdInNoiseUnits; // Pixel cluster seed in units of noise
float theClusterThresholdInNoiseUnits; // Cluster threshold in units of noise

const int thePixelThreshold; // Pixel threshold in electrons
const int theSeedThreshold; // Seed threshold in electrons
const float theClusterThreshold; // Cluster threshold in electrons
const int theConversionFactor; // adc to electron conversion factor
const int theOffset; // adc to electron conversion offset
const int thePixelThreshold; // Pixel threshold in electrons
const int theSeedThreshold; // Seed threshold in electrons
const int theClusterThreshold; // Cluster threshold in electrons
const int theClusterThreshold_L1; // Cluster threshold in electrons for Layer 1
const int theConversionFactor; // adc to electron conversion factor
const int theConversionFactor_L1; // adc to electron conversion factor for Layer 1
const int theOffset; // adc to electron conversion offset
const int theOffset_L1; // adc to electron conversion offset for Layer 1

const int theStackADC_; // The maximum ADC count for the stack layers
const int theFirstStack_; // The index of the first stack layer
Expand Down
Expand Up @@ -9,20 +9,33 @@
ChannelThreshold = cms.int32(1000),
MissCalibrate = cms.untracked.bool(True),
SplitClusters = cms.bool(False),
VCaltoElectronGain = cms.int32(65),
VCaltoElectronOffset = cms.int32(-414),
VCaltoElectronGain = cms.int32(65),
VCaltoElectronGain_L1 = cms.int32(65),
VCaltoElectronOffset = cms.int32(-414),
VCaltoElectronOffset_L1 = cms.int32(-414),
# **************************************
# **** payLoadType Options ****
# **** HLT - column granularity ****
# **** Offline - gain:col/ped:pix ****
# **************************************
payloadType = cms.string('Offline'),
SeedThreshold = cms.int32(1000),
ClusterThreshold = cms.double(4000.0),
ClusterThreshold = cms.int32(4000),
ClusterThreshold_L1 = cms.int32(4000),
# **************************************
maxNumberOfClusters = cms.int32(-1), # -1 means no limit.
)

# phase1 pixel
from Configuration.Eras.Modifier_phase1Pixel_cff import phase1Pixel
phase1Pixel.toModify(siPixelClusters,
VCaltoElectronGain = cms.int32(47), # L2-4: 47 +- 4.7
VCaltoElectronGain_L1 = cms.int32(50), # L1: 49.6 +- 2.6
VCaltoElectronOffset = cms.int32(-60), # L2-4: -60 +- 130
VCaltoElectronOffset_L1 = cms.int32(-670), # L1: -670 +- 220
ChannelThreshold = cms.int32(250),
ClusterThreshold_L1 = cms.int32(2000)
)

# Need these until phase2 pixel templates are used
from Configuration.Eras.Modifier_phase2_tracker_cff import phase2_tracker
Expand Down
7 changes: 7 additions & 0 deletions SimGeneral/MixingModule/python/SiPixelSimParameters_cfi.py
Expand Up @@ -29,6 +29,11 @@ def _modifyPixelDigitizerForPhase1Pixel( digitizer ) :
digitizer.BPix_SignalResponse_p1 = cms.double(0.711)
digitizer.BPix_SignalResponse_p2 = cms.double(203.)
digitizer.BPix_SignalResponse_p3 = cms.double(148.)
# gains and offsets are ints in the Clusterizer, so round to the same value
digitizer.ElectronsPerVcal = cms.double(47) # L2-4: 47 +- 4.7
digitizer.ElectronsPerVcal_L1 = cms.double(50) # L1: 49.6 +- 2.6
digitizer.ElectronsPerVcal_Offset = cms.double(-60) # L2-4: -60 +- 130
digitizer.ElectronsPerVcal_L1_Offset = cms.double(-670) # L1: -670 +- 220


SiPixelSimBlock = cms.PSet(
Expand Down Expand Up @@ -59,7 +64,9 @@ def _modifyPixelDigitizerForPhase1Pixel( digitizer ) :
BPix_SignalResponse_p2 = cms.double(97.4),
BPix_SignalResponse_p3 = cms.double(126.5),
ElectronsPerVcal = cms.double(65.5),
ElectronsPerVcal_L1 = cms.double(65.5),
ElectronsPerVcal_Offset = cms.double(-414.0),
ElectronsPerVcal_L1_Offset = cms.double(-414.0),
ElectronPerAdc = cms.double(135.0),
TofUpperCut = cms.double(12.5),
AdcFullScale = cms.int32(255),
Expand Down
13 changes: 9 additions & 4 deletions SimTracker/SiPixelDigitizer/plugins/SiPixelDigitizerAlgorithm.cc
Expand Up @@ -184,6 +184,8 @@ SiPixelDigitizerAlgorithm::SiPixelDigitizerAlgorithm(const edm::ParameterSet& co
// electrons to VCAL conversion needed in misscalibrate()
electronsPerVCAL(conf.getParameter<double>("ElectronsPerVcal")),
electronsPerVCAL_Offset(conf.getParameter<double>("ElectronsPerVcal_Offset")),
electronsPerVCAL_L1(conf.exists("ElectronsPerVcal_L1")?conf.getParameter<double>("ElectronsPerVcal_L1"):electronsPerVCAL),
electronsPerVCAL_L1_Offset(conf.exists("ElectronsPerVcal_L1_Offset")?conf.getParameter<double>("ElectronsPerVcal_L1_Offset"):electronsPerVCAL_Offset),

//theTofCut 12.5, cut in particle TOD +/- 12.5ns
//theTofCut(conf.getUntrackedParameter<double>("TofCut",12.5)),
Expand Down Expand Up @@ -1287,7 +1289,7 @@ void SiPixelDigitizerAlgorithm::make_digis(float thePixelThresholdInE,
if(doMissCalibrate) {
int row = ip.first; // X in row
int col = ip.second; // Y is in col
adc = int(missCalibrate(detID, pixdet, col, row, signalInElectrons)); //full misscalib.
adc = int(missCalibrate(detID, tTopo, pixdet, col, row, signalInElectrons)); //full misscalib.
} else { // Just do a simple electron->adc conversion
adc = int( signalInElectrons / theElectronPerADC ); // calibrate gain
}
Expand Down Expand Up @@ -1662,7 +1664,7 @@ float SiPixelDigitizerAlgorithm::pixel_aging(const PixelAging& aging,
//float offset = RandGaussQ::shoot(0.,theOffsetSmearing);
//float newAmp = amp * gain + offset;
// More complex misscalibration
float SiPixelDigitizerAlgorithm::missCalibrate(uint32_t detID, const PixelGeomDetUnit* pixdet, int col,int row,
float SiPixelDigitizerAlgorithm::missCalibrate(uint32_t detID, const TrackerTopology *tTopo, const PixelGeomDetUnit* pixdet, int col,int row,
const float signalInElectrons) const {
// Central values
//const float p0=0.00352, p1=0.868, p2=112., p3=113.; // pix(0,0,0)
Expand Down Expand Up @@ -1693,13 +1695,16 @@ float SiPixelDigitizerAlgorithm::missCalibrate(uint32_t detID, const PixelGeomDe
throw cms::Exception("NotAPixelGeomDetUnit") << "Not a pixel geomdet unit" << detID;
}

// const float electronsPerVCAL = 65.5; // our present VCAL calibration (feb 2009)
// const float electronsPerVCAL_Offset = -414.0; // our present VCAL calibration (feb 2009)
float newAmp = 0.; //Modified signal

// Convert electrons to VCAL units
float signal = (signalInElectrons-electronsPerVCAL_Offset)/electronsPerVCAL;

// New gains/offsets are needed for phase1 L1
int layer = 0;
if (DetId(detID).subdetId()==1) layer = tTopo->pxbLayer(detID);
if (layer==1) signal = (signalInElectrons-electronsPerVCAL_L1_Offset)/electronsPerVCAL_L1;

// Simulate the analog response with fixed parametrization
newAmp = p3 + p2 * tanh(p0*signal - p1);

Expand Down
Expand Up @@ -322,8 +322,10 @@ class SiPixelDigitizerAlgorithm {
const double theThresholdSmearing_BPix;
const double theThresholdSmearing_BPix_L1;

const double electronsPerVCAL; // for electrons - VCAL conversion
const double electronsPerVCAL_Offset; // in misscalibrate()
const float electronsPerVCAL; // for electrons - VCAL conversion
const float electronsPerVCAL_Offset; // in misscalibrate()
const float electronsPerVCAL_L1; // same for Layer 1
const float electronsPerVCAL_L1_Offset;// same for Layer 1

const float theTofLowerCut; // Cut on the particle TOF
const float theTofUpperCut; // Cut on the particle TOF
Expand Down Expand Up @@ -415,7 +417,7 @@ class SiPixelDigitizerAlgorithm {

// access to the gain calibration payloads in the db. Only gets initialized if check_dead_pixels_ is set to true.
const std::unique_ptr<SiPixelGainCalibrationOfflineSimService> theSiPixelGainCalibrationService_;
float missCalibrate(uint32_t detID, const PixelGeomDetUnit* pixdet, int col, int row, float amp) const;
float missCalibrate(uint32_t detID, const TrackerTopology *tTopo, const PixelGeomDetUnit* pixdet, int col, int row, float amp) const;
LocalVector DriftDirection(const PixelGeomDetUnit* pixdet,
const GlobalVector& bfield,
const DetId& detId) const;
Expand Down