Conversation
…s and made a seperate doc for 2d contour plots
|
Huma, I think this looks good. I think more documentation and comments would help to understand what each block of code is doing. Perhaps we can encapsulate some of the functions into sub-classes or at most other function to avoid large blocks. This is just style preference though - I can see a lot of great work here :) |
There was a problem hiding this comment.
Pull request overview
This PR implements enhancements to the Mu2e Source Calibration Analysis framework, focusing on three main areas: MC truth analysis for validation, an intelligent fitting retry mechanism to handle convergence failures, and 2D likelihood contour visualization for parameter correlation studies.
Changes:
- Added MC truth cut-and-count analysis module (
mcinfo.cc) to extract peak fractions from Monte Carlo truth histograms - Implemented adaptive retry logic in the fitter with three levels: initial fit, covariance-guided refit, and randomized parameter search
- Created standalone 2D contour plotting module (
2dcontour.cc) for visualizing parameter correlations and fit quality - Extended plotter to support asymmetric errors, LYSO/CsI crystal separation, and additional diagnostic plots
- Updated main executable to support MC mode and optional contour generation via command-line flags
Reviewed changes
Copilot reviewed 12 out of 12 changed files in this pull request and generated 17 comments.
Show a summary per file
| File | Description |
|---|---|
| SourceCalib/src/mcinfo.cc | New file implementing MC truth peak extraction and summary statistics |
| SourceCalib/src/SourceFitter.cc | Major refactor: added retry logic, asymmetric errors, improved error handling, message suppression |
| SourceCalib/src/2dcontour.cc | New file implementing 2D likelihood/chi2 contour scanning with auto-ranging |
| SourceCalib/src/SourcePlotter.cc | Enhanced plots: asymmetric error bars, CsI/LYSO splitting, odd/even comparisons, 45° projections |
| SourceCalib/src/MakeAnalysisTree_main.cc | Added MC mode block, improved flag parsing, convergence summary reporting |
| SourceCalib/inc/mcinfo.hh | Header for MC truth analysis class |
| SourceCalib/inc/SourceFitter.hh | Added static members for tracking retry statistics |
| SourceCalib/inc/SourcePlotter.hh | Added includes for new plot types |
| SourceCalib/inc/2dcontour.hh | Header for contour plotting function |
| SourceCalib/src/classes.h | Added includes for new modules |
| SourceCalib/src/classes_def.xml | Registered mcinfo class |
| SourceCalib/README.md | Comprehensive documentation update with usage examples |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| //add to counter to calculae avg of all values | ||
| covarAcc.count++; | ||
| covarAcc.sumPeak += fpeak; | ||
| covarAcc.sumWidth += fsigma; | ||
| covarAcc.sumAlpha += fcbalphaparam; | ||
| covarAcc.sumBeta += combetaparam; | ||
| covarAcc.sumEvtFull += fr_fullparam*nEvents; | ||
| covarAcc.sumEvtFst += fr_frstparam*nEvents; | ||
| covarAcc.sumEvtScd += fr_scndparam*nEvents; | ||
| covarAcc.sumEvtBkg += fr_bkgparam*nEvents; |
There was a problem hiding this comment.
The covariance accumulator is updated using uninitialized variables. Lines 224-231 use fpeak, fsigma, fcbalphaparam, combetaparam, fr_fullparam, fr_frstparam, fr_scndparam, and fr_bkgparam before they are assigned values. These parameters are only set much later (lines 355-364), after the fit has completed and the plot has been drawn. This means the accumulator is collecting garbage values on the first fit, which will corrupt the average values used for retry fits.
| //add to counter to calculae avg of all values | |
| covarAcc.count++; | |
| covarAcc.sumPeak += fpeak; | |
| covarAcc.sumWidth += fsigma; | |
| covarAcc.sumAlpha += fcbalphaparam; | |
| covarAcc.sumBeta += combetaparam; | |
| covarAcc.sumEvtFull += fr_fullparam*nEvents; | |
| covarAcc.sumEvtFst += fr_frstparam*nEvents; | |
| covarAcc.sumEvtScd += fr_scndparam*nEvents; | |
| covarAcc.sumEvtBkg += fr_bkgparam*nEvents; | |
| //add to counter to calculate avg of all values | |
| covarAcc.count++; | |
| covarAcc.sumPeak += fpeak; | |
| covarAcc.sumWidth += fsigma; | |
| // For parameters that are not initialized in run_one_fit(), avoid | |
| // reading them on the very first fit to prevent using undefined values. | |
| if (covarAcc.count > 1) { | |
| covarAcc.sumAlpha += fcbalphaparam; | |
| covarAcc.sumBeta += combetaparam; | |
| covarAcc.sumEvtFull += fr_fullparam*nEvents; | |
| covarAcc.sumEvtFst += fr_frstparam*nEvents; | |
| covarAcc.sumEvtScd += fr_scndparam*nEvents; | |
| covarAcc.sumEvtBkg += fr_bkgparam*nEvents; | |
| } |
| truthtable->cd(); | ||
| trueinfo->Write(); // Write the tree explicitly | ||
| truthtable->Write(); // Write any histograms created by FinalizeMCSummary | ||
| truthtable->Close(); |
There was a problem hiding this comment.
Memory leak: The truthtable TFile pointer created on line 109 is closed but never deleted. Add delete truthtable; after truthtable->Close(); on line 163.
| truthtable->Close(); | |
| truthtable->Close(); | |
| delete truthtable; |
| **1. Standard Analysis (Chi2 Fit)** | ||
| Accumulate and fit events for crystals 0 to 674 on disk 0 using Chi-Square minimization: | ||
|
|
||
| ./build/al9-prof-e28-p056/CaloCalibration/bin/MakeAnalysisTree 0 674 "chi2" 0 |
There was a problem hiding this comment.
The example command on line 88 is missing backticks for proper markdown code formatting, unlike the other examples. It should be wrapped in triple backticks like the other examples for consistency.
| // 1) First attempt: start from init values | ||
| // ----------------------- | ||
| bool first_ok = run_one_fit(); //run with init guesses | ||
| //add to counter to calculae avg of all values |
There was a problem hiding this comment.
Typo in comment: "calculae" should be "calculate".
| //add to counter to calculae avg of all values | |
| //add to counter to calculate avg of all values |
| errbarhigh = mparam*(peakerrorhigh/fpeak); | ||
| errbarlo = mparam*(peakerrorlo/fpeak); |
There was a problem hiding this comment.
Potential division by zero: If fpeak is zero (which could happen in pathological fit failures), lines 369-370 will divide by zero when computing errbarhigh and errbarlo. Add a check to ensure fpeak != 0 before performing these divisions.
| errbarhigh = mparam*(peakerrorhigh/fpeak); | |
| errbarlo = mparam*(peakerrorlo/fpeak); | |
| if (fpeak != 0.0) { | |
| errbarhigh = mparam * (peakerrorhigh / fpeak); | |
| errbarlo = mparam * (peakerrorlo / fpeak); | |
| } else { | |
| // Pathological case: peak height is zero; cannot compute scaled errors. | |
| errbarhigh = 0.0; | |
| errbarlo = 0.0; | |
| } |
| // Width of the first bin (in keV if the axis is in keV) | ||
| double binWidth = h_spec->GetXaxis()->GetBinWidth(500); | ||
| std::cout << "Bin width = " << binWidth << " keV" << std::endl; | ||
| //main function that assignes the values to the params |
There was a problem hiding this comment.
Typo in comment: "assignes" should be "assigns".
| //main function that assignes the values to the params | |
| //main function that assigns the values to the params |
| // Lambda: randomize all parameters | ||
| auto randomize_all_parameters = [&]() { | ||
| auto randomDouble = [](double min, double max) { | ||
| return min + (max - min) * ((double)rand() / RAND_MAX); |
There was a problem hiding this comment.
Using rand() for random number generation is not recommended for modern C++. The header <random> is included but not used. Consider using std::mt19937 and std::uniform_real_distribution for better randomness and thread safety, which is especially important since this code may be run in parallel contexts.
| return min + (max - min) * ((double)rand() / RAND_MAX); | |
| static thread_local std::mt19937 gen{std::random_device{}()}; | |
| std::uniform_real_distribution<double> dist(min, max); | |
| return dist(gen); |
| std::cout << " peak std dev " << Peaks_stddev << std::endl; | ||
| std::cout << " ------------ " << std::endl; | ||
|
|
||
| //----create loop to count and print how many crys/sipms have pear error larger than 2 stdev---- |
There was a problem hiding this comment.
Typo in comment: "pear" should be "peak".
| //----create loop to count and print how many crys/sipms have pear error larger than 2 stdev---- | |
| //----create loop to count and print how many crys/sipms have peak error larger than 2 stdev---- |
| }(); | ||
| void SourceFitter::FitCrystal(TH1F* h_spec, TString opt, int crystalNo, TTree *covar, Int_t &nEvents,Int_t &convergencestatus, Float_t &fpeak, Float_t &peakerrorhigh,Float_t &peakerrorlo, Float_t &fsigma,Float_t &widtherrorhigh,Float_t &widtherrorlo, Float_t &chiSq, Float_t &fstpeak, Float_t &scdpeak,Float_t &fcbalphaparam,Float_t &fcbndegparam,Float_t &comCnstparam, Float_t &combetaparam, Float_t &fr_fullparam, Float_t &fr_frstparam,Float_t &fr_scndparam,Float_t &crystalNoparam,Float_t &fr_bkgparam, Float_t &pval,Float_t &h_means,Float_t &h_stddevs, Float_t &unreducedchi2,Float_t &fval,Float_t &mparam,Float_t &etaparam,Int_t &ndof,bool contour, Float_t &errbarhigh, Float_t &errbarlo,Float_t &evtfullerrorhigh,Float_t &evtfullerrorlo,Float_t &eventsFull,Float_t &Esparam ){ | ||
|
|
||
| // set stlye optionsr |
There was a problem hiding this comment.
Typo in comment: "stlye optionsr" should be "style options".
| // set stlye optionsr | |
| // set style options |
Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
…input values from command line
Updated code logic to re-try fits for failing sipms, added cut and count analysis for mc truth, created a separate file for 2d contour and updated the readme