Skip to content

Commit

Permalink
new parametrization and geometry v10
Browse files Browse the repository at this point in the history
  • Loading branch information
@deguio
deguio committed May 16, 2019
1 parent e6e9f7b commit 625bdda
Show file tree
Hide file tree
Showing 8 changed files with 190 additions and 94 deletions.
6 changes: 5 additions & 1 deletion RecoLocalCalo/HGCalRecProducers/interface/HGCalImagingAlgo.h
View file
Expand Up @@ -114,7 +114,11 @@ static void fillPSetDescription(edm::ParameterSetDescription& iDesc) {
descNestedNoises.add<std::vector<double> >("values", {});
iDesc.add<edm::ParameterSetDescription>("noises", descNestedNoises);
edm::ParameterSetDescription descNestedNoiseMIP;
descNestedNoiseMIP.add<double>("value", 0 );
descNestedNoiseMIP.add<bool>("scaleByDose", false );
iDesc.add<edm::ParameterSetDescription>("scaleByDose", descNestedNoiseMIP);
descNestedNoiseMIP.add<std::string>("doseMap", "" );
iDesc.add<edm::ParameterSetDescription>("doseMap", descNestedNoiseMIP);
descNestedNoiseMIP.add<double>("noise_MIP", 1./100. );
iDesc.add<edm::ParameterSetDescription>("noiseMip", descNestedNoiseMIP);
}

Expand Down
16 changes: 0 additions & 16 deletions SimCalorimetry/HGCalSimProducers/data/doseParams_3000fb.txt
View file

This file was deleted.

14 changes: 14 additions & 0 deletions SimCalorimetry/HGCalSimProducers/data/doseParams_3000fb_fluka-3.5.15.9.txt
View file
@@ -0,0 +1,14 @@
9 3.985908 -0.0233164 0.000479769 -8.90993e-06 5.74269e-08 16.460299 -0.0282926 0.000176 -1.24226e-06 3.33661e-09
10 3.865145 -0.0157737 1.4102e-05 -2.59944e-08 5.18596e-09 16.074041 -0.0142934 -2.3325e-05 -9.91417e-08 1.05114e-09
11 3.797408 -0.0110597 -0.000150187 1.70552e-06 -9.14724e-10 15.856707 -0.00698478 -0.000129121 5.11442e-07 -1.75663e-10
12 3.851879 -0.00268387 -0.000502753 6.91762e-06 -2.62058e-08 15.514695 0.00534509 -0.000304267 1.51454e-06 -2.17969e-09
13 3.637403 -0.0102584 -0.000280915 4.39397e-06 -1.67187e-08 15.235543 0.0140316 -0.000430826 2.25819e-06 -3.69355e-09
14 3.502875 -0.00930028 -0.000349823 5.39975e-06 -2.10104e-08 15.364650 0.00885852 -0.000386321 2.09248e-06 -3.47966e-09
15 3.408253 -0.00713129 -0.000358892 4.95141e-06 -1.79598e-08 15.591748 0.00056066 -0.000302279 1.70913e-06 -2.84797e-09
16 3.308403 -0.010831 -0.000302115 4.5388e-06 -1.67407e-08 15.353311 0.00786759 -0.000395923 2.15057e-06 -3.55091e-09
17 3.274672 -0.00658889 -0.000350272 4.3462e-06 -1.44797e-08 15.555453 0.00110278 -0.000332765 1.8694e-06 -3.08188e-09
18 3.083276 -0.0102271 -0.00028511 3.5795e-06 -1.08933e-08 15.822784 -0.00815451 -0.000231267 1.35734e-06 -2.15243e-09
19 3.125504 -0.0073409 -0.000282136 3.16049e-06 -9.34773e-09 15.909287 -0.0112805 -0.000205724 1.24919e-06 -1.99359e-09
20 2.995494 -0.0193198 -6.97825e-06 6.81663e-07 -1.76918e-09 16.115456 -0.0185292 -0.000131819 9.26615e-07 -1.49273e-09
21 2.973583 -0.00764613 -0.000184954 1.73145e-06 -4.03636e-09 16.356612 -0.0288606 -1.22903e-05 3.72635e-07 -6.06651e-10
22 2.821734 -0.0158869 1.38538e-06 3.83847e-07 -1.34869e-09 16.758791 -0.0481384 0.000220859 -6.4673e-07 8.8165e-10
19 changes: 10 additions & 9 deletions SimCalorimetry/HGCalSimProducers/interface/HGCHEbackDigitizer.h
View file
Expand Up @@ -11,8 +11,9 @@ class HGCHEbackSignalScaler
public:

struct DoseParameters {
DoseParameters(): a_(0.), b_(0.), c_(0.), d_(0.), e_(0.), f_(0.) {}
float a_, b_, c_, d_, e_, f_;
DoseParameters(): a_(0.), b_(0.), c_(0.), d_(0.), e_(0.),
f_(0.), g_(0.), h_(0.), i_(0.), j_(0.) {}
float a_, b_, c_, d_, e_, f_, g_, h_, i_, j_;
};

HGCHEbackSignalScaler() {};
Expand All @@ -21,19 +22,19 @@ class HGCHEbackSignalScaler
void setGeometry(const CaloSubdetectorGeometry*);
void setDoseMap(const std::string&);

float scaleByArea(const HGCScintillatorDetId&, const std::pair<double, double>&);
std::pair<float, float> scaleByDose(const HGCScintillatorDetId&, const std::pair<double, double>&);
double getDoseValue(const int, const std::pair<double, double>&);
double getFluenceValue(const int, const std::pair<double, double>&);
std::pair<double, double> computeRadius(const HGCScintillatorDetId&);
float scaleByArea(const HGCScintillatorDetId&, const std::array<double, 8>&);
std::pair<float, float> scaleByDose(const HGCScintillatorDetId&, const std::array<double, 8>&);
double getDoseValue(const int, const std::array<double, 8>&);
double getFluenceValue(const int, const std::array<double, 8>&);
std::array<double, 8> computeRadius(const HGCScintillatorDetId&);

private:
std::map<int, DoseParameters> readDosePars(const std::string&);

const HGCalGeometry* hgcalGeom_;
std::map<int, DoseParameters> doseMap_;
static constexpr double radToKrad_ = 1/1000.;
const float refEdge_ = 0.03; //3 cm
static constexpr double greyToKrad_ = 0.1;
const float refEdge_ = 3; //3 cm

bool verbose_ = false;

Expand Down
6 changes: 3 additions & 3 deletions SimCalorimetry/HGCalSimProducers/python/hgcalDigitizer_cfi.py
View file
Expand Up @@ -180,9 +180,9 @@
calibDigis = cms.bool(True),
keV2MIP = cms.double(1./500.0),
doTimeSamples = cms.bool(False),
nPEperMIP = cms.double(18.0), #conservative
nPEperMIP = cms.double(21.0),
nTotalPE = cms.double(7500),
xTalk = cms.double(0.2),
xTalk = cms.double(0.01),
sdPixels = cms.double(1e-6), # this is additional photostatistics noise (as implemented), not sure why it's here...
feCfg = cms.PSet(
# 0 only ADC, 1 ADC with pulse shape, 2 ADC+TDC with pulse shape
Expand Down Expand Up @@ -281,7 +281,7 @@ def HGCal_setEndOfLifeNoise(process):
)
process.HGCAL_noise_heback = cms.PSet(
scaleByDose = cms.bool(True),
doseMap = cms.string("SimCalorimetry/HGCalSimProducers/data/doseParams_3000fb.txt"),
doseMap = cms.string("SimCalorimetry/HGCalSimProducers/data/doseParams_3000fb_fluka-3.5.15.9.txt"),
noise_MIP = cms.double(1./5.) #uses noise map
)
process.HGCAL_noises = cms.PSet(
Expand Down
43 changes: 27 additions & 16 deletions SimCalorimetry/HGCalSimProducers/src/HGCHEbackDigitizer.cc
View file
Expand Up @@ -44,26 +44,26 @@ std::map<int, HGCHEbackSignalScaler::DoseParameters> HGCHEbackSignalScaler::read

//space-separated
std::stringstream linestream(line);
linestream >> layer >> dosePars.a_ >> dosePars.b_ >> dosePars.c_ >> dosePars.d_ >> dosePars.e_ >> dosePars.f_;
linestream >> layer >> dosePars.a_ >> dosePars.b_ >> dosePars.c_ >> dosePars.d_ >> dosePars.e_ >> dosePars.f_ >> dosePars.g_ >> dosePars.h_ >> dosePars.i_ >> dosePars.j_;

result[layer] = dosePars;
}
return result;
}

double HGCHEbackSignalScaler::getDoseValue(const int layer, const std::pair<double, double>& radius)
double HGCHEbackSignalScaler::getDoseValue(const int layer, const std::array<double, 8>& radius)
{
double cellDose = std::pow(10, doseMap_[layer].a_ + doseMap_[layer].b_*radius.first + doseMap_[layer].c_*radius.second); //dose in rad
return cellDose * radToKrad_; //convert to kRad
double cellDose = std::pow(10, doseMap_[layer].a_ + doseMap_[layer].b_*radius[4] + doseMap_[layer].c_*radius[5] + doseMap_[layer].d_*radius[6] + doseMap_[layer].e_*radius[7]); //dose in grey
return cellDose * greyToKrad_; //convert to kRad
}

double HGCHEbackSignalScaler::getFluenceValue(const int layer, const std::pair<double, double>& radius)
double HGCHEbackSignalScaler::getFluenceValue(const int layer, const std::array<double, 8>& radius)
{
double cellFluence = std::pow(10, doseMap_[layer].d_ + doseMap_[layer].e_*radius.first + doseMap_[layer].f_*radius.second); //dose in rad
double cellFluence = std::pow(10, doseMap_[layer].f_ + doseMap_[layer].g_*radius[0] + doseMap_[layer].h_*radius[1] + doseMap_[layer].i_*radius[2] + doseMap_[layer].j_*radius[3]); //dose in grey
return cellFluence;
}

std::pair<float, float> HGCHEbackSignalScaler::scaleByDose(const HGCScintillatorDetId& cellId, const std::pair<double, double>& radius)
std::pair<float, float> HGCHEbackSignalScaler::scaleByDose(const HGCScintillatorDetId& cellId, const std::array<double, 8>& radius)
{
if(doseMap_.empty())
return std::make_pair(1., 0.);
Expand All @@ -75,7 +75,7 @@ std::pair<float, float> HGCHEbackSignalScaler::scaleByDose(const HGCScintillator

double cellFluence = getFluenceValue(layer, radius); //in 1-Mev-equivalent neutrons per cm2

constexpr double factor = 1. / (2*1e13);
constexpr double factor = 2. / (2*1e13); //SiPM area = 2mm^2
double noise = 2.18 * sqrt(cellFluence * factor);

if(verbose_)
Expand All @@ -97,9 +97,9 @@ std::pair<float, float> HGCHEbackSignalScaler::scaleByDose(const HGCScintillator
return std::make_pair(scaleFactor, noise);
}

float HGCHEbackSignalScaler::scaleByArea(const HGCScintillatorDetId& cellId, const std::pair<double, double>& radius)
float HGCHEbackSignalScaler::scaleByArea(const HGCScintillatorDetId& cellId, const std::array<double, 8>& radius)
{
float circ = 2 * M_PI * radius.first;
float circ = 2 * M_PI * radius[0];

float edge(refEdge_);
if(cellId.type() == 0)
Expand All @@ -117,22 +117,33 @@ float HGCHEbackSignalScaler::scaleByArea(const HGCScintillatorDetId& cellId, con

if(verbose_)
{
LogDebug("HGCHEbackSignalScaler") << "HGCHEbackSignalScaler::scaleByArea - Type, layer, edge, radius: "
LogDebug("HGCHEbackSignalScaler") << "HGCHEbackSignalScaler::scaleByArea - Type, layer, edge, radius, SF: "
<< cellId.type() << " "
<< cellId.layer() << " "
<< edge << " "
<< radius;
<< radius[0] << " "
<< scaleFactor << std::endl;
}

return scaleFactor;
}

std::pair<double, double> HGCHEbackSignalScaler::computeRadius(const HGCScintillatorDetId& cellId)
std::array<double, 8> HGCHEbackSignalScaler::computeRadius(const HGCScintillatorDetId& cellId)
{
GlobalPoint global = hgcalGeom_->getPosition(cellId);
double radius2 = (std::pow(global.x(), 2) + std::pow(global.y(), 2)) * 0.0001; //in m

double radius2 = std::pow(global.x(), 2) + std::pow(global.y(), 2); //in cm
double radius4 = std::pow(radius2, 2);
double radius = sqrt(radius2);
return std::make_pair(radius, radius2);
double radius3 = std::pow(radius, 3);

double radius_m100 = radius-100;
double radius_m100_2 = std::pow(radius_m100, 2);
double radius_m100_3 = std::pow(radius_m100, 3);
double radius_m100_4 = std::pow(radius_m100_2, 2);

std::array<double, 8> radii { {radius, radius2, radius3, radius4, radius_m100, radius_m100_2, radius_m100_3, radius_m100_4} };
return radii;
}


Expand Down Expand Up @@ -236,7 +247,7 @@ void HGCHEbackDigitizer::runRealisticDigitizer(std::unique_ptr<HGCalDigiCollecti

if(id.det() == DetId::HGCalHSc) //skip those geometries that have HE used as BH
{
std::pair<double, double> radius;
std::array<double, 8> radius;
if(scaleByArea_ or scaleByDose_)
radius = scal_.computeRadius(id);

Expand Down

0 comments on commit 625bdda

Please sign in to comment.