Skip to content

Modify calculations for the Lorentz angle #17

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

Merged
mmusich merged 2 commits into from
Nov 24, 2021
Merged

Modify calculations for the Lorentz angle #17

Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
dc0f559
Modify calculations for the Lorentz angle
wweiphy Nov 17, 2021
43fcd1e
Remove magic numbers
wweiphy Nov 19, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
37 changes: 33 additions & 4 deletions CalibTracker/SiPixelLorentzAngle/src/SiPixelLorentzAnglePCLHarvester.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -214,7 +214,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 +229,7 @@ void SiPixelLorentzAnglePCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQMS
f1->SetParName(5, "quintic term");

double p1_simul_newmodule = 0.294044;
double half_width = 0.0285 / 2 * 10000; // pixel width in units of micro meter
Copy link

@mmusich mmusich Nov 17, 2021

Choose a reason for hiding this comment

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

No magic numbers, please fetch from geometry

Copy link

@tvami tvami Nov 18, 2021

Choose a reason for hiding this comment

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

@wweiphy can you please do this?


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 +270,25 @@ 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)) *
pow(tan_LA, 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 +350,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)) *
pow(tan_LA, 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 +380,7 @@ void SiPixelLorentzAnglePCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQMS
// if the fit quality is OK
if (prob > fitProbCut_) {
Copy link

@mmusich mmusich Nov 17, 2021

Choose a reason for hiding this comment

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

I just want to draw your attention to the fact that right now fitProbCut_ is set to 0.5

cmssw/CalibTracker/SiPixelLorentzAngle/src/SiPixelLorentzAnglePCLHarvester.cc

Line 498 in dc0f559

desc.add<double>("fitProbCut", 0.5)->setComment("cut on fit chi2 probabiblity to accept measurement");

which is not tuned at all, it's just a random value I picked.

thinking in the long run for the multi-run harvesting applications, a part from the minimum number of events in the MRH configuration, I think it would be good to have some sort of internal logic from within the SiPixel LA harvester to veto the creation of the payload (as it is done for all the other Tracker PCL workflows)
At the moment it just checks that the single "module-group" fit probability is higher than this given cut and if that condition is not satisfied it copies the values from the current value stored in DB (not stopping the payload creation).
If a more refined logic is introduced the input can be balanced in a more refined way, since the MRH will retry aggregating more statistics until a valid payload is produced.

@wweiphy @tvami @ferencek @OzAmram

Copy link
Member

@ferencek ferencek Nov 22, 2021

Choose a reason for hiding this comment

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

We have a similar issue with the SiPixelQuality payloads where the old payloads gets appended when the occupancy threshold is not passed. So it looks like we need a similar fix for both types of payloads.

Concerning the fitProbCut, I guess it would make sense to have a look at the distribution and see if some other cut value would make more sense. However, the question of a more refined internal logic still remains. Once the required minimum number of events is reached, is it possible to continue collecting statistics until some other conditions is met (and only then the payload is created). I really don't know much about MRH and how complicated the logic can be.

Copy link

@mmusich mmusich Nov 22, 2021

Choose a reason for hiding this comment

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

However, the question of a more refined internal logic still remains. Once the required minimum number of events is reached, is it possible to continue collecting statistics until some other conditions is met (and only then the payload is created). I really don't know much about MRH and how complicated the logic can be.

the logic better be implemented on the cmssw side, where you have full control and NOT in the MRH backend.

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