Skip to content

Commit

Permalink
Merge pull request #36266 from tvami/LANewMethod
Browse files Browse the repository at this point in the history 
[121X] Modify calculations for the Lorentz angle
  • Loading branch information
@cmsbuild
cmsbuild committed Nov 28, 2021
2 parents 357a595 + 17b8322 commit 6213a97
Show file tree
Hide file tree
Showing 2 changed files with 36 additions and 6 deletions.
40 changes: 35 additions & 5 deletions CalibTracker/SiPixelLorentzAngle/src/SiPixelLorentzAnglePCLHarvester.cc
View file
Expand Up @@ -44,6 +44,7 @@ class SiPixelLorentzAnglePCLHarvester : public DQMEDHarvester {
const double fitProbCut_;
const std::string recordName_;
std::unique_ptr<TF1> f1;
float width_;

SiPixelLorentzAngleCalibrationHistograms hists;
const SiPixelLorentzAngle* currentLorentzAngle;
Expand Down Expand Up @@ -106,6 +107,7 @@ void SiPixelLorentzAnglePCLHarvester::beginRun(const edm::Run& iRun, const edm::

for (auto det : geom->detsPXB()) {
const PixelGeomDetUnit* pixelDet = dynamic_cast<const PixelGeomDetUnit*>(det);
width_ = pixelDet->surface().bounds().thickness();
const auto& layer = tTopo->pxbLayer(pixelDet->geographicalId());
const auto& module = tTopo->pxbModule(pixelDet->geographicalId());
int i_index = module + (layer - 1) * hists.nModules_[layer - 1];
Expand Down Expand Up @@ -214,7 +216,8 @@ void SiPixelLorentzAnglePCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQMS
<< "p4" << "\t" << "e4" << "\t"
<< "p5" << "\t" << "e5" << "\t"
<< "chi2" << "\t" << "prob" << "\t"
<< "newDetId" << std::endl;
<< "newDetId" << "\t" << "tan(LA)" << "\t"
<< "Error(LA)" << std::endl;
// clang-format on

std::unique_ptr<SiPixelLorentzAngle> LorentzAngle = std::make_unique<SiPixelLorentzAngle>();
Expand All @@ -228,6 +231,7 @@ void SiPixelLorentzAnglePCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQMS
f1->SetParName(5, "quintic term");

double p1_simul_newmodule = 0.294044;
double half_width = width_ * 10000 / 2; // pixel half thickness in units of micro meter

for (int j = 0; j < (int)hists.BPixnewDetIds_.size(); j++) {
int new_index = j + 1 + hists.nModules_[hists.nlay - 1] + (hists.nlay - 1) * hists.nModules_[hists.nlay - 1];
Expand Down Expand Up @@ -268,12 +272,24 @@ void SiPixelLorentzAnglePCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQMS
double chi2 = f1->GetChisquare();
double prob = f1->GetProb();

double f1_halfwidth = p0 + p1 * half_width + p2 * pow(half_width, 2) + p3 * pow(half_width, 3) +
p4 * pow(half_width, 4) + p5 * pow(half_width, 5);

double f1_zerowidth = p0;

double tan_LA =
(f1_halfwidth - f1_zerowidth) / half_width; // tan_LA = (f1(x = half_width) - f1(x = 0)) / (half_width - 0)
double errsq_LA = (pow(e1, 2) + pow((half_width * e2), 2) + pow((half_width * half_width * e3), 2) +
pow((half_width * half_width * half_width * e4), 2) +
pow((half_width * half_width * half_width * half_width * e5), 2)); // Propagation of uncertainty
double error_LA = sqrt(errsq_LA);

edm::LogPrint("LorentzAngle") << std::setprecision(4) << hists.BPixnewModule_[j] << "\t" << hists.BPixnewLayer_[j]
<< "\t" << p0 << "\t" << e0 << "\t" << p1 << std::setprecision(3) << "\t" << e1
<< "\t" << e1 / p1 * 100. << "\t" << (p1 - p1_simul_newmodule) / e1 << "\t" << p2
<< "\t" << e2 << "\t" << p3 << "\t" << e3 << "\t" << p4 << "\t" << e4 << "\t" << p5
<< "\t" << e5 << "\t" << chi2 << "\t" << prob << "\t" << hists.BPixnewDetIds_[j]
<< std::endl;
<< "\t" << tan_LA << "\t" << error_LA << std::endl;
}

double p1_simul[hists.nlay][hists.nModules_[hists.nlay - 1]];
Expand Down Expand Up @@ -335,12 +351,26 @@ void SiPixelLorentzAnglePCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQMS
double chi2 = f1->GetChisquare();
double prob = f1->GetProb();

double f1_halfwidth = p0 + p1 * half_width + p2 * pow(half_width, 2) + p3 * pow(half_width, 3) +
p4 * pow(half_width, 4) + p5 * pow(half_width, 5);

double f1_zerowidth = p0;

double tan_LA =
(f1_halfwidth - f1_zerowidth) / half_width; // tan_LA = (f1(x = half_width) - f1(x = 0)) / (half_width - 0)
double errsq_LA =
(pow(e1, 2) + pow((half_width * e2), 2) + pow((half_width * half_width * e3), 2) +
pow((half_width * half_width * half_width * e4), 2) +
pow((half_width * half_width * half_width * half_width * e5), 2)); // Propagation of uncertainty
double error_LA = sqrt(errsq_LA);

edm::LogPrint("LorentzAngle") << std::setprecision(4) << i_module << "\t" << i_layer << "\t" << p0 << "\t" << e0
<< "\t" << p1 << std::setprecision(3) << "\t" << e1 << "\t" << e1 / p1 * 100.
<< "\t" << (p1 - p1_simul[i_layer - 1][i_module - 1]) / e1 << "\t" << p2 << "\t"
<< e2 << "\t" << p3 << "\t" << e3 << "\t" << p4 << "\t" << e4 << "\t" << p5 << "\t"
<< e5 << "\t" << chi2 << "\t" << prob << "\t"
<< "null" << std::endl;
<< "null"
<< "\t" << tan_LA << "\t" << error_LA << std::endl;

const auto& detIdsToFill = hists.detIdsList.at(i_index);

Expand All @@ -351,7 +381,7 @@ void SiPixelLorentzAnglePCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQMS
// if the fit quality is OK
if (prob > fitProbCut_) {
for (const auto& id : detIdsToFill) {
bPixLorentzAnglePerTesla_ = p1 / theMagField;
bPixLorentzAnglePerTesla_ = tan_LA / theMagField;
if (!LorentzAngle->putLorentzAngle(id, bPixLorentzAnglePerTesla_)) {
edm::LogError("SiPixelLorentzAnglePCLHarvester")
<< "[SiPixelLorentzAnglePCLHarvester::dqmEndRun] detid already exists" << std::endl;
Expand Down Expand Up @@ -412,7 +442,7 @@ void SiPixelLorentzAnglePCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQMS
edm::Service<cond::service::PoolDBOutputService> mydbservice;
if (mydbservice.isAvailable()) {
try {
mydbservice->writeOneIOV(LorentzAngle.get(), mydbservice->currentTime(), recordName_);
mydbservice->writeOneIOV(*LorentzAngle, mydbservice->currentTime(), recordName_);
} catch (const cond::Exception& er) {
edm::LogError("SiPixelLorentzAngleDB") << er.what() << std::endl;
} catch (const std::exception& er) {
Expand Down
2 changes: 1 addition & 1 deletion CalibTracker/SiPixelLorentzAngle/src/SiPixelLorentzAnglePCLWorker.cc
View file
Expand Up @@ -379,7 +379,7 @@ void SiPixelLorentzAnglePCLWorker::analyze(edm::Event const& iEvent, edm::EventS
const PixelTopology* topol = &(theGeomDet->specificTopology());

float ypitch_ = topol->pitch().second;
float width_ = 0.0285;
float width_ = theGeomDet->surface().bounds().thickness();

if (!topol)
continue;
Expand Down

0 comments on commit 6213a97

Please sign in to comment.