Skip to content

Commit 3df9956

Browse files
cmsbuildcmsbuild
authored
Merge pull request #48555 from bsunanda/Run3-alca257
Run3-alca257 Change FitPlot and PlotPropertiies to chck differences between EGamma and AlCaRaw
2 parents 6d48a7b + a2780ef commit 3df9956

File tree

2 files changed

+433
-22
lines changed

2 files changed

+433
-22
lines changed

‎Calibration/HcalCalibAlgos/macros/CalibFitPlots.C‎

Lines changed: 112 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -568,6 +568,94 @@ results fitOneGauss(TH1D* hist, bool fitTwice, bool debug) {
568568
return results(value, error, width, werror);
569569
}
570570

571+
double DoubleSidedCrystalballFunction(double* x, double* par) {
572+
double alpha_l = par[0];
573+
double alpha_h = par[1];
574+
double n_l = par[2];
575+
double n_h = par[3];
576+
double mean = par[4];
577+
double sigma = par[5];
578+
double N = par[6];
579+
float t = (x[0] - mean) / sigma;
580+
double result;
581+
double fact1TLessMinosAlphaL = alpha_l / n_l;
582+
double fact2TLessMinosAlphaL = (n_l / alpha_l) - alpha_l - t;
583+
double fact1THihgerAlphaH = alpha_h / n_h;
584+
// double fact2THigherAlphaH = (n_h/alpha_h) - alpha_h + t;
585+
586+
if (-alpha_l <= t && alpha_h >= t) {
587+
result = exp(-0.5 * t * t);
588+
} else if (t < -alpha_l) {
589+
result = exp(-0.5 * alpha_l * alpha_l) * pow(fact1TLessMinosAlphaL * fact2TLessMinosAlphaL, -n_l);
590+
} else if (t > alpha_h) {
591+
result = exp(-0.5 * alpha_l * alpha_l) * pow(fact1THihgerAlphaH * fact1THihgerAlphaH, -n_h);
592+
}
593+
return N * result;
594+
}
595+
596+
results fitDoubleSidedCrystalball(TH1D* hist, bool /* fitTwice */, bool debug) {
597+
double rms;
598+
std::pair<double, double> mrms = GetMean(hist, 0.2, 2.0, rms);
599+
double mean = mrms.first;
600+
double LowEdge = ((mean - fitrangeFactor * rms) < 0.5) ? 0.5 : (mean - fitrangeFactor * rms);
601+
double diff = mean - LowEdge;
602+
double HighEdge = (diff > fitrangeFactor1 * rms) ? (mean + fitrangeFactor1 * rms) : (mean + diff);
603+
if (debug)
604+
std::cout << hist->GetName() << " Mean " << mean << " RMS " << rms << " Range " << LowEdge << ":" << HighEdge
605+
<< "\n";
606+
std::string option = (hist->GetEntries() > 100) ? "QRS" : "QRWLS";
607+
TObject* ob = gROOT->FindObject("g1");
608+
if (ob)
609+
ob->Delete();
610+
TF1* g1 = new TF1("g1", "gaus", LowEdge, HighEdge);
611+
g1->SetLineColor(kGreen);
612+
TFitResultPtr Fit1 = hist->Fit(g1, option.c_str(), "");
613+
double value = Fit1->Value(1);
614+
double error = Fit1->FitResult::Error(1);
615+
double width = Fit1->Value(2);
616+
double werror = Fit1->FitResult::Error(2);
617+
618+
LowEdge = Fit1->Value(1) - fitrangeFactor * Fit1->Value(2);
619+
if (LowEdge < 0.5)
620+
LowEdge = 0.5;
621+
diff = Fit1->Value(1) - LowEdge;
622+
HighEdge = (diff > fitrangeFactor1 * rms) ? (Fit1->Value(1) + fitrangeFactor1 * Fit1->Value(2))
623+
: (Fit1->Value(1) + 1.5 * diff);
624+
if (HighEdge > 5.0)
625+
HighEdge = 5.0;
626+
if (debug)
627+
std::cout << " Range for second Fit " << LowEdge << ":" << HighEdge << std::endl;
628+
TObject* ob2 = gROOT->FindObject("fitDSCB");
629+
if (ob2 != nullptr)
630+
ob2->Delete();
631+
TF1* fitDSCB = new TF1("fitDSCB", DoubleSidedCrystalballFunction, LowEdge, HighEdge, 7);
632+
fitDSCB->SetParameters(1, 2, 2, 1, value, width, hist->Integral(LowEdge, HighEdge));
633+
fitDSCB->SetParNames("alpha_{low}", "alpha_{high}", "n_{low}", "n_{high}", "mean", "sigma", "Norm");
634+
TFitResultPtr Fit = hist->Fit(fitDSCB, option.c_str(), "", LowEdge, HighEdge);
635+
value = Fit->Value(4);
636+
error = Fit->FitResult::Error(4);
637+
width = Fit->Value(5);
638+
werror = Fit->FitResult::Error(5);
639+
640+
std::pair<double, double> meaner = GetMean(hist, 0.2, 2.0, rms);
641+
std::pair<double, double> rmserr = GetWidth(hist, 0.2, 2.0);
642+
if (debug) {
643+
std::cout << "FitDSCB " << value << ":" << error << ":" << hist->GetMeanError() << " Mean " << meaner.first << ":"
644+
<< meaner.second << " Width " << rmserr.first << ":" << rmserr.second;
645+
}
646+
double minvalue(0.30);
647+
if (value < minvalue || value > 2.0 || error > 0.5) {
648+
value = meaner.first;
649+
error = meaner.second;
650+
width = rmserr.first;
651+
werror = rmserr.second;
652+
}
653+
if (debug) {
654+
std::cout << " Final " << value << "+-" << error << ":" << width << "+-" << werror << std::endl;
655+
}
656+
return results(value, error, width, werror);
657+
}
658+
571659
void readCorrFactors(char* infile,
572660
double scale,
573661
std::map<int, cfactors>& cfacs,
@@ -640,7 +728,7 @@ void FitHistStandard(std::string infile,
640728
std::string xname[4] = {"i#eta", "i#eta", "d_{L1}", "# Vertex"};
641729
int numb[4] = {10, 8, 8, 5};
642730

643-
if (type == 0) {
731+
if ((type % 10) == 0) {
644732
numb[0] = 8;
645733
for (int i = 0; i < 9; ++i)
646734
xbins[i] = xbin0[i];
@@ -683,7 +771,8 @@ void FitHistStandard(std::string infile,
683771
}
684772
if (hist->GetEntries() > 4) {
685773
bool flag = (j == 0) ? true : false;
686-
results meaner = fitOneGauss(hist, flag, debug);
774+
results meaner = (((type / 10) % 10) == 0) ? fitOneGauss(hist, flag, debug)
775+
: fitDoubleSidedCrystalball(hist, flag, debug);
687776
value = meaner.mean;
688777
error = meaner.errmean;
689778
width = meaner.width;
@@ -756,11 +845,11 @@ void FitHistExtended(const char* infile,
756845
double xbins[99];
757846
double xbin[23] = {-23.0, -21.0, -19.0, -17.0, -15.0, -13.0, -11.0, -9.0, -7.0, -5.0, -3.0, 0.0,
758847
3.0, 5.0, 7.0, 9.0, 11.0, 13.0, 15.0, 17.0, 19.0, 21.0, 23.0};
759-
if (type == 2) {
848+
if ((type % 10) == 2) {
760849
numb = 22;
761850
for (int k = 0; k <= numb; ++k)
762851
xbins[k] = xbin[k];
763-
} else if (type == 1) {
852+
} else if ((type % 10) == 1) {
764853
numb = 1;
765854
xbins[0] = -25;
766855
xbins[1] = 25;
@@ -802,7 +891,8 @@ void FitHistExtended(const char* infile,
802891
}
803892
if (hist0->GetEntries() > 10) {
804893
double rms;
805-
results meaner0 = fitOneGauss(hist0, true, debug);
894+
results meaner0 =
895+
(((type / 10) % 10) == 0) ? fitOneGauss(hist0, true, debug) : fitDoubleSidedCrystalball(hist0, true, debug);
806896
std::pair<double, double> meaner1 = GetMean(hist0, 0.2, 2.0, rms);
807897
std::pair<double, double> meaner2 = GetWidth(hist0, 0.2, 2.0);
808898
if (debug) {
@@ -846,13 +936,15 @@ void FitHistExtended(const char* infile,
846936
TH1D* hist2 = (TH1D*)hist1->Clone(name);
847937
fitOneGauss(hist2, true, debug);
848938
hists.push_back(hist2);
849-
results meaner = fitOneGauss(hist, true, debug);
939+
results meaner = (((type / 10) % 10) == 0) ? fitOneGauss(hist, true, debug)
940+
: fitDoubleSidedCrystalball(hist, true, debug);
850941
value = meaner.mean;
851942
error = meaner.errmean;
852943
width = meaner.width;
853944
werror = meaner.errwidth;
854945
} else {
855-
results meaner = fitOneGauss(hist, true, debug);
946+
results meaner = (((type / 10) % 10) == 0) ? fitOneGauss(hist, true, debug)
947+
: fitDoubleSidedCrystalball(hist, true, debug);
856948
value = meaner.mean;
857949
error = meaner.errmean;
858950
width = meaner.width;
@@ -923,7 +1015,8 @@ void FitHistExtended(const char* infile,
9231015
if (total > 4) {
9241016
sprintf(name, "%sOne", hist1->GetName());
9251017
TH1D* hist2 = (TH1D*)hist1->Clone(name);
926-
results meanerr = fitOneGauss(hist2, true, debug);
1018+
results meanerr = (((type / 10) % 10) == 0) ? fitOneGauss(hist2, true, debug)
1019+
: fitDoubleSidedCrystalball(hist2, true, debug);
9271020
value = meanerr.mean;
9281021
error = meanerr.errmean;
9291022
width = meanerr.width;
@@ -940,7 +1033,8 @@ void FitHistExtended(const char* infile,
9401033
hists.push_back(hist3);
9411034
}
9421035
// results meaner0 = fitTwoGauss(hist, debug);
943-
results meaner0 = fitOneGauss(hist, true, debug);
1036+
results meaner0 = (((type / 10) % 10) == 0) ? fitOneGauss(hist, true, debug)
1037+
: fitDoubleSidedCrystalball(hist, true, debug);
9441038
value = meaner0.mean;
9451039
error = meaner0.errmean;
9461040
double rms;
@@ -984,6 +1078,7 @@ void FitHistExtended2(const char* infile,
9841078
const char* outfile,
9851079
std::string prefix,
9861080
int numb = 50,
1081+
int type = 0,
9871082
bool append = true,
9881083
int iname = 3,
9891084
bool debug = false) {
@@ -1057,13 +1152,15 @@ void FitHistExtended2(const char* infile,
10571152
TH1D* hist2 = (TH1D*)hist1->Clone(name);
10581153
fitOneGauss(hist2, true, debug);
10591154
hists.push_back(hist2);
1060-
results meaner = fitOneGauss(hist, true, debug);
1155+
results meaner = (((type / 10) % 10) == 0) ? fitOneGauss(hist, true, debug)
1156+
: fitDoubleSidedCrystalball(hist, true, debug);
10611157
value = meaner.mean;
10621158
error = meaner.errmean;
10631159
width = meaner.width;
10641160
werror = meaner.errwidth;
10651161
} else {
1066-
results meaner = fitOneGauss(hist, true, debug);
1162+
results meaner = (((type / 10) % 10) == 0) ? fitOneGauss(hist, true, debug)
1163+
: fitDoubleSidedCrystalball(hist, true, debug);
10671164
value = meaner.mean;
10681165
error = meaner.errmean;
10691166
width = meaner.width;
@@ -1139,7 +1236,8 @@ void FitHistExtended2(const char* infile,
11391236
if (total > 4) {
11401237
sprintf(name, "%sOne", hist1->GetName());
11411238
TH1D* hist2 = (TH1D*)hist1->Clone(name);
1142-
results meanerr = fitOneGauss(hist2, true, debug);
1239+
results meanerr = (((type / 10) % 10) == 0) ? fitOneGauss(hist2, true, debug)
1240+
: fitDoubleSidedCrystalball(hist2, true, debug);
11431241
value = meanerr.mean;
11441242
error = meanerr.errmean;
11451243
width = meanerr.width;
@@ -1149,7 +1247,8 @@ void FitHistExtended2(const char* infile,
11491247
std::cout << hist2->GetName() << " MPV " << value << " +- " << error << " Width " << width << " +- "
11501248
<< werror << " W/M " << wbyv << " +- " << wverr << std::endl;
11511249
hists.push_back(hist2);
1152-
results meaner0 = fitOneGauss(hist, true, debug);
1250+
results meaner0 = (((type / 10) % 10) == 0) ? fitOneGauss(hist, true, debug)
1251+
: fitDoubleSidedCrystalball(hist, true, debug);
11531252
value = meaner0.mean;
11541253
error = meaner0.errmean;
11551254
double rms;

0 commit comments

Comments
 (0)