Skip to content

GitLab

Commit d1e58991 authored by Jérémie Dudouet's avatar Jérémie Dudouet
Browse files

correct a bug when calibraing mono energetic sources

parent 6aa4388f
Pipeline #146443 passed with stage
in 10 minutes and 12 seconds
Hide whitespace changes
Inline Side-by-side
Showing with 124 additions and 120 deletions
src/root/gui/cubix/src/CXHist1DCalib.cpp
View file @ d1e58991
......@@ -275,7 +275,9 @@ void CXHist1DCalib::Calibrate()
vector < Fitted > FitResults = fRecalEnergy->GetFitResults();
if(fVerboseLevel->GetNumber()==0) {
cout<<"FitResults.size(): "<<FitResults.size()<<endl;
if(fVerboseLevel->GetNumber()==0 && fRecalEnergy->fCalibFunction) {
cout<< left << scientific << setprecision(6);
cout<< hist->GetName()<<": ";
cout << setw(14) << fRecalEnergy->fCalibFunction->GetParameter(0);
......@@ -283,14 +285,12 @@ void CXHist1DCalib::Calibrate()
cout<<endl;
}
double xmin=-1;
double xmax=-1;
int NGoodPeak=0;
for(size_t i=0 ; i<FitResults.size() ; i++)
{
for(size_t i=0 ; i<FitResults.size() ; i++) {
Fitted FitRes = FitResults[i];
if(!FitRes.good) continue;
......@@ -345,20 +345,24 @@ void CXHist1DCalib::Calibrate()
i += NSubPeaks-1;
}
hist->GetXaxis()->SetRangeUser(xmin-(xmax-xmin)*0.1,xmax+(xmax-xmin)*0.1);
if(xmin!=-1) hist->GetXaxis()->SetRangeUser(xmin-(xmax-xmin)*0.1,xmax+(xmax-xmin)*0.1);
fMainWindow->RefreshPads();
if(fCalibCanvas != nullptr) delete fCalibCanvas;
fCalibCanvas = new TCanvas;
fCalibCanvas->SetName("CalibrationResults");
fCalibCanvas->SetTitle("Calibration Results");
fCalibCanvas->Divide(1,2,0.0001,0.0001);
fCalibCanvas->cd(1);
fRecalEnergy->fCalibGraph->Draw("ap");
fRecalEnergy->fCalibFunction->Draw("same");
fCalibCanvas->cd(2);
fRecalEnergy->fResidueGraph->Draw("ape");
if(fRecalEnergy->fCalibFunction) {
if(fCalibCanvas != nullptr) delete fCalibCanvas;
fCalibCanvas = new TCanvas;
fCalibCanvas->SetName("CalibrationResults");
fCalibCanvas->SetTitle("Calibration Results");
fCalibCanvas->Divide(1,2,0.0001,0.0001);
fCalibCanvas->cd(1);
fRecalEnergy->fCalibGraph->Draw("ap");
fRecalEnergy->fCalibFunction->Draw("same");
fCalibCanvas->cd(2);
fRecalEnergy->fResidueGraph->Draw("ape");
fCalibCanvas->Update();
fCalibCanvas->Modified();
}
fMainWindow->GetCanvas()->cd();
}
......
src/root/tools/GwRecalEnergy.cpp
View file @ d1e58991
......@@ -139,6 +139,8 @@ void GwRecalEnergy::StartCalib()
else
np = PeakSearch2(specData, specFromDef, specToDef, specFWHMdef, specAMPLdef, specPeaks);
cout<<"yo :" <<np<<endl;
if(np) {
int fp = FitPeaks( Verbosity );
if(fp) {
......@@ -160,6 +162,7 @@ void GwRecalEnergy::StartCalib()
Peaks.clear();
if(Verbosity > 0)
cout<<Form("\n");
return;
}
if(fCalibOrder) {
......@@ -821,9 +824,9 @@ int GwRecalEnergy::xP_Next1(float *yDat, int yLen)
}
else { // transition to negative
if(widthP >= xP_nMinWidthP && // positive lobe wide enough
//widthP < 4*xP_nMinWidthP && // but not too wide
xP_1 > xP_0 && // and has a maximum
yMax > xP_minAmpl) { // which is large enough
//widthP < 4*xP_nMinWidthP && // but not too wide
xP_1 > xP_0 && // and has a maximum
yMax > xP_minAmpl) { // which is large enough
state = stateN; // initiate a negative lobe sequence
xP_3 = xP_4 = xP_5 = nn;
yMin = yInp;
......@@ -842,8 +845,8 @@ int GwRecalEnergy::xP_Next1(float *yDat, int yLen)
xP_5 = nn;
widthN++; // incr negative width and stay in this state
if(widthN > xP_nMinWidthN && // second positive lobe wide enough for a valid sequence
yInp > yMin && // and signal is beyond minimum
yMin < -xP_minAmpl) { // and minimum large enough
yInp > yMin && // and signal is beyond minimum
yMin < -xP_minAmpl) { // and minimum large enough
float cmp = abs(yMax/yMin);
if(cmp < 0.2f || cmp > 2.f)
state = stateUnknown;
......@@ -905,7 +908,7 @@ int GwRecalEnergy::xP_Next2(float *yDat, int yLen)
}
else { // transition to negative
if(widthP1 >= xP_nMinWidthP1 && // positive lobe wide enough
xP_1 > xP_0) { // and has a maximum
xP_1 > xP_0) { // and has a maximum
state = stateN; // initiate a negative lobe sequence
xP_3 = xP_4 = nn;
yMin = yInp;
......@@ -926,9 +929,9 @@ int GwRecalEnergy::xP_Next2(float *yDat, int yLen)
}
else { // transition to positive
if(widthN >= xP_nMinWidthN && // negative lobe wide enough
yMin < -xP_minAmpl && // its amplitude large enough
yMin < -1.1*yMax1 && // and minimum-to-maximum is acceptable
yMin > -4.5*yMax1) { // (exact value is -0.5*exp(1.5)) = -2.2408)
yMin < -xP_minAmpl && // its amplitude large enough
yMin < -1.1*yMax1 && // and minimum-to-maximum is acceptable
yMin > -4.5*yMax1) { // (exact value is -0.5*exp(1.5)) = -2.2408)
state = stateP2; // initiate the second positive lobe sequence
xP_5 = xP_6 = nn;
yMax2 = yInp;
......@@ -951,7 +954,7 @@ int GwRecalEnergy::xP_Next2(float *yDat, int yLen)
xP_6 = nn;
widthP2++; // incr positive width and stay in this state
if(widthP2 > xP_nMinWidthP2 && // second positive lobe wide enough for a valid sequence
yInp < yMax2) { // and signal is beyond maximum
yInp < yMax2) { // and signal is beyond maximum
return nn; // this is a good trigger
}
}
......@@ -1050,6 +1053,10 @@ Int_t GwRecalEnergy::EROOTCalibration()
{
vector <Fitted>::iterator Iter;
double bestSlope = 1.;
Float_t Xmin = 0;
Float_t Xmax = 0.;
if(Energies.size() == 1) {
// select the peak with the largest area
int pn_max = 0;
......@@ -1058,118 +1065,112 @@ Int_t GwRecalEnergy::EROOTCalibration()
pn_max = nn;
}
}
double slope_ = Energies[0] / Peaks[pn_max].posi;
bestSlope = Energies[0] / Peaks[pn_max].posi;
Peaks[pn_max].good = true;
Peaks[pn_max].erefindex = 0;
return 0;
Xmin = Peaks[pn_max].posi-100;
Xmax = Peaks[pn_max].posi+100;
}
if( Peaks.size() < 2 || Energies.size() < 2) {
else if( Peaks.size() < 2 || Energies.size() < 2){
if(Verbosity> 1) cout << "# Number of peaks (" << Peaks.size() << ") or number of energies (" << Energies.size() << ") too small" << endl;
return 0;
}
else {
#if 0
// all peaks used as pivot
int npmax = Peaks.size();
for(int np = 0; np < npmax; np++)
Peaks[np].good = true;
#else
// only the peak with largest area used as pivot
sort( Peaks.begin( ), Peaks.end( ), largerAmplitude() );
int npmax = min( min(Peaks.size(), size_t(4)), Energies.size() ); // limit the number of peaks to min(np,ne,4)
for(int np = 0; np < npmax; np++)
Peaks[np].good = true;
for(size_t np = npmax; np < Peaks.size(); np++)
Peaks[np].good = false;
#endif
// (re)order the peaks
sort( Peaks.begin( ), Peaks.end( ), smallerPosi() );
if(Verbosity> 2) {
cout<<Form("# Sorted Peak Positions = (");
for ( Iter = Peaks.begin( ) ; Iter != Peaks.end( ) ; Iter++ ) cout<<Form( " %d", int((*Iter).posi) );
cout<<Form(" )\n");
}
// only the peak with largest area used as pivot
sort( Peaks.begin( ), Peaks.end( ), largerAmplitude() );
int npmax = min( min(Peaks.size(), size_t(4)), Energies.size() ); // limit the number of peaks to min(np,ne,4)
for(int np = 0; np < npmax; np++)
Peaks[np].good = true;
for(size_t np = npmax; np < Peaks.size(); np++)
Peaks[np].good = false;
// (re)order the peaks
sort( Peaks.begin( ), Peaks.end( ), smallerPosi() );
if(Verbosity> 2) {
cout<<Form("# Sorted Peak Positions = (");
for ( Iter = Peaks.begin( ) ; Iter != Peaks.end( ) ; Iter++ ) cout<<Form( " %d", int((*Iter).posi) );
cout<<Form(" )\n");
}
// Connect every (large) peak with every energy and check how many other (peak,energy) pairs agree with this slope
// The combination with the largest number of agreeing pairs provides the reference gain for the final average
// Equal number of matches are ordered by chi2
int np_max = 0; // number of matches
int np_ind = 0; // reference peak
int np_ref = 0; // reference energy
double np_chi = 1.e40; // to distinguish among equal matchers
double dpos2max = specFWHMdef*specFWHMdef;
for(size_t np = 0; np < Peaks.size(); np++) {
if(!Peaks[np].good)
continue;
int ec_max = 0; // number of matches for peak np
int ec_ref = 0; // reference energy
double ec_chi = 1.e30; // its chi2
for(size_t ne = 0; ne < Energies.size(); ne++) { // connecting peak np with energy ne
double lgain = (Peaks[np].posi)/Energies[ne]; // gives this gain
int match = 0;
double chi2 = 0;
for(size_t nee = 0; nee < Energies.size(); nee++) { // count the number of matches for this combination
double epos = Energies[nee]*lgain; // expected position
for(size_t ip = 0; ip < Peaks.size(); ip++) {
double dpos = Peaks[ip].posi-epos;
double dpos2 = dpos*dpos;
if(dpos2 < dpos2max) {
match++;
chi2 += dpos2;
break;
// Connect every (large) peak with every energy and check how many other (peak,energy) pairs agree with this slope
// The combination with the largest number of agreeing pairs provides the reference gain for the final average
// Equal number of matches are ordered by chi2
int np_max = 0; // number of matches
int np_ind = 0; // reference peak
int np_ref = 0; // reference energy
double np_chi = 1.e40; // to distinguish among equal matchers
double dpos2max = specFWHMdef*specFWHMdef;
for(size_t np = 0; np < Peaks.size(); np++) {
if(!Peaks[np].good)
continue;
int ec_max = 0; // number of matches for peak np
int ec_ref = 0; // reference energy
double ec_chi = 1.e30; // its chi2
for(size_t ne = 0; ne < Energies.size(); ne++) { // connecting peak np with energy ne
double lgain = (Peaks[np].posi)/Energies[ne]; // gives this gain
int match = 0;
double chi2 = 0;
for(size_t nee = 0; nee < Energies.size(); nee++) { // count the number of matches for this combination
double epos = Energies[nee]*lgain; // expected position
for(size_t ip = 0; ip < Peaks.size(); ip++) {
double dpos = Peaks[ip].posi-epos;
double dpos2 = dpos*dpos;
if(dpos2 < dpos2max) {
match++;
chi2 += dpos2;
break;
}
}
}
}
if( (match > ec_max) || // more matches
if( (match > ec_max) || // more matches
((match == ec_max) && (chi2 < ec_chi)) ) { // or better chi2
ec_max = match; // record the best combination so far
ec_ref = ne;
ec_chi = chi2;
ec_max = match; // record the best combination so far
ec_ref = ne;
ec_chi = chi2;
}
}
}
if( (ec_max > np_max) || // more matches
if( (ec_max > np_max) || // more matches
((ec_max == np_max) && (ec_chi < np_chi)) ) { // or better chi2
np_max = ec_max; // record the best combination so far
np_ind = np;
np_ref = ec_ref;
np_chi = ec_chi;
np_max = ec_max; // record the best combination so far
np_ind = np;
np_ref = ec_ref;
np_chi = ec_chi;
}
//if(np_max == Energies.size())
// break; // cannot do better than this ?
}
//if(np_max == Energies.size())
// break; // cannot do better than this ?
}
double bestSlope = Energies[np_ref]/(Peaks[np_ind].posi);
if(Verbosity > 2) cout<<Form("# Best-match slope %g [p=%d e=%d] with %d values\n", bestSlope, np_ind, np_ref, np_max);
bestSlope = Energies[np_ref]/(Peaks[np_ind].posi);
if(Verbosity > 2) cout<<Form("# Best-match slope %g [p=%d e=%d] with %d values\n", bestSlope, np_ind, np_ref, np_max);
// find the good peaks
for(size_t np = 0; np < Peaks.size(); np++) {
Peaks[np].erefindex = -1;
Peaks[np].eref = 0.;
Peaks[np].good = false;
}
Float_t Xmin = 0;
Float_t Xmax = 0.;
int match = 0;
for(size_t ne = 0; ne < Energies.size(); ne++) {
double epos = Energies[ne]/bestSlope;
// find the good peaks
for(size_t np = 0; np < Peaks.size(); np++) {
if(abs(Peaks[np].posi-epos) < specFWHMdef && !Peaks[np].good) {
match++;
Peaks[np].erefindex = ne;
Peaks[np].eref = Energies[ne];
Peaks[np].good = true;
// Peaks[np].errposi = pow(Peaks[np].fwhm/s2fwhm, 2.)/Peaks[np].area;
if(Xmin==0. || Peaks[np].posi < Xmin) Xmin = Peaks[np].posi;
if(Peaks[np].posi > Xmax) Xmax = Peaks[np].posi;
break;
Peaks[np].erefindex = -1;
Peaks[np].eref = 0.;
Peaks[np].good = false;
}
int match = 0;
for(size_t ne = 0; ne < Energies.size(); ne++) {
double epos = Energies[ne]/bestSlope;
for(size_t np = 0; np < Peaks.size(); np++) {
if(abs(Peaks[np].posi-epos) < specFWHMdef && !Peaks[np].good) {
match++;
Peaks[np].erefindex = ne;
Peaks[np].eref = Energies[ne];
Peaks[np].good = true;
// Peaks[np].errposi = pow(Peaks[np].fwhm/s2fwhm, 2.)/Peaks[np].area;
if(Xmin==0. || Peaks[np].posi < Xmin) Xmin = Peaks[np].posi;
if(Peaks[np].posi > Xmax) Xmax = Peaks[np].posi;
break;
}
}
}
Xmin-=100;
Xmax+=100;
}
Xmin-=100;
Xmax+=100;
if(fCalibFunction) delete fCalibFunction;
if(fCalibGraph) delete fCalibGraph;
......@@ -1186,7 +1187,6 @@ Int_t GwRecalEnergy::EROOTCalibration()
fCalibFunction->SetParameter(1,bestSlope);
for(int i=2 ; i<=fCalibOrder ; i++) fCalibFunction->SetParameter(i,0);
fCalibGraph = new TGraphErrors;
fCalibGraph->SetName("CalibrationGraph");
fCalibGraph->SetMarkerColor(kRed);
......@@ -1205,7 +1205,7 @@ Int_t GwRecalEnergy::EROOTCalibration()
for(size_t np = 0; np < Peaks.size(); np++) {
if(Peaks[np].good) {
fCalibGraph->SetPoint(fCalibGraph->GetN(), Peaks[np].posi, Energies[Peaks[np].erefindex]);
// fCalibGraph->SetPointError(fCalibGraph->GetN()-1, Peaks[np].errposi, 0.);
// fCalibGraph->SetPointError(fCalibGraph->GetN()-1, Peaks[np].errposi, 0.);
}
}
......@@ -1240,7 +1240,7 @@ Int_t GwRecalEnergy::EROOTCalibration()
double fwhm_e = fCalibFunction->Eval(Peaks[np].fwhm)-fCalibFunction->GetParameter(0);
double res = ecal-Energies[Peaks[np].erefindex];
fResidueGraph->SetPoint(fResidueGraph->GetN(),Energies[Peaks[np].erefindex],res);
// fResidueGraph->SetPointError(fResidueGraph->GetN()-1,0.,r->GetConfidenceIntervals()[fResidueGraph->GetN()-1]);
// fResidueGraph->SetPointError(fResidueGraph->GetN()-1,0.,r->GetConfidenceIntervals()[fResidueGraph->GetN()-1]);
if(Verbosity > 1) cout<<"#2" << setw(12) << setprecision(1) << Peaks[np].area << setw(12) << setprecision(2) << Peaks[np].posi << setw(12) << setprecision(3) << Peaks[np].fwhm << setw(12) << setprecision(3) << fwhm_e << setw(13) << setprecision(3) << ecal << setw(15) << setprecision(3) << res <<endl;
ngood++;
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment